Hi fellow PETSc users!
I am attempting to use a DMCOMPOSITE alongside TS and have run into some
trouble. I'm attaching a MWE demonstrating the problem. The goal is to
combine a DMDA3d for spatial data and a DMREDUNDANT for non-spatial,
time-dependent fields. In the attached example, this additional field is
temperature.
The DMDA data is currently just temperature-dependent, and the temperature
is supposed to increase linearly. Unfortunately, no matter how I integrate
(explicitly, using Euler or RK, or implicitly, using backward Euler or
BDF), I get the wrong final temperature. The integration proceeds for 100
timesteps and then stops (verified by -ts_monitor). At a heating rate of 1,
and an initial temperature of 150, I should get a final temperature of 250
very easily. However, I get something closer to 200 (not exact, which makes
me think this isn't a simple missing-a-factor-of-2 error).
I'm currently using TSSetDM with the composite DM to let PETSc compute the
Jacobian. Having checked all the inputs and outputs as well as I can, this
is the only step that seems like it may be causing a problem; yet even
explicit integration doesn't work, and that shouldn't need the Jacobian at
all. I'm at a loss.
Run using:
./mwe -implicit yes -ts_type beuler -ts_monitor OR ./mwe -implicit no
-ts_type euler -ts_monitor
Thanks in advance for any help,
Ellen Price
#include
#include
#include
#include
#include
#define ZERO(0.)
#define ONE (1.)
#define AUXTEMP (0)
#define TEMP0 (150.)
#define HEATING (1.)
#define NU (0.1)
#define TFINAL (100.)
#define DELTA (1.)
#define unused(x) ((void) (x))
PetscErrorCode rhsfunc_petsc(TS ts, PetscReal t, Vec y, Vec ydot, void *ptr);
PetscErrorCode ifunc_petsc(TS ts, PetscReal t, Vec y, Vec ydot, Vec f, void *ptr);
void func(double t, Vec y, Vec aux, Vec ydot, Vec auxdot, DM da);
PetscErrorCode rhsfunc_petsc(TS ts, PetscReal t, Vec y, Vec ydot, void *ptr)
{
DM pack, da;
Vec yda, yred, ydotda, ydotred;
unused(ptr);
TSGetDM(ts, &pack);
DMCompositeGetEntries(pack, &da, NULL);
DMCompositeGetLocalVectors(pack, &yda, &yred);
DMCompositeGetLocalVectors(pack, &ydotda, &ydotred);
DMCompositeScatter(pack, y, yda, yred);
VecSet(ydotda, ZERO);
VecSet(ydotred, ZERO);
func(t, yda, yred, ydotda, ydotred, da);
DMCompositeGather(pack, INSERT_VALUES, ydot, ydotda, ydotred);
DMCompositeRestoreLocalVectors(pack, &ydotda, &ydotred);
DMCompositeRestoreLocalVectors(pack, &yda, &yred);
VecView(ydot, PETSC_VIEWER_STDOUT_WORLD);
return 0;
}
PetscErrorCode ifunc_petsc(TS ts, PetscReal t, Vec y, Vec ydot, Vec f, void *ptr)
{
DM pack, da;
Vec yda, yred, ydotda, ydotred;
unused(ptr);
TSGetDM(ts, &pack);
DMCompositeGetEntries(pack, &da, NULL);
VecCopy(ydot, f);
DMCompositeGetLocalVectors(pack, &yda, &yred);
DMCompositeGetLocalVectors(pack, &ydotda, &ydotred);
DMCompositeScatter(pack, y, yda, yred);
VecSet(ydotda, ZERO);
VecSet(ydotred, ZERO);
func(t, yda, yred, ydotda, ydotred, da);
VecScale(ydotda, -ONE);
VecScale(ydotred, -ONE);
DMCompositeGather(pack, ADD_VALUES, f, ydotda, ydotred);
DMCompositeRestoreLocalVectors(pack, &ydotda, &ydotred);
DMCompositeRestoreLocalVectors(pack, &yda, &yred);
return 0;
}
void func(double t, Vec y, Vec aux, Vec ydot, Vec auxdot, DM da)
{
int rank;
PetscScalar ydotarr, *auxdotarr;
const PetscScalar yarr, *auxarr;
PetscInt i, j, k, is, js, ks, n, m, p, ys, ye;
unused(t);
MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
VecGetArrayRead(aux, &auxarr);
DMDAGetCorners(da, &is, &js, &ks, &n, &m, &p);
DMDAVecGetArrayDOFRead(da, y, &yarr);
DMDAVecGetArrayDOF(da, ydot, &ydotarr);
for (k = ks; k < ks + p; ++k)
{
for (j = js; j < js + m; ++j)
{
for (i = is; i < is + n; ++i)
{
double temp = auxarr[AUXTEMP];
double energy = 100.;
double R = NU * exp(-energy / temp);
ydotarr[k][j][i][0] += R;
ydotarr[k][j][i][1] -= R;
}
}
}
DMDAVecRestoreArrayDOFRead(da, y, &yarr);
DMDAVecRestoreArrayDOF(da, ydot, &ydotarr);
VecGetOwnershipRange(auxdot, &ys, &ye);
VecGetArray(auxdot, &auxdotarr);
if (rank == 0)
{
for (i = ys; i < ye; ++i)
{
switch (i)
{
case AUXTEMP:
auxdotarr[i-ys] = HEATING;
break;
}
}
}
VecRestoreArray(auxdot, &auxdotarr);
VecRestoreArrayRead(aux, &auxarr);
}
int main(int argc, char *argv[])
{
DM da, red, pack;
Vec y, yda, yred;
PetscScalar yarr, *auxarr;
PetscInt i, j, k, is, js, ks, n, m, p, as, ae;
double t = ZERO, abstol = 1.e-6, reltol = 1.e-3;
int rank, tt
