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 <stdlib.h>
#include <petscdmda.h>
#include <petscdmredundant.h>
#include <petscdmcomposite.h>
#include <petscts.h>
#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 = 0, ndof = 2, naux = 1, nx = 3, ny = 3, nz = 3;
PetscBool implicit = PETSC_TRUE;
TS ts;
TSAdapt adapt;
TSConvergedReason reason;
PetscInitialize(&argc, &argv, NULL, NULL);
PetscOptionsGetBool(NULL, NULL, "-implicit", &implicit, NULL);
MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
DMDACreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,
DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, nx, ny, nz, PETSC_DECIDE,
PETSC_DECIDE, PETSC_DECIDE, ndof, 2, NULL, NULL, NULL, &da);
DMSetFromOptions(da);
DMSetUp(da);
DMRedundantCreate(PETSC_COMM_WORLD, 0, naux, &red);
DMSetFromOptions(red);
DMSetUp(red);
DMCompositeCreate(PETSC_COMM_WORLD, &pack);
DMSetFromOptions(pack);
DMCompositeAddDM(pack, da);
DMCompositeAddDM(pack, red);
DMSetUp(pack);
DMCreateGlobalVector(pack, &y);
DMCompositeGetLocalVectors(pack, &yda, &yred);
VecSet(yda, ZERO);
DMDAGetCorners(da, &is, &js, &ks, &n, &m, &p);
DMDAVecGetArrayDOF(da, yda, &yarr);
for (k = ks; k < ks + p; ++k)
{
for (j = js; j < js + m; ++j)
{
for (i = is; i < is + n; ++i)
{
yarr[k][j][i][0] = ZERO;
yarr[k][j][i][1] = ONE;
}
}
}
DMDAVecRestoreArrayDOF(da, yda, &yarr);
VecGetOwnershipRange(yred, &as, &ae);
VecGetArray(yred, &auxarr);
if (rank == 0)
{
for (i = as; i < ae; ++i)
{
switch (i)
{
case AUXTEMP:
auxarr[i-as] = TEMP0;
break;
}
}
}
VecRestoreArray(yred, &auxarr);
DMCompositeGather(pack, INSERT_VALUES, y, yda, yred);
DMCompositeRestoreLocalVectors(pack, &yda, &yred);
TSCreate(PETSC_COMM_WORLD, &ts);
TSSetDM(ts, pack);
TSSetTolerances(ts, abstol, NULL, reltol, NULL);
TSSetExactFinalTime(ts, TS_EXACTFINALTIME_INTERPOLATE);
TSGetAdapt(ts, &adapt);
TSAdaptSetType(adapt, TSADAPTBASIC);
TSAdaptSetStepLimits(adapt, ZERO, DELTA);
TSSetFromOptions(ts);
if (implicit == PETSC_TRUE)
{
TSSetIFunction(ts, NULL, ifunc_petsc, NULL);
}
else
{
TSSetRHSFunction(ts, NULL, rhsfunc_petsc, NULL);
}
do
{
TSSetMaxTime(ts, t + DELTA);
TSSolve(ts, y);
TSGetConvergedReason(ts, &reason);
if (reason < 0) break;
TSGetSolveTime(ts, &t);
tt++;
}
while (t < TFINAL);
VecView(y, PETSC_VIEWER_STDOUT_WORLD);
TSDestroy(&ts);
VecDestroy(&yred);
VecDestroy(&yda);
VecDestroy(&y);
DMDestroy(&da);
DMDestroy(&red);
DMDestroy(&pack);
PetscFinalize();
return 0;
}
/* vim: set ft=c.doxygen: */
