Hi all,
I adapted ex13.c from .../ts/examples/tutorials/ex13.c to solve a 1D TD PDE
(swift-hohenberg).
*Long story short*: the function which sets the Jacobian is never being
called despite being assigned by TSSetRHSJacobian().
I create the TS environment with the following commands (omitting some code
before and after -- Full code attached):
---------------------------------------------------------------------------------------------------------------------------------------------------------------------
ierr = DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 15, 1, 1, NULL,
&da); CHKERRQ(ierr);
/*==================================================================
Let DMDA create
vector
==================================================================*/
ierr = DMCreateGlobalVector(da, &x);CHKERRQ(ierr);
/*==================================================================
Create timestepper context which will use Euler
Method
==================================================================*/
ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr);
ierr = TSSetDM(ts, da);CHKERRQ(ierr);
ierr = TSSetType(ts, TSEULER);CHKERRQ(ierr);
ierr = TSSetRHSFunction(ts, x, RHSFunction, NULL);CHKERRQ(ierr);
/*=================================================================
Create Jacobian for Euler
Method
=================================================================*/
ierr = DMSetMatType(da, MATAIJ);CHKERRQ(ierr);
ierr = DMCreateMatrix(da,&J);CHKERRQ(ierr);
ierr = TSSetRHSJacobian(ts,J,J,RHSJacobian,NULL); CHKERRQ(ierr)
---------------------------------------------------------------------------------------------------------------------------------------------------------------------
I define the function RHSJacobian later, PetscPrintf(PETSC_COMM_WORLD,
"Jacobian function called") is used to see if the function is ever called.
Neither is the message ever called, nor is the Jacobian matrix ever
non-zero. The J matrix is checked with MatView( ). No error messages are
produced either at compiling nor runtime.
Has anyone else run into such an issue or have any advice of how best to
probe the problem?
Thanks so much!
-Scott Dossa
static char help[] = "Solving swift hohenberg equation in 1D \n\
u_t=(r-1)u-2u_xx+u_xxxx in 2D.\n";
#include <petscts.h>
#include <petscdm.h>
#include <petscdmda.h>
extern PetscErrorCode RHSFunction(TS, PetscReal, Vec, Vec, void*);
extern PetscErrorCode InitialValues(DM, Vec, void*);
extern PetscErrorCode RHSJacobian(TS, PetscReal, Vec, Mat, Mat, void*);
extern PetscErrorCode Monitor(TS,PetscInt,PetscReal,Vec,void*);
int main(int argc, char **argv)
{
Vec x;
Mat J;
TS ts;
PetscReal a,b,c,h,ftime, dt;
DM da;
PetscErrorCode ierr;
PetscInt steps, maxsteps=100;
PetscInitialize(&argc, &argv, (char*)0,help);
/*==================================================================
Create 1D distributed array for managing parallel grid and vectors
==================================================================*/
ierr = DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 15, 1, 1, NULL, &da); CHKERRQ(ierr);
/*==================================================================
Let DMDA create vector
==================================================================*/
ierr = DMCreateGlobalVector(da, &x);CHKERRQ(ierr);
/*==================================================================
Create timestepper context to use Euler Method
==================================================================*/
ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr);
ierr = TSSetDM(ts, da);CHKERRQ(ierr);
ierr = TSSetType(ts, TSEULER);CHKERRQ(ierr);
ierr = TSMonitorSet(ts, Monitor, NULL, NULL);
ierr = TSSetRHSFunction(ts, x, RHSFunction, NULL);CHKERRQ(ierr);
/*=================================================================
Create Jacobian for Euler Method
=================================================================*/
ierr = DMSetMatType(da, MATAIJ);CHKERRQ(ierr);
ierr = DMCreateMatrix(da,&J);CHKERRQ(ierr);
ierr = TSSetRHSJacobian(ts,J,J,RHSJacobian,NULL); CHKERRQ(ierr);
ftime = 1.0;
ierr = TSSetDuration(ts, maxsteps, ftime);CHKERRQ(ierr);
ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr);
/*==================================================================
Set initial conditions (see user function InitialValues();)
==================================================================*/
ierr = InitialValues(da, x, NULL);CHKERRQ(ierr);
// ierr = VecView(x, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
dt = 0.0001;
ierr = TSSetInitialTimeStep(ts, 0.0, dt); CHKERRQ(ierr);
/*==================================================================
Set runtime options for TS
==================================================================*/
ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
//ierr = MatView(J, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
/*==================================================================
Solve PDE at each time step
==================================================================*/
ierr = PetscPrintf(PETSC_COMM_WORLD, "time stepping\n\n");CHKERRQ(ierr);
ierr = TSSolve(ts, x);CHKERRQ(ierr);
ierr = TSGetSolveTime(ts, &ftime);CHKERRQ(ierr);
ierr = TSGetTimeStepNumber(ts, &steps);CHKERRQ(ierr);
ierr = MatView(J, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
/*==================================================================
Clean up created objects
==================================================================*/
ierr = VecDestroy(&x);CHKERRQ(ierr);
ierr = MatDestroy(&J);CHKERRQ(ierr);
ierr = TSDestroy(&ts);CHKERRQ(ierr);
ierr = DMDestroy(&da);CHKERRQ(ierr);
PetscFinalize();
PetscFunctionReturn(0);
}
/*==================================================================
Function to evaluate RHS of each time step
==================================================================*/
#undef __FUNCT__
#define __FUNCT__ "RHSFunction"
PetscErrorCode RHSFunction(TS ts, PetscReal ftime, Vec X, Vec F, void *unused)
{
DM da;
PetscErrorCode ierr;
PetscInt i, Mx, xs, xm;
PetscReal h;
PetscScalar u, uxx, uxxxx, *uarray, *f;
Vec localX;
/* Always begin user defined functions with this call */
PetscFunctionBeginUser;
/* Access DM for TS */
ierr = TSGetDM(ts, &da); CHKERRQ(ierr);
/* Get a local part of the vecor x */
ierr = DMGetLocalVector(da, &localX);CHKERRQ(ierr);
/* Access info about da */
ierr = DMDAGetInfo(da, PETSC_IGNORE, &Mx, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE); CHKERRQ(ierr);
/* Scatter ghost points */
ierr = DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX);CHKERRQ(ierr);
ierr = DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX);CHKERRQ(ierr);
/* Give local DM access to local vector data */
ierr = DMDAVecGetArrayRead(da, localX, &uarray);CHKERRQ(ierr);
ierr = DMDAVecGetArray(da, F, &f);CHKERRQ(ierr);
/* Get where the local grid boundaries are */
ierr = DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL);CHKERRQ(ierr);
/* Compute point separation in discretization */
h = 1/(Mx-1);
// ierr = PetscPrintf(PETSC_COMM_WORLD, "Calculating RHS Functino\n\n");CHKERRQ(ierr);
/* Compute the function over the local grid points */
for (i=xs; i<xs+xm; i++) {
if(i==0 || i==Mx-1) {
f[i]=uarray[i];
continue;
}
u = uarray[i];
uxx = (uarray[i+1]-2.0*uarray[i]+uarray[i-1])/(h*h);
uxxxx = (uarray[i+2]-4.0*uarray[i+1]+6.0*uarray[i]-4.0*uarray[i-1]+uarray[i-2])/(h*h*h*h);
f[i] =u-2.0*uxx+uxxxx;
}
/* Restore the vectors */
ierr = DMDAVecRestoreArrayRead(da,localX, &uarray);CHKERRQ(ierr);
ierr = DMDAVecRestoreArray(da, F, &f);CHKERRQ(ierr);
ierr = DMRestoreLocalVector(da, &localX);CHKERRQ(ierr);
/* end user defined function */
PetscFunctionReturn(0);
}
/*==================================================================
User defined function to define initial values
==================================================================*/
#undef __FUNCT__
#define __FUNCT__ "InitialValues"
PetscErrorCode InitialValues(DM da, Vec X, void *unused)
{
PetscErrorCode ierr;
PetscInt i, xs, xm, Mx;
PetscScalar *x;
PetscReal h;
/* Always start user defined functions with this call */
PetscFunctionBeginUser;
/* Get DMDA info for global and local vector lengths */
ierr = DMDAGetInfo(da, PETSC_IGNORE, &Mx, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE);
/* Get vector data for local operations */
ierr = DMDAVecGetArray(da, X, &x); CHKERRQ(ierr);
/* Get corner data to know length of local vector */
ierr = DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL); CHKERRQ(ierr);
/* Set initial values */
h = 1./(Mx-1);
for (i=xs; i<xs+xm; i++) {
x[i] =100.*sin(2*3.1415*i*h);
}
/* Restore array to update values */
ierr = DMDAVecRestoreArray(da, X, &x);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
/*==================================================================
User provided routine for RHS's Jacobian
==================================================================*/
#undef __FUNCT__
#define __FUNCT__ "RHSJacobian"
PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec X, Mat J, Mat Jpre, void *unused)
{
PetscErrorCode ierr;
DM da;
DMDALocalInfo info;
PetscInt i;
PetscReal h;
/* Always start user defined functions with this call */
PetscFunctionBeginUser;
/* Get TS info and DMDA info */
ierr = TSGetDM(ts, &da);CHKERRQ(ierr);
ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr);
/* Set local grid spacing */
h = 1.0/(info.mx-1);
PetscPrintf(PETSC_COMM_WORLD,"Definining Jacobian");
/*for loop for assigning matrix values */
for (i=info.xs; i<info.xs+info.xm;i++) {
PetscInt nc=0;
MatStencil row,col[5];
PetscScalar val[5];
row.i=i;
if (i==0 || i==info.mx-1) {
col[nc].i=i; val[nc++]=1.0;
} else {
col[nc].i=i-2; val[nc++]=1./(h*h*h*h);
col[nc].i=i-1; val[nc++]=-2./(h*h)-4./(h*h*h*h);
col[nc].i=i; val[nc++]=1+4./(h*h)+6./(h*h*h*h);
col[nc].i=i+1; val[nc++]=-2./(h*h)-4./(h*h*h*h);
col[nc].i=i+2; val[nc++]=1./(h*h*h*h);
}
ierr = MatSetValuesStencil(Jpre, 1, &row, nc, col, val, INSERT_VALUES);CHKERRQ(ierr);
}
ierr = MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
if (J!=Jpre) {
ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
}
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "Monitor"
PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal time, Vec X, void* unused){
PetscFunctionBeginUser;
VecView(X, PETSC_VIEWER_DRAW_WORLD);
PetscFunctionReturn(0);
}
