I checked out branch barry/fix-huge-flaw-in-ts
<https://bitbucket.org/petsc/petsc/branch/barry/fix-huge-flaw-in-ts> (see
pull request #655) and reconfigured and rebuilt.
No, the ts/.../ex2.c example is not fixed. It gives same error:
~/petsc/src/ts/examples/tutorials[barry/fix-huge-flaw-in-ts*]$ ./ex2 -ts_dt
0.001 -ts_final_time 0.001 -snes_monitor -snes_mf_operator
0 SNES Function norm 1.059316422854e+01
[0]PETSC ERROR: --------------------- Error Message
--------------------------------------------------------------
[0]PETSC ERROR: No support for this operation for this object type
[0]PETSC ERROR: Mat type mffd
[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for
trouble shooting.
[0]PETSC ERROR: Petsc Development GIT revision: v3.7.5-3346-g41a1d4d GIT
Date: 2017-03-24 14:39:28 -0500
[0]PETSC ERROR: ./ex2 on a linux-c-dbg named bueler-leopard by ed Sun Apr
2 18:31:55 2017
[0]PETSC ERROR: Configure options --download-mpich --with-debugging=1
[0]PETSC ERROR: #1 MatZeroEntries() line 5471 in
/home/ed/petsc/src/mat/interface/matrix.c
[0]PETSC ERROR: #2 TSComputeIJacobian() line 942 in
/home/ed/petsc/src/ts/interface/ts.c
[0]PETSC ERROR: #3 SNESTSFormJacobian_Theta() line 515 in
/home/ed/petsc/src/ts/impls/implicit/theta/theta.c
[0]PETSC ERROR: #4 SNESTSFormJacobian() line 5055 in
/home/ed/petsc/src/ts/interface/ts.c
[0]PETSC ERROR: #5 SNESComputeJacobian() line 2276 in
/home/ed/petsc/src/snes/interface/snes.c
[0]PETSC ERROR: #6 SNESSolve_NEWTONLS() line 222 in
/home/ed/petsc/src/snes/impls/ls/ls.c
[0]PETSC ERROR: #7 SNESSolve() line 3967 in
/home/ed/petsc/src/snes/interface/snes.c
[0]PETSC ERROR: #8 TS_SNESSolve() line 171 in
/home/ed/petsc/src/ts/impls/implicit/theta/theta.c
[0]PETSC ERROR: #9 TSStep_Theta() line 211 in
/home/ed/petsc/src/ts/impls/implicit/theta/theta.c
[0]PETSC ERROR: #10 TSStep() line 3820 in
/home/ed/petsc/src/ts/interface/ts.c
[0]PETSC ERROR: #11 TSSolve() line 4065 in
/home/ed/petsc/src/ts/interface/ts.c
[0]PETSC ERROR: #12 main() line 194 in
/home/ed/petsc/src/ts/examples/tutorials/ex2.c
[0]PETSC ERROR: PETSc Option Table entries:
[0]PETSC ERROR: -snes_mf_operator
[0]PETSC ERROR: -snes_monitor
[0]PETSC ERROR: -ts_dt 0.001
[0]PETSC ERROR: -ts_final_time 0.001
[0]PETSC ERROR: ----------------End of Error Message -------send entire
error message to [email protected]
application called MPI_Abort(MPI_COMM_WORLD, 56) - process 0
[unset]: aborting job:
application called MPI_Abort(MPI_COMM_WORLD, 56) - process 0
Curiously, my own code (attached), which differs from ex2.c in that it uses
DMDA, *does* change behavior to a DIVERGED_NONLINEAR_SOLVE error:
$ ./advect -da_refine 0 -ts_monitor -adv_circlewind -adv_conex 0.3
-adv_coney 0.3 -ts_type beuler -snes_monitor_short -ts_final_time 0.01
-ts_dt 0.01 -snes_rtol 1.0e-4 -adv_firstorder -snes_mf_operator
solving on 5 x 5 grid with dx=0.2 x dy=0.2 cells, t0=0.,
and initial step dt=0.01 ...
0 TS dt 0.01 time 0.
0 SNES Function norm 2.1277
[0]PETSC ERROR: --------------------- Error Message
--------------------------------------------------------------
[0]PETSC ERROR:
[0]PETSC ERROR: TSStep has failed due to DIVERGED_NONLINEAR_SOLVE, increase
-ts_max_snes_failures or make negative to attempt recovery
[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for
trouble shooting.
[0]PETSC ERROR: Petsc Development GIT revision: v3.7.5-3346-g41a1d4d GIT
Date: 2017-03-24 14:39:28 -0500
[0]PETSC ERROR: ./advect on a linux-c-dbg named bueler-leopard by ed Sun
Apr 2 18:34:18 2017
[0]PETSC ERROR: Configure options --download-mpich --with-debugging=1
[0]PETSC ERROR: #1 TSStep() line 3829 in
/home/ed/petsc/src/ts/interface/ts.c
[0]PETSC ERROR: #2 TSSolve() line 4065 in
/home/ed/petsc/src/ts/interface/ts.c
[0]PETSC ERROR: #3 main() line 314 in /home/ed/repos/p4pdes/c/ch9/advect.c
[0]PETSC ERROR: PETSc Option Table entries:
[0]PETSC ERROR: -adv_circlewind
[0]PETSC ERROR: -adv_conex 0.3
[0]PETSC ERROR: -adv_coney 0.3
[0]PETSC ERROR: -adv_firstorder
[0]PETSC ERROR: -da_refine 0
[0]PETSC ERROR: -snes_mf_operator
[0]PETSC ERROR: -snes_monitor_short
[0]PETSC ERROR: -snes_rtol 1.0e-4
[0]PETSC ERROR: -ts_dt 0.01
[0]PETSC ERROR: -ts_final_time 0.01
[0]PETSC ERROR: -ts_monitor
[0]PETSC ERROR: -ts_type beuler
[0]PETSC ERROR: ----------------End of Error Message -------send entire
error message to [email protected]
application called MPI_Abort(MPI_COMM_WORLD, 91) - process 0
[unset]: aborting job:
application called MPI_Abort(MPI_COMM_WORLD, 91) - process 0
Previously it gave a "No support for this operation for this object type
... Mat type mffd" error like ex2.c. Progress?
Note that my code works fine (i.e. no errors) with any of: the analytical
Jacobian or -snes_fd or -snes_fd_color or -snes_mf.
Ed
On Sun, Apr 2, 2017 at 6:23 PM, Jed Brown <[email protected]> wrote:
> The issue is that we need to create
>
> a*I - Jrhs
>
> and this is currently done by creating a*I first when we have separate
> matrices for the left and right hand sides. There is code to just scale
> and shift Jrhs when there is no IJacobian, but the creation logic got
> messed up at some point (or at least for some TS configurations that are
> common). We were discussing this recently in a performance context and
> this branch is supposed to fix that logic. Does this branch work for
> you?
>
> https://bitbucket.org/petsc/petsc/pull-requests/655/fix-
> flaw-with-tssetrhsjacobian-and-no/diff
>
> Ed Bueler <[email protected]> writes:
>
> > Dear PETSc --
> >
> > I have a TS-using and DMDA-using code in which I want to set a
> RHSJacobian
> > which is only approximate. (The Jacobian uses first-order upwinding MOL
> > while the RHSFunction uses a flux-limited MOL.) While it works with the
> > analytical Jacobian, and -snes_fd, and -snes_fd_color, and -snes_mf, I
> get
> > a "No support ..." message for -snes_mf_operator.
> >
> > It suffices to show the problem with
> >
> > src/ts/examples/tutorials/ex2.c
> >
> > (Note ex2.c does not use DMDA so ..._color is not available for ex2.c.)
> >
> > Error as follows:
> >
> > ~/petsc/src/ts/examples/tutorials[master*]$ ./ex2 -ts_dt 0.001
> > -ts_final_time 0.001 -snes_monitor
> > 0 SNES Function norm 1.059316422854e+01
> > 1 SNES Function norm 1.035505461114e-05
> > 2 SNES Function norm 5.498223366328e-12
> > ~/petsc/src/ts/examples/tutorials[master*]$ ./ex2 -ts_dt 0.001
> > -ts_final_time 0.001 -snes_monitor -snes_fd
> > 0 SNES Function norm 1.059316422854e+01
> > 1 SNES Function norm 1.208245988550e-05
> > 2 SNES Function norm 6.022374930788e-12
> > ~/petsc/src/ts/examples/tutorials[master*]$ ./ex2 -ts_dt 0.001
> > -ts_final_time 0.001 -snes_monitor -snes_mf
> > 0 SNES Function norm 1.059316422854e+01
> > 1 SNES Function norm 6.136984336801e-05
> > 2 SNES Function norm 5.355730806625e-10
> > ~/petsc/src/ts/examples/tutorials[master*]$ ./ex2 -ts_dt 0.001
> > -ts_final_time 0.001 -snes_monitor -snes_mf_operator
> > 0 SNES Function norm 1.059316422854e+01
> > [0]PETSC ERROR: --------------------- Error Message
> > --------------------------------------------------------------
> > [0]PETSC ERROR: No support for this operation for this object type
> > [0]PETSC ERROR: Mat type mffd
> > [0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html
> for
> > trouble shooting.
> > [0]PETSC ERROR: Petsc Development GIT revision: v3.7.5-3426-g0c7851c GIT
> > Date: 2017-04-01 18:40:06 -0600
> > [0]PETSC ERROR: ./ex2 on a linux-c-dbg named ed-lemur by ed Sun Apr 2
> > 17:02:47 2017
> > [0]PETSC ERROR: Configure options --download-mpich --with-debugging=1
> > [0]PETSC ERROR: #1 MatZeroEntries() line 5471 in
> > /home/ed/petsc/src/mat/interface/matrix.c
> > [0]PETSC ERROR: #2 TSComputeIJacobian() line 965 in
> > /home/ed/petsc/src/ts/interface/ts.c
> > [0]PETSC ERROR: #3 SNESTSFormJacobian_Theta() line 515 in
> > /home/ed/petsc/src/ts/impls/implicit/theta/theta.c
> > [0]PETSC ERROR: #4 SNESTSFormJacobian() line 5078 in
> > /home/ed/petsc/src/ts/interface/ts.c
> > [0]PETSC ERROR: #5 SNESComputeJacobian() line 2276 in
> > /home/ed/petsc/src/snes/interface/snes.c
> > [0]PETSC ERROR: #6 SNESSolve_NEWTONLS() line 222 in
> > /home/ed/petsc/src/snes/impls/ls/ls.c
> > [0]PETSC ERROR: #7 SNESSolve() line 3967 in
> > /home/ed/petsc/src/snes/interface/snes.c
> > [0]PETSC ERROR: #8 TS_SNESSolve() line 171 in
> > /home/ed/petsc/src/ts/impls/implicit/theta/theta.c
> > [0]PETSC ERROR: #9 TSStep_Theta() line 211 in
> > /home/ed/petsc/src/ts/impls/implicit/theta/theta.c
> > [0]PETSC ERROR: #10 TSStep() line 3843 in
> > /home/ed/petsc/src/ts/interface/ts.c
> > [0]PETSC ERROR: #11 TSSolve() line 4088 in
> > /home/ed/petsc/src/ts/interface/ts.c
> > [0]PETSC ERROR: #12 main() line 194 in
> > /home/ed/petsc/src/ts/examples/tutorials/ex2.c
> > [0]PETSC ERROR: PETSc Option Table entries:
> > [0]PETSC ERROR: -snes_mf_operator
> > [0]PETSC ERROR: -snes_monitor
> > [0]PETSC ERROR: -ts_dt 0.001
> > [0]PETSC ERROR: -ts_final_time 0.001
> > [0]PETSC ERROR: ----------------End of Error Message -------send entire
> > error message to [email protected]
> > application called MPI_Abort(MPI_COMM_WORLD, 56) - process 0
> > [unset]: aborting job:
> > application called MPI_Abort(MPI_COMM_WORLD, 56) - process 0
> >
> > Is this intended or a bug? If it is intended, what is the issue?
> >
> > I am using current master branch (commit
> > 0c7851c55cba8e40da5083f79ba1ff846acd45b2).
> >
> > Thanks for your help and awesome library!
> >
> > Ed
> >
> >
> >
> >
> >
> > --
> > Ed Bueler
> > Dept of Math and Stat and Geophysical Institute
> > University of Alaska Fairbanks
> > Fairbanks, AK 99775-6660
> > 301C Chapman
>
--
Ed Bueler
Dept of Math and Stat and Geophysical Institute
University of Alaska Fairbanks
Fairbanks, AK 99775-6660
301C Chapman
static char help[] =
"Time-dependent pure-advection equation in 2D using TS. Option prefix -adv_.\n"
"Domain is (0,1) x (0,1). Equation is\n"
" u_t + div(a(x,y) u) = g(x,y,u).\n"
"Boundary conditions are periodic in x and y. Cells are grid-point centered.\n"
"Uses van Leer (1974) flux-limited (non-oscillatory) method-of-lines\n"
"discretization [default], or first-order upwind.\n";
#include <petsc.h>
// try:
// ./advect -da_refine 5 -ts_monitor_solution draw -ts_monitor -ts_rk_type 5dp
// one lap of circular motion, computed in parallel:
// mpiexec -n 4 ./advect -da_refine 5 -adv_circlewind -adv_conex 0.3 -adv_coney 0.3 -ts_final_time 3.1416 -ts_monitor -ts_monitor_solution draw -ksp_rtol 1.0e-10
// implicit:
// mpiexec -n 4 ./advect -ts_monitor_solution draw -ts_monitor -adv_circlewind -ts_final_time 0.5 -adv_conex 0.3 -adv_coney 0.3 -ts_type cn -da_refine 6 -snes_monitor -ts_dt 0.05 -ksp_rtol 1.0e-10
// with -adv_firstorder, -snes_type test suggests Jacobian is correct
// testing Jacobian options (fails with XX=-snes_mf_operator; use -ksp_view_mat ::ascii_matlab etc.):
// ./advect -da_refine 0 -ts_monitor -adv_circlewind -adv_conex 0.3 -adv_coney 0.3 -ts_type beuler -snes_monitor_short -ts_final_time 0.01 -ts_dt 0.01 -snes_rtol 1.0e-4 -adv_firstorder XX
typedef struct {
PetscBool firstorder, // if true, use first-order upwinding
circlewind; // if true, wind is equivalent to rigid rotation
double windx, windy, // x,y components of wind (if not circular)
conex, coney, coner, coneh; // parameters for cone initial cond.
} AdvectCtx;
PetscErrorCode FormInitial(DMDALocalInfo *info, Vec u, AdvectCtx* user) {
PetscErrorCode ierr;
DMDACoor2d **coords;
int i, j;
double x, y, r, **au;
ierr = VecSet(u,0.0); CHKERRQ(ierr); // clear it first
ierr = DMDAGetCoordinateArray(info->da, &coords); CHKERRQ(ierr);
ierr = DMDAVecGetArray(info->da, u, &au); CHKERRQ(ierr);
for (j=info->ys; j<info->ys+info->ym; j++) {
for (i=info->xs; i<info->xs+info->xm; i++) {
x = coords[j][i].x - user->conex;
y = coords[j][i].y - user->coney;
r = PetscSqrtReal(x * x + y * y);
if (r < user->coner) {
au[j][i] = user->coneh * (1.0 - r / user->coner);
}
}
}
ierr = DMDAVecRestoreArray(info->da, u, &au); CHKERRQ(ierr);
ierr = DMDARestoreCoordinateArray(info->da, &coords); CHKERRQ(ierr);
return 0;
}
// velocity a(x,y) = ( a^x(x,y), a^y(x,y) )
static double a_wind(double x, double y, int dir, AdvectCtx* user) {
if (user->circlewind) {
return (dir == 0) ? - 2.0 * (y - 0.5) : 2.0 * (x - 0.5);
} else {
return (dir == 0) ? user->windx : user->windy;
}
}
// source g(x,y,u)
static double g_source(double x, double y, double u, AdvectCtx* user) {
return 0.0;
}
// d g(x,y,u) / d u
static double dg_source(double x, double y, double u, AdvectCtx* user) {
return 0.0;
}
/* the van Leer (1974) limiter is formula (1.11) in section III.1 of
Hundsdorfer & Verwer */
//STARTLIMITER
static double limiter(double th) {
return 0.5 * (th + PetscAbsReal(th)) / (1.0 + PetscAbsReal(th));
}
//ENDLIMITER
/* method-of-lines discretization gives ODE system u' = G(t,u) */
//STARTFUNCTION
PetscErrorCode FormRHSFunctionLocal(DMDALocalInfo *info, double t,
double **au, double **aG, AdvectCtx *user) {
int i, j, l;
const int dir[4] = {0, 1, 0, 1}, // use x (0) or y (1) component
xsh[4] = { 1, 0,-1, 0}, ysh[4] = { 0, 1, 0,-1},
aposi[4] = { 0, 0,-1, 0}, anegi[4] = { 1, 0, 0, 0},
aposj[4] = { 0, 0, 0,-1}, anegj[4] = { 0, 1, 0, 0},
posfari[4] = {-1, 0,-2, 0}, negfari[4] = { 2, 0, 1, 0},
posfarj[4] = { 0,-1, 0,-2}, negfarj[4] = { 0, 2, 0, 1};
double hx, hy, halfx, halfy, x, y, a, u_up, u_dn, u_far, theta,
flux[4];
// fluxes on cell boundaries are traversed in this order (=l): ,-1-,
// cell center at * has coordinates (x,y): 2 * 0
// '-3-'
hx = 1.0 / info->mx; hy = 1.0 / info->my;
halfx = hx / 2.0; halfy = hy / 2.0;
for (j = info->ys; j < info->ys + info->ym; j++) {
y = (j + 0.5) * hy;
for (i = info->xs; i < info->xs + info->xm; i++) {
x = (i + 0.5) * hx;
for (l = 0; l < 4; l++) { // loop over cell boundaries
a = a_wind(x + halfx*xsh[l],y + halfy*ysh[l],dir[l],user);
if (user->firstorder) {
u_up = (a >= 0) ? au[j + aposj[l]][i + aposi[l]]
: au[j + anegj[l]][i + anegi[l]];
flux[l] = a * u_up;
} else { // use formulas (1.2), (1.3), (1.6) on
// pp 216--217 of Hundsdorfer & Verwer
if (a >= 0.0) {
u_dn = au[j + anegj[l]][i + anegi[l]];
u_up = au[j + aposj[l]][i + aposi[l]];
u_far = au[j + posfarj[l]][i + posfari[l]];
} else {
u_dn = au[j + aposj[l]][i + aposi[l]];
u_up = au[j + anegj[l]][i + anegi[l]];
u_far = au[j + negfarj[l]][i + negfari[l]];
}
flux[l] = a * u_up;
if (u_dn != u_up) {
theta = (u_up - u_far) / (u_dn - u_up);
flux[l] += a * limiter(theta) * (u_dn - u_up);
}
}
}
aG[j][i] = - (flux[0] - flux[2])/hx - (flux[1] - flux[3])/hy
+ g_source(x,y,au[j][i],user);
}
}
return 0;
}
//ENDFUNCTION
PetscErrorCode FormRHSJacobianLocal(DMDALocalInfo *info, double t,
double **au, Mat J, Mat P, AdvectCtx *user) {
PetscErrorCode ierr;
const int dir[4] = {0, 1, 0, 1}, // use x (0) or y (1) component
xsh[4] = { 1, 0,-1, 0}, ysh[4] = { 0, 1, 0,-1};
int i, j, l, ncols;
double hx, hy, halfx, halfy, x, y, a, v[5];
MatStencil col[5],row;
ierr = MatZeroEntries(P); CHKERRQ(ierr);
hx = 1.0 / info->mx; hy = 1.0 / info->my;
halfx = hx / 2.0; halfy = hy / 2.0;
for (j = info->ys; j < info->ys+info->ym; j++) {
y = (j + 0.5) * hy;
row.j = j; col[0].j = j;
for (i = info->xs; i < info->xs+info->xm; i++) {
x = (i + 0.5) * hx;
row.i = i; col[0].i = i;
v[0] = dg_source(x,y,au[j][i],user);
ncols = 1;
for (l = 0; l < 4; l++) { // loop over cell boundaries
a = a_wind(x + halfx*xsh[l],y + halfy*ysh[l],dir[l],user);
// u_up = (a >= 0) ? au[j + aposj[l]][i + aposi[l]]
// : au[j + anegj[l]][i + anegi[l]];
// flux[l] = a * u_up;
if (a >= 0.0) {
switch (l) {
case 0:
col[ncols].j = j; col[ncols].i = i;
v[ncols++] = - a / hx;
break;
case 1:
col[ncols].j = j; col[ncols].i = i;
v[ncols++] = - a / hy;
break;
case 2:
col[ncols].j = j; col[ncols].i = i-1;
v[ncols++] = a / hx;
break;
case 3:
col[ncols].j = j-1; col[ncols].i = i;
v[ncols++] = a / hy;
break;
}
} else { // a < 0
switch (l) {
case 0:
col[ncols].j = j; col[ncols].i = i+1;
v[ncols++] = - a / hx;
break;
case 1:
col[ncols].j = j+1; col[ncols].i = i;
v[ncols++] = - a / hy;
break;
case 2:
col[ncols].j = j; col[ncols].i = i;
v[ncols++] = a / hx;
break;
case 3:
col[ncols].j = j; col[ncols].i = i;
v[ncols++] = a / hy;
break;
}
}
}
ierr = MatSetValuesStencil(P,1,&row,ncols,col,v,ADD_VALUES); CHKERRQ(ierr);
}
}
ierr = MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
ierr = MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
if (J != P) {
ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
}
return 0;
}
int main(int argc,char **argv) {
PetscErrorCode ierr;
AdvectCtx user;
TS ts;
DM da;
Vec u;
DMDALocalInfo info;
double hx, hy, t0, dt;
PetscBool dump = PETSC_FALSE;
char fileroot[PETSC_MAX_PATH_LEN] = "";
PetscInitialize(&argc,&argv,(char*)0,help);
// the wind, and the cone initial condition, are from
// Hundsdorfer & Verwer, "Numerical Solution of Time-Dependent Advection-
// Diffusion-Reaction Equations", Springer 2003, page 303
user.windx = 1.0;
user.windy = 1.0;
user.conex = 0.2;
user.coney = 0.2;
user.coner = 0.1;
user.coneh = 1.0;
user.circlewind = PETSC_FALSE;
user.firstorder = PETSC_FALSE;
ierr = PetscOptionsBegin(PETSC_COMM_WORLD,
"adv_", "options for advect.c", ""); CHKERRQ(ierr);
ierr = PetscOptionsBool("-circlewind","if true, wind is rigid rotation",
"advect.c",user.circlewind,&user.circlewind,NULL);CHKERRQ(ierr);
ierr = PetscOptionsReal("-coneh","cone height",
"advect.c",user.coneh,&user.coneh,NULL);CHKERRQ(ierr);
ierr = PetscOptionsReal("-coner","cone radius",
"advect.c",user.coner,&user.coner,NULL);CHKERRQ(ierr);
ierr = PetscOptionsReal("-conex","x component of cone center",
"advect.c",user.conex,&user.conex,NULL);CHKERRQ(ierr);
ierr = PetscOptionsReal("-coney","y component of cone center",
"advect.c",user.coney,&user.coney,NULL);CHKERRQ(ierr);
ierr = PetscOptionsString("-dumpto","filename root for initial/final state",
"advect.c",fileroot,fileroot,PETSC_MAX_PATH_LEN,&dump);CHKERRQ(ierr);
ierr = PetscOptionsBool("-firstorder","if true, use first-order upwinding",
"advect.c",user.firstorder,&user.firstorder,NULL);CHKERRQ(ierr);
ierr = PetscOptionsReal("-windx","x component of wind (if not circular)",
"advect.c",user.windx,&user.windx,NULL);CHKERRQ(ierr);
ierr = PetscOptionsReal("-windy","y component of wind (if not circular)",
"advect.c",user.windy,&user.windy,NULL);CHKERRQ(ierr);
ierr = PetscOptionsEnd(); CHKERRQ(ierr);
ierr = DMDACreate2d(PETSC_COMM_WORLD,
DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC,
DMDA_STENCIL_STAR, // no diagonal differencing
5,5,PETSC_DECIDE,PETSC_DECIDE, // default to hx=hx=0.2 grid
// (mx=my=5 allows -snes_fd_color)
1, // degrees of freedom
2, // stencil width (flux-limiting)
NULL,NULL,&da); CHKERRQ(ierr);
ierr = DMSetFromOptions(da); CHKERRQ(ierr);
ierr = DMSetUp(da); CHKERRQ(ierr);
ierr = DMSetApplicationContext(da,&user); CHKERRQ(ierr);
// grid is cell-centered
ierr = DMDAGetLocalInfo(da,&info); CHKERRQ(ierr);
hx = 1.0 / info.mx; hy = 1.0 / info.my;
ierr = DMDASetUniformCoordinates(da,
0.0+hx/2.0,1.0-hx/2.0,0.0+hy/2.0,1.0-hy/2.0,0.0,1.0);CHKERRQ(ierr);
ierr = TSCreate(PETSC_COMM_WORLD,&ts); CHKERRQ(ierr);
ierr = TSSetProblemType(ts,TS_NONLINEAR); CHKERRQ(ierr);
ierr = TSSetDM(ts,da); CHKERRQ(ierr);
ierr = DMDATSSetRHSFunctionLocal(da,INSERT_VALUES,
(DMDATSRHSFunctionLocal)FormRHSFunctionLocal,&user); CHKERRQ(ierr);
ierr = DMDATSSetRHSJacobianLocal(da,
(DMDATSRHSJacobianLocal)FormRHSJacobianLocal,&user); CHKERRQ(ierr);
ierr = TSSetType(ts,TSRK); CHKERRQ(ierr);
ierr = TSRKSetType(ts,TSRK2A); CHKERRQ(ierr);
ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP); CHKERRQ(ierr);
ierr = TSSetInitialTimeStep(ts,0.0,0.1); CHKERRQ(ierr);
ierr = TSSetDuration(ts,1000000,0.6); CHKERRQ(ierr);
ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
ierr = TSGetTime(ts,&t0); CHKERRQ(ierr);
ierr = TSGetTimeStep(ts,&dt); CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,
"solving on %d x %d grid with dx=%g x dy=%g cells, t0=%g,\n"
"and initial step dt=%g ...\n",
info.mx,info.my,hx,hy,t0,dt); CHKERRQ(ierr);
ierr = DMCreateGlobalVector(da,&u); CHKERRQ(ierr);
ierr = FormInitial(&info,u,&user); CHKERRQ(ierr);
if (dump) {
PetscViewer viewer;
char filename[PETSC_MAX_PATH_LEN] = "";
ierr = sprintf(filename,"%s_initial.dat",fileroot);
ierr = PetscPrintf(PETSC_COMM_WORLD,
"writing PETSC binary file %s ...\n",filename); CHKERRQ(ierr);
ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,
FILE_MODE_WRITE,&viewer); CHKERRQ(ierr);
ierr = VecView(u,viewer); CHKERRQ(ierr);
ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr);
}
ierr = TSSolve(ts,u); CHKERRQ(ierr);
if (dump) {
PetscViewer viewer;
char filename[PETSC_MAX_PATH_LEN] = "";
ierr = sprintf(filename,"%s_final.dat",fileroot);
ierr = PetscPrintf(PETSC_COMM_WORLD,
"writing PETSC binary file %s ...\n",filename); CHKERRQ(ierr);
ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,
FILE_MODE_WRITE,&viewer); CHKERRQ(ierr);
ierr = VecView(u,viewer); CHKERRQ(ierr);
ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr);
}
VecDestroy(&u); TSDestroy(&ts); DMDestroy(&da);
PetscFinalize();
return 0;
}
