1*c4762a1bSJed Brown 2*c4762a1bSJed Brown static char help[] ="Nonlinear Radiative Transport PDE with multigrid in 2d.\n\ 3*c4762a1bSJed Brown Uses 2-dimensional distributed arrays.\n\ 4*c4762a1bSJed Brown A 2-dim simplified Radiative Transport test problem is used, with analytic Jacobian. \n\ 5*c4762a1bSJed Brown \n\ 6*c4762a1bSJed Brown Solves the linear systems via multilevel methods \n\ 7*c4762a1bSJed Brown \n\ 8*c4762a1bSJed Brown The command line\n\ 9*c4762a1bSJed Brown options are:\n\ 10*c4762a1bSJed Brown -tleft <tl>, where <tl> indicates the left Diriclet BC \n\ 11*c4762a1bSJed Brown -tright <tr>, where <tr> indicates the right Diriclet BC \n\ 12*c4762a1bSJed Brown -beta <beta>, where <beta> indicates the exponent in T \n\n"; 13*c4762a1bSJed Brown 14*c4762a1bSJed Brown /*T 15*c4762a1bSJed Brown Concepts: SNES^solving a system of nonlinear equations 16*c4762a1bSJed Brown Concepts: DMDA^using distributed arrays 17*c4762a1bSJed Brown Concepts: multigrid; 18*c4762a1bSJed Brown Processors: n 19*c4762a1bSJed Brown T*/ 20*c4762a1bSJed Brown 21*c4762a1bSJed Brown 22*c4762a1bSJed Brown 23*c4762a1bSJed Brown /* 24*c4762a1bSJed Brown 25*c4762a1bSJed Brown This example models the partial differential equation 26*c4762a1bSJed Brown 27*c4762a1bSJed Brown - Div(alpha* T^beta (GRAD T)) = 0. 28*c4762a1bSJed Brown 29*c4762a1bSJed Brown where beta = 2.5 and alpha = 1.0 30*c4762a1bSJed Brown 31*c4762a1bSJed Brown BC: T_left = 1.0, T_right = 0.1, dT/dn_top = dTdn_bottom = 0. 32*c4762a1bSJed Brown 33*c4762a1bSJed Brown in the unit square, which is uniformly discretized in each of x and 34*c4762a1bSJed Brown y in this simple encoding. The degrees of freedom are cell centered. 35*c4762a1bSJed Brown 36*c4762a1bSJed Brown A finite volume approximation with the usual 5-point stencil 37*c4762a1bSJed Brown is used to discretize the boundary value problem to obtain a 38*c4762a1bSJed Brown nonlinear system of equations. 39*c4762a1bSJed Brown 40*c4762a1bSJed Brown This code was contributed by David Keyes 41*c4762a1bSJed Brown 42*c4762a1bSJed Brown */ 43*c4762a1bSJed Brown 44*c4762a1bSJed Brown #include <petscsnes.h> 45*c4762a1bSJed Brown #include <petscdm.h> 46*c4762a1bSJed Brown #include <petscdmda.h> 47*c4762a1bSJed Brown 48*c4762a1bSJed Brown /* User-defined application context */ 49*c4762a1bSJed Brown 50*c4762a1bSJed Brown typedef struct { 51*c4762a1bSJed Brown PetscReal tleft,tright; /* Dirichlet boundary conditions */ 52*c4762a1bSJed Brown PetscReal beta,bm1,coef; /* nonlinear diffusivity parameterizations */ 53*c4762a1bSJed Brown } AppCtx; 54*c4762a1bSJed Brown 55*c4762a1bSJed Brown #define POWFLOP 5 /* assume a pow() takes five flops */ 56*c4762a1bSJed Brown 57*c4762a1bSJed Brown extern PetscErrorCode FormInitialGuess(SNES,Vec,void*); 58*c4762a1bSJed Brown extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*); 59*c4762a1bSJed Brown extern PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*); 60*c4762a1bSJed Brown 61*c4762a1bSJed Brown int main(int argc,char **argv) 62*c4762a1bSJed Brown { 63*c4762a1bSJed Brown SNES snes; 64*c4762a1bSJed Brown AppCtx user; 65*c4762a1bSJed Brown PetscErrorCode ierr; 66*c4762a1bSJed Brown PetscInt its,lits; 67*c4762a1bSJed Brown PetscReal litspit; 68*c4762a1bSJed Brown DM da; 69*c4762a1bSJed Brown 70*c4762a1bSJed Brown ierr = PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr; 71*c4762a1bSJed Brown 72*c4762a1bSJed Brown /* set problem parameters */ 73*c4762a1bSJed Brown user.tleft = 1.0; 74*c4762a1bSJed Brown user.tright = 0.1; 75*c4762a1bSJed Brown user.beta = 2.5; 76*c4762a1bSJed Brown ierr = PetscOptionsGetReal(NULL,NULL,"-tleft",&user.tleft,NULL);CHKERRQ(ierr); 77*c4762a1bSJed Brown ierr = PetscOptionsGetReal(NULL,NULL,"-tright",&user.tright,NULL);CHKERRQ(ierr); 78*c4762a1bSJed Brown ierr = PetscOptionsGetReal(NULL,NULL,"-beta",&user.beta,NULL);CHKERRQ(ierr); 79*c4762a1bSJed Brown user.bm1 = user.beta - 1.0; 80*c4762a1bSJed Brown user.coef = user.beta/2.0; 81*c4762a1bSJed Brown 82*c4762a1bSJed Brown /* 83*c4762a1bSJed Brown Create the multilevel DM data structure 84*c4762a1bSJed Brown */ 85*c4762a1bSJed Brown ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr); 86*c4762a1bSJed Brown 87*c4762a1bSJed Brown /* 88*c4762a1bSJed Brown Set the DMDA (grid structure) for the grids. 89*c4762a1bSJed Brown */ 90*c4762a1bSJed Brown ierr = DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,5,5,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&da);CHKERRQ(ierr); 91*c4762a1bSJed Brown ierr = DMSetFromOptions(da);CHKERRQ(ierr); 92*c4762a1bSJed Brown ierr = DMSetUp(da);CHKERRQ(ierr); 93*c4762a1bSJed Brown ierr = DMSetApplicationContext(da,&user);CHKERRQ(ierr); 94*c4762a1bSJed Brown ierr = SNESSetDM(snes,(DM)da);CHKERRQ(ierr); 95*c4762a1bSJed Brown 96*c4762a1bSJed Brown /* 97*c4762a1bSJed Brown Create the nonlinear solver, and tell it the functions to use 98*c4762a1bSJed Brown */ 99*c4762a1bSJed Brown ierr = SNESSetFunction(snes,NULL,FormFunction,&user);CHKERRQ(ierr); 100*c4762a1bSJed Brown ierr = SNESSetJacobian(snes,NULL,NULL,FormJacobian,&user);CHKERRQ(ierr); 101*c4762a1bSJed Brown ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 102*c4762a1bSJed Brown ierr = SNESSetComputeInitialGuess(snes,FormInitialGuess,NULL);CHKERRQ(ierr); 103*c4762a1bSJed Brown 104*c4762a1bSJed Brown ierr = SNESSolve(snes,NULL,NULL);CHKERRQ(ierr); 105*c4762a1bSJed Brown ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr); 106*c4762a1bSJed Brown ierr = SNESGetLinearSolveIterations(snes,&lits);CHKERRQ(ierr); 107*c4762a1bSJed Brown litspit = ((PetscReal)lits)/((PetscReal)its); 108*c4762a1bSJed Brown ierr = PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n",its);CHKERRQ(ierr); 109*c4762a1bSJed Brown ierr = PetscPrintf(PETSC_COMM_WORLD,"Number of Linear iterations = %D\n",lits);CHKERRQ(ierr); 110*c4762a1bSJed Brown ierr = PetscPrintf(PETSC_COMM_WORLD,"Average Linear its / SNES = %e\n",(double)litspit);CHKERRQ(ierr); 111*c4762a1bSJed Brown 112*c4762a1bSJed Brown ierr = DMDestroy(&da);CHKERRQ(ierr); 113*c4762a1bSJed Brown ierr = SNESDestroy(&snes);CHKERRQ(ierr); 114*c4762a1bSJed Brown ierr = PetscFinalize(); 115*c4762a1bSJed Brown return ierr; 116*c4762a1bSJed Brown } 117*c4762a1bSJed Brown /* -------------------- Form initial approximation ----------------- */ 118*c4762a1bSJed Brown PetscErrorCode FormInitialGuess(SNES snes,Vec X,void *ctx) 119*c4762a1bSJed Brown { 120*c4762a1bSJed Brown AppCtx *user; 121*c4762a1bSJed Brown PetscInt i,j,xs,ys,xm,ym; 122*c4762a1bSJed Brown PetscErrorCode ierr; 123*c4762a1bSJed Brown PetscReal tleft; 124*c4762a1bSJed Brown PetscScalar **x; 125*c4762a1bSJed Brown DM da; 126*c4762a1bSJed Brown 127*c4762a1bSJed Brown PetscFunctionBeginUser; 128*c4762a1bSJed Brown ierr = SNESGetDM(snes,&da);CHKERRQ(ierr); 129*c4762a1bSJed Brown ierr = DMGetApplicationContext(da,&user);CHKERRQ(ierr); 130*c4762a1bSJed Brown tleft = user->tleft; 131*c4762a1bSJed Brown /* Get ghost points */ 132*c4762a1bSJed Brown ierr = DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);CHKERRQ(ierr); 133*c4762a1bSJed Brown ierr = DMDAVecGetArray(da,X,&x);CHKERRQ(ierr); 134*c4762a1bSJed Brown 135*c4762a1bSJed Brown /* Compute initial guess */ 136*c4762a1bSJed Brown for (j=ys; j<ys+ym; j++) { 137*c4762a1bSJed Brown for (i=xs; i<xs+xm; i++) { 138*c4762a1bSJed Brown x[j][i] = tleft; 139*c4762a1bSJed Brown } 140*c4762a1bSJed Brown } 141*c4762a1bSJed Brown ierr = DMDAVecRestoreArray(da,X,&x);CHKERRQ(ierr); 142*c4762a1bSJed Brown PetscFunctionReturn(0); 143*c4762a1bSJed Brown } 144*c4762a1bSJed Brown /* -------------------- Evaluate Function F(x) --------------------- */ 145*c4762a1bSJed Brown PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ptr) 146*c4762a1bSJed Brown { 147*c4762a1bSJed Brown AppCtx *user = (AppCtx*)ptr; 148*c4762a1bSJed Brown PetscErrorCode ierr; 149*c4762a1bSJed Brown PetscInt i,j,mx,my,xs,ys,xm,ym; 150*c4762a1bSJed Brown PetscScalar zero = 0.0,one = 1.0; 151*c4762a1bSJed Brown PetscScalar hx,hy,hxdhy,hydhx; 152*c4762a1bSJed Brown PetscScalar t0,tn,ts,te,tw,an,as,ae,aw,dn,ds,de,dw,fn = 0.0,fs = 0.0,fe =0.0,fw = 0.0; 153*c4762a1bSJed Brown PetscScalar tleft,tright,beta; 154*c4762a1bSJed Brown PetscScalar **x,**f; 155*c4762a1bSJed Brown Vec localX; 156*c4762a1bSJed Brown DM da; 157*c4762a1bSJed Brown 158*c4762a1bSJed Brown PetscFunctionBeginUser; 159*c4762a1bSJed Brown ierr = SNESGetDM(snes,&da);CHKERRQ(ierr); 160*c4762a1bSJed Brown ierr = DMGetLocalVector(da,&localX);CHKERRQ(ierr); 161*c4762a1bSJed Brown ierr = DMDAGetInfo(da,NULL,&mx,&my,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 162*c4762a1bSJed Brown hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1); 163*c4762a1bSJed Brown hxdhy = hx/hy; hydhx = hy/hx; 164*c4762a1bSJed Brown tleft = user->tleft; tright = user->tright; 165*c4762a1bSJed Brown beta = user->beta; 166*c4762a1bSJed Brown 167*c4762a1bSJed Brown /* Get ghost points */ 168*c4762a1bSJed Brown ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);CHKERRQ(ierr); 169*c4762a1bSJed Brown ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);CHKERRQ(ierr); 170*c4762a1bSJed Brown ierr = DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);CHKERRQ(ierr); 171*c4762a1bSJed Brown ierr = DMDAVecGetArray(da,localX,&x);CHKERRQ(ierr); 172*c4762a1bSJed Brown ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr); 173*c4762a1bSJed Brown 174*c4762a1bSJed Brown /* Evaluate function */ 175*c4762a1bSJed Brown for (j=ys; j<ys+ym; j++) { 176*c4762a1bSJed Brown for (i=xs; i<xs+xm; i++) { 177*c4762a1bSJed Brown t0 = x[j][i]; 178*c4762a1bSJed Brown 179*c4762a1bSJed Brown if (i > 0 && i < mx-1 && j > 0 && j < my-1) { 180*c4762a1bSJed Brown 181*c4762a1bSJed Brown /* general interior volume */ 182*c4762a1bSJed Brown 183*c4762a1bSJed Brown tw = x[j][i-1]; 184*c4762a1bSJed Brown aw = 0.5*(t0 + tw); 185*c4762a1bSJed Brown dw = PetscPowScalar(aw,beta); 186*c4762a1bSJed Brown fw = dw*(t0 - tw); 187*c4762a1bSJed Brown 188*c4762a1bSJed Brown te = x[j][i+1]; 189*c4762a1bSJed Brown ae = 0.5*(t0 + te); 190*c4762a1bSJed Brown de = PetscPowScalar(ae,beta); 191*c4762a1bSJed Brown fe = de*(te - t0); 192*c4762a1bSJed Brown 193*c4762a1bSJed Brown ts = x[j-1][i]; 194*c4762a1bSJed Brown as = 0.5*(t0 + ts); 195*c4762a1bSJed Brown ds = PetscPowScalar(as,beta); 196*c4762a1bSJed Brown fs = ds*(t0 - ts); 197*c4762a1bSJed Brown 198*c4762a1bSJed Brown tn = x[j+1][i]; 199*c4762a1bSJed Brown an = 0.5*(t0 + tn); 200*c4762a1bSJed Brown dn = PetscPowScalar(an,beta); 201*c4762a1bSJed Brown fn = dn*(tn - t0); 202*c4762a1bSJed Brown 203*c4762a1bSJed Brown } else if (i == 0) { 204*c4762a1bSJed Brown 205*c4762a1bSJed Brown /* left-hand boundary */ 206*c4762a1bSJed Brown tw = tleft; 207*c4762a1bSJed Brown aw = 0.5*(t0 + tw); 208*c4762a1bSJed Brown dw = PetscPowScalar(aw,beta); 209*c4762a1bSJed Brown fw = dw*(t0 - tw); 210*c4762a1bSJed Brown 211*c4762a1bSJed Brown te = x[j][i+1]; 212*c4762a1bSJed Brown ae = 0.5*(t0 + te); 213*c4762a1bSJed Brown de = PetscPowScalar(ae,beta); 214*c4762a1bSJed Brown fe = de*(te - t0); 215*c4762a1bSJed Brown 216*c4762a1bSJed Brown if (j > 0) { 217*c4762a1bSJed Brown ts = x[j-1][i]; 218*c4762a1bSJed Brown as = 0.5*(t0 + ts); 219*c4762a1bSJed Brown ds = PetscPowScalar(as,beta); 220*c4762a1bSJed Brown fs = ds*(t0 - ts); 221*c4762a1bSJed Brown } else fs = zero; 222*c4762a1bSJed Brown 223*c4762a1bSJed Brown if (j < my-1) { 224*c4762a1bSJed Brown tn = x[j+1][i]; 225*c4762a1bSJed Brown an = 0.5*(t0 + tn); 226*c4762a1bSJed Brown dn = PetscPowScalar(an,beta); 227*c4762a1bSJed Brown fn = dn*(tn - t0); 228*c4762a1bSJed Brown } else fn = zero; 229*c4762a1bSJed Brown 230*c4762a1bSJed Brown } else if (i == mx-1) { 231*c4762a1bSJed Brown 232*c4762a1bSJed Brown /* right-hand boundary */ 233*c4762a1bSJed Brown tw = x[j][i-1]; 234*c4762a1bSJed Brown aw = 0.5*(t0 + tw); 235*c4762a1bSJed Brown dw = PetscPowScalar(aw,beta); 236*c4762a1bSJed Brown fw = dw*(t0 - tw); 237*c4762a1bSJed Brown 238*c4762a1bSJed Brown te = tright; 239*c4762a1bSJed Brown ae = 0.5*(t0 + te); 240*c4762a1bSJed Brown de = PetscPowScalar(ae,beta); 241*c4762a1bSJed Brown fe = de*(te - t0); 242*c4762a1bSJed Brown 243*c4762a1bSJed Brown if (j > 0) { 244*c4762a1bSJed Brown ts = x[j-1][i]; 245*c4762a1bSJed Brown as = 0.5*(t0 + ts); 246*c4762a1bSJed Brown ds = PetscPowScalar(as,beta); 247*c4762a1bSJed Brown fs = ds*(t0 - ts); 248*c4762a1bSJed Brown } else fs = zero; 249*c4762a1bSJed Brown 250*c4762a1bSJed Brown if (j < my-1) { 251*c4762a1bSJed Brown tn = x[j+1][i]; 252*c4762a1bSJed Brown an = 0.5*(t0 + tn); 253*c4762a1bSJed Brown dn = PetscPowScalar(an,beta); 254*c4762a1bSJed Brown fn = dn*(tn - t0); 255*c4762a1bSJed Brown } else fn = zero; 256*c4762a1bSJed Brown 257*c4762a1bSJed Brown } else if (j == 0) { 258*c4762a1bSJed Brown 259*c4762a1bSJed Brown /* bottom boundary,and i <> 0 or mx-1 */ 260*c4762a1bSJed Brown tw = x[j][i-1]; 261*c4762a1bSJed Brown aw = 0.5*(t0 + tw); 262*c4762a1bSJed Brown dw = PetscPowScalar(aw,beta); 263*c4762a1bSJed Brown fw = dw*(t0 - tw); 264*c4762a1bSJed Brown 265*c4762a1bSJed Brown te = x[j][i+1]; 266*c4762a1bSJed Brown ae = 0.5*(t0 + te); 267*c4762a1bSJed Brown de = PetscPowScalar(ae,beta); 268*c4762a1bSJed Brown fe = de*(te - t0); 269*c4762a1bSJed Brown 270*c4762a1bSJed Brown fs = zero; 271*c4762a1bSJed Brown 272*c4762a1bSJed Brown tn = x[j+1][i]; 273*c4762a1bSJed Brown an = 0.5*(t0 + tn); 274*c4762a1bSJed Brown dn = PetscPowScalar(an,beta); 275*c4762a1bSJed Brown fn = dn*(tn - t0); 276*c4762a1bSJed Brown 277*c4762a1bSJed Brown } else if (j == my-1) { 278*c4762a1bSJed Brown 279*c4762a1bSJed Brown /* top boundary,and i <> 0 or mx-1 */ 280*c4762a1bSJed Brown tw = x[j][i-1]; 281*c4762a1bSJed Brown aw = 0.5*(t0 + tw); 282*c4762a1bSJed Brown dw = PetscPowScalar(aw,beta); 283*c4762a1bSJed Brown fw = dw*(t0 - tw); 284*c4762a1bSJed Brown 285*c4762a1bSJed Brown te = x[j][i+1]; 286*c4762a1bSJed Brown ae = 0.5*(t0 + te); 287*c4762a1bSJed Brown de = PetscPowScalar(ae,beta); 288*c4762a1bSJed Brown fe = de*(te - t0); 289*c4762a1bSJed Brown 290*c4762a1bSJed Brown ts = x[j-1][i]; 291*c4762a1bSJed Brown as = 0.5*(t0 + ts); 292*c4762a1bSJed Brown ds = PetscPowScalar(as,beta); 293*c4762a1bSJed Brown fs = ds*(t0 - ts); 294*c4762a1bSJed Brown 295*c4762a1bSJed Brown fn = zero; 296*c4762a1bSJed Brown 297*c4762a1bSJed Brown } 298*c4762a1bSJed Brown 299*c4762a1bSJed Brown f[j][i] = -hydhx*(fe-fw) - hxdhy*(fn-fs); 300*c4762a1bSJed Brown 301*c4762a1bSJed Brown } 302*c4762a1bSJed Brown } 303*c4762a1bSJed Brown ierr = DMDAVecRestoreArray(da,localX,&x);CHKERRQ(ierr); 304*c4762a1bSJed Brown ierr = DMDAVecRestoreArray(da,F,&f);CHKERRQ(ierr); 305*c4762a1bSJed Brown ierr = DMRestoreLocalVector(da,&localX);CHKERRQ(ierr); 306*c4762a1bSJed Brown ierr = PetscLogFlops((22.0 + 4.0*POWFLOP)*ym*xm);CHKERRQ(ierr); 307*c4762a1bSJed Brown PetscFunctionReturn(0); 308*c4762a1bSJed Brown } 309*c4762a1bSJed Brown /* -------------------- Evaluate Jacobian F(x) --------------------- */ 310*c4762a1bSJed Brown PetscErrorCode FormJacobian(SNES snes,Vec X,Mat jac,Mat B,void *ptr) 311*c4762a1bSJed Brown { 312*c4762a1bSJed Brown AppCtx *user = (AppCtx*)ptr; 313*c4762a1bSJed Brown PetscErrorCode ierr; 314*c4762a1bSJed Brown PetscInt i,j,mx,my,xs,ys,xm,ym; 315*c4762a1bSJed Brown PetscScalar one = 1.0,hx,hy,hxdhy,hydhx,t0,tn,ts,te,tw; 316*c4762a1bSJed Brown PetscScalar dn,ds,de,dw,an,as,ae,aw,bn,bs,be,bw,gn,gs,ge,gw; 317*c4762a1bSJed Brown PetscScalar tleft,tright,beta,bm1,coef; 318*c4762a1bSJed Brown PetscScalar v[5],**x; 319*c4762a1bSJed Brown Vec localX; 320*c4762a1bSJed Brown MatStencil col[5],row; 321*c4762a1bSJed Brown DM da; 322*c4762a1bSJed Brown 323*c4762a1bSJed Brown PetscFunctionBeginUser; 324*c4762a1bSJed Brown ierr = SNESGetDM(snes,&da);CHKERRQ(ierr); 325*c4762a1bSJed Brown ierr = DMGetLocalVector(da,&localX);CHKERRQ(ierr); 326*c4762a1bSJed Brown ierr = DMDAGetInfo(da,NULL,&mx,&my,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 327*c4762a1bSJed Brown hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1); 328*c4762a1bSJed Brown hxdhy = hx/hy; hydhx = hy/hx; 329*c4762a1bSJed Brown tleft = user->tleft; tright = user->tright; 330*c4762a1bSJed Brown beta = user->beta; bm1 = user->bm1; coef = user->coef; 331*c4762a1bSJed Brown 332*c4762a1bSJed Brown /* Get ghost points */ 333*c4762a1bSJed Brown ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);CHKERRQ(ierr); 334*c4762a1bSJed Brown ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);CHKERRQ(ierr); 335*c4762a1bSJed Brown ierr = DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);CHKERRQ(ierr); 336*c4762a1bSJed Brown ierr = DMDAVecGetArray(da,localX,&x);CHKERRQ(ierr); 337*c4762a1bSJed Brown 338*c4762a1bSJed Brown /* Evaluate Jacobian of function */ 339*c4762a1bSJed Brown for (j=ys; j<ys+ym; j++) { 340*c4762a1bSJed Brown for (i=xs; i<xs+xm; i++) { 341*c4762a1bSJed Brown t0 = x[j][i]; 342*c4762a1bSJed Brown 343*c4762a1bSJed Brown if (i > 0 && i < mx-1 && j > 0 && j < my-1) { 344*c4762a1bSJed Brown 345*c4762a1bSJed Brown /* general interior volume */ 346*c4762a1bSJed Brown 347*c4762a1bSJed Brown tw = x[j][i-1]; 348*c4762a1bSJed Brown aw = 0.5*(t0 + tw); 349*c4762a1bSJed Brown bw = PetscPowScalar(aw,bm1); 350*c4762a1bSJed Brown /* dw = bw * aw */ 351*c4762a1bSJed Brown dw = PetscPowScalar(aw,beta); 352*c4762a1bSJed Brown gw = coef*bw*(t0 - tw); 353*c4762a1bSJed Brown 354*c4762a1bSJed Brown te = x[j][i+1]; 355*c4762a1bSJed Brown ae = 0.5*(t0 + te); 356*c4762a1bSJed Brown be = PetscPowScalar(ae,bm1); 357*c4762a1bSJed Brown /* de = be * ae; */ 358*c4762a1bSJed Brown de = PetscPowScalar(ae,beta); 359*c4762a1bSJed Brown ge = coef*be*(te - t0); 360*c4762a1bSJed Brown 361*c4762a1bSJed Brown ts = x[j-1][i]; 362*c4762a1bSJed Brown as = 0.5*(t0 + ts); 363*c4762a1bSJed Brown bs = PetscPowScalar(as,bm1); 364*c4762a1bSJed Brown /* ds = bs * as; */ 365*c4762a1bSJed Brown ds = PetscPowScalar(as,beta); 366*c4762a1bSJed Brown gs = coef*bs*(t0 - ts); 367*c4762a1bSJed Brown 368*c4762a1bSJed Brown tn = x[j+1][i]; 369*c4762a1bSJed Brown an = 0.5*(t0 + tn); 370*c4762a1bSJed Brown bn = PetscPowScalar(an,bm1); 371*c4762a1bSJed Brown /* dn = bn * an; */ 372*c4762a1bSJed Brown dn = PetscPowScalar(an,beta); 373*c4762a1bSJed Brown gn = coef*bn*(tn - t0); 374*c4762a1bSJed Brown 375*c4762a1bSJed Brown v[0] = -hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i; 376*c4762a1bSJed Brown v[1] = -hydhx*(dw - gw); col[1].j = j; col[1].i = i-1; 377*c4762a1bSJed Brown v[2] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge); col[2].j = row.j = j; col[2].i = row.i = i; 378*c4762a1bSJed Brown v[3] = -hydhx*(de + ge); col[3].j = j; col[3].i = i+1; 379*c4762a1bSJed Brown v[4] = -hxdhy*(dn + gn); col[4].j = j+1; col[4].i = i; 380*c4762a1bSJed Brown ierr = MatSetValuesStencil(B,1,&row,5,col,v,INSERT_VALUES);CHKERRQ(ierr); 381*c4762a1bSJed Brown 382*c4762a1bSJed Brown } else if (i == 0) { 383*c4762a1bSJed Brown 384*c4762a1bSJed Brown /* left-hand boundary */ 385*c4762a1bSJed Brown tw = tleft; 386*c4762a1bSJed Brown aw = 0.5*(t0 + tw); 387*c4762a1bSJed Brown bw = PetscPowScalar(aw,bm1); 388*c4762a1bSJed Brown /* dw = bw * aw */ 389*c4762a1bSJed Brown dw = PetscPowScalar(aw,beta); 390*c4762a1bSJed Brown gw = coef*bw*(t0 - tw); 391*c4762a1bSJed Brown 392*c4762a1bSJed Brown te = x[j][i + 1]; 393*c4762a1bSJed Brown ae = 0.5*(t0 + te); 394*c4762a1bSJed Brown be = PetscPowScalar(ae,bm1); 395*c4762a1bSJed Brown /* de = be * ae; */ 396*c4762a1bSJed Brown de = PetscPowScalar(ae,beta); 397*c4762a1bSJed Brown ge = coef*be*(te - t0); 398*c4762a1bSJed Brown 399*c4762a1bSJed Brown /* left-hand bottom boundary */ 400*c4762a1bSJed Brown if (j == 0) { 401*c4762a1bSJed Brown 402*c4762a1bSJed Brown tn = x[j+1][i]; 403*c4762a1bSJed Brown an = 0.5*(t0 + tn); 404*c4762a1bSJed Brown bn = PetscPowScalar(an,bm1); 405*c4762a1bSJed Brown /* dn = bn * an; */ 406*c4762a1bSJed Brown dn = PetscPowScalar(an,beta); 407*c4762a1bSJed Brown gn = coef*bn*(tn - t0); 408*c4762a1bSJed Brown 409*c4762a1bSJed Brown v[0] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[0].j = row.j = j; col[0].i = row.i = i; 410*c4762a1bSJed Brown v[1] = -hydhx*(de + ge); col[1].j = j; col[1].i = i+1; 411*c4762a1bSJed Brown v[2] = -hxdhy*(dn + gn); col[2].j = j+1; col[2].i = i; 412*c4762a1bSJed Brown ierr = MatSetValuesStencil(B,1,&row,3,col,v,INSERT_VALUES);CHKERRQ(ierr); 413*c4762a1bSJed Brown 414*c4762a1bSJed Brown /* left-hand interior boundary */ 415*c4762a1bSJed Brown } else if (j < my-1) { 416*c4762a1bSJed Brown 417*c4762a1bSJed Brown ts = x[j-1][i]; 418*c4762a1bSJed Brown as = 0.5*(t0 + ts); 419*c4762a1bSJed Brown bs = PetscPowScalar(as,bm1); 420*c4762a1bSJed Brown /* ds = bs * as; */ 421*c4762a1bSJed Brown ds = PetscPowScalar(as,beta); 422*c4762a1bSJed Brown gs = coef*bs*(t0 - ts); 423*c4762a1bSJed Brown 424*c4762a1bSJed Brown tn = x[j+1][i]; 425*c4762a1bSJed Brown an = 0.5*(t0 + tn); 426*c4762a1bSJed Brown bn = PetscPowScalar(an,bm1); 427*c4762a1bSJed Brown /* dn = bn * an; */ 428*c4762a1bSJed Brown dn = PetscPowScalar(an,beta); 429*c4762a1bSJed Brown gn = coef*bn*(tn - t0); 430*c4762a1bSJed Brown 431*c4762a1bSJed Brown v[0] = -hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i; 432*c4762a1bSJed Brown v[1] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i; 433*c4762a1bSJed Brown v[2] = -hydhx*(de + ge); col[2].j = j; col[2].i = i+1; 434*c4762a1bSJed Brown v[3] = -hxdhy*(dn + gn); col[3].j = j+1; col[3].i = i; 435*c4762a1bSJed Brown ierr = MatSetValuesStencil(B,1,&row,4,col,v,INSERT_VALUES);CHKERRQ(ierr); 436*c4762a1bSJed Brown /* left-hand top boundary */ 437*c4762a1bSJed Brown } else { 438*c4762a1bSJed Brown 439*c4762a1bSJed Brown ts = x[j-1][i]; 440*c4762a1bSJed Brown as = 0.5*(t0 + ts); 441*c4762a1bSJed Brown bs = PetscPowScalar(as,bm1); 442*c4762a1bSJed Brown /* ds = bs * as; */ 443*c4762a1bSJed Brown ds = PetscPowScalar(as,beta); 444*c4762a1bSJed Brown gs = coef*bs*(t0 - ts); 445*c4762a1bSJed Brown 446*c4762a1bSJed Brown v[0] = -hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i; 447*c4762a1bSJed Brown v[1] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i; 448*c4762a1bSJed Brown v[2] = -hydhx*(de + ge); col[2].j = j; col[2].i = i+1; 449*c4762a1bSJed Brown ierr = MatSetValuesStencil(B,1,&row,3,col,v,INSERT_VALUES);CHKERRQ(ierr); 450*c4762a1bSJed Brown } 451*c4762a1bSJed Brown 452*c4762a1bSJed Brown } else if (i == mx-1) { 453*c4762a1bSJed Brown 454*c4762a1bSJed Brown /* right-hand boundary */ 455*c4762a1bSJed Brown tw = x[j][i-1]; 456*c4762a1bSJed Brown aw = 0.5*(t0 + tw); 457*c4762a1bSJed Brown bw = PetscPowScalar(aw,bm1); 458*c4762a1bSJed Brown /* dw = bw * aw */ 459*c4762a1bSJed Brown dw = PetscPowScalar(aw,beta); 460*c4762a1bSJed Brown gw = coef*bw*(t0 - tw); 461*c4762a1bSJed Brown 462*c4762a1bSJed Brown te = tright; 463*c4762a1bSJed Brown ae = 0.5*(t0 + te); 464*c4762a1bSJed Brown be = PetscPowScalar(ae,bm1); 465*c4762a1bSJed Brown /* de = be * ae; */ 466*c4762a1bSJed Brown de = PetscPowScalar(ae,beta); 467*c4762a1bSJed Brown ge = coef*be*(te - t0); 468*c4762a1bSJed Brown 469*c4762a1bSJed Brown /* right-hand bottom boundary */ 470*c4762a1bSJed Brown if (j == 0) { 471*c4762a1bSJed Brown 472*c4762a1bSJed Brown tn = x[j+1][i]; 473*c4762a1bSJed Brown an = 0.5*(t0 + tn); 474*c4762a1bSJed Brown bn = PetscPowScalar(an,bm1); 475*c4762a1bSJed Brown /* dn = bn * an; */ 476*c4762a1bSJed Brown dn = PetscPowScalar(an,beta); 477*c4762a1bSJed Brown gn = coef*bn*(tn - t0); 478*c4762a1bSJed Brown 479*c4762a1bSJed Brown v[0] = -hydhx*(dw - gw); col[0].j = j; col[0].i = i-1; 480*c4762a1bSJed Brown v[1] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i; 481*c4762a1bSJed Brown v[2] = -hxdhy*(dn + gn); col[2].j = j+1; col[2].i = i; 482*c4762a1bSJed Brown ierr = MatSetValuesStencil(B,1,&row,3,col,v,INSERT_VALUES);CHKERRQ(ierr); 483*c4762a1bSJed Brown 484*c4762a1bSJed Brown /* right-hand interior boundary */ 485*c4762a1bSJed Brown } else if (j < my-1) { 486*c4762a1bSJed Brown 487*c4762a1bSJed Brown ts = x[j-1][i]; 488*c4762a1bSJed Brown as = 0.5*(t0 + ts); 489*c4762a1bSJed Brown bs = PetscPowScalar(as,bm1); 490*c4762a1bSJed Brown /* ds = bs * as; */ 491*c4762a1bSJed Brown ds = PetscPowScalar(as,beta); 492*c4762a1bSJed Brown gs = coef*bs*(t0 - ts); 493*c4762a1bSJed Brown 494*c4762a1bSJed Brown tn = x[j+1][i]; 495*c4762a1bSJed Brown an = 0.5*(t0 + tn); 496*c4762a1bSJed Brown bn = PetscPowScalar(an,bm1); 497*c4762a1bSJed Brown /* dn = bn * an; */ 498*c4762a1bSJed Brown dn = PetscPowScalar(an,beta); 499*c4762a1bSJed Brown gn = coef*bn*(tn - t0); 500*c4762a1bSJed Brown 501*c4762a1bSJed Brown v[0] = -hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i; 502*c4762a1bSJed Brown v[1] = -hydhx*(dw - gw); col[1].j = j; col[1].i = i-1; 503*c4762a1bSJed Brown v[2] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge); col[2].j = row.j = j; col[2].i = row.i = i; 504*c4762a1bSJed Brown v[3] = -hxdhy*(dn + gn); col[3].j = j+1; col[3].i = i; 505*c4762a1bSJed Brown ierr = MatSetValuesStencil(B,1,&row,4,col,v,INSERT_VALUES);CHKERRQ(ierr); 506*c4762a1bSJed Brown /* right-hand top boundary */ 507*c4762a1bSJed Brown } else { 508*c4762a1bSJed Brown 509*c4762a1bSJed Brown ts = x[j-1][i]; 510*c4762a1bSJed Brown as = 0.5*(t0 + ts); 511*c4762a1bSJed Brown bs = PetscPowScalar(as,bm1); 512*c4762a1bSJed Brown /* ds = bs * as; */ 513*c4762a1bSJed Brown ds = PetscPowScalar(as,beta); 514*c4762a1bSJed Brown gs = coef*bs*(t0 - ts); 515*c4762a1bSJed Brown 516*c4762a1bSJed Brown v[0] = -hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i; 517*c4762a1bSJed Brown v[1] = -hydhx*(dw - gw); col[1].j = j; col[1].i = i-1; 518*c4762a1bSJed Brown v[2] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge); col[2].j = row.j = j; col[2].i = row.i = i; 519*c4762a1bSJed Brown ierr = MatSetValuesStencil(B,1,&row,3,col,v,INSERT_VALUES);CHKERRQ(ierr); 520*c4762a1bSJed Brown } 521*c4762a1bSJed Brown 522*c4762a1bSJed Brown /* bottom boundary,and i <> 0 or mx-1 */ 523*c4762a1bSJed Brown } else if (j == 0) { 524*c4762a1bSJed Brown 525*c4762a1bSJed Brown tw = x[j][i-1]; 526*c4762a1bSJed Brown aw = 0.5*(t0 + tw); 527*c4762a1bSJed Brown bw = PetscPowScalar(aw,bm1); 528*c4762a1bSJed Brown /* dw = bw * aw */ 529*c4762a1bSJed Brown dw = PetscPowScalar(aw,beta); 530*c4762a1bSJed Brown gw = coef*bw*(t0 - tw); 531*c4762a1bSJed Brown 532*c4762a1bSJed Brown te = x[j][i+1]; 533*c4762a1bSJed Brown ae = 0.5*(t0 + te); 534*c4762a1bSJed Brown be = PetscPowScalar(ae,bm1); 535*c4762a1bSJed Brown /* de = be * ae; */ 536*c4762a1bSJed Brown de = PetscPowScalar(ae,beta); 537*c4762a1bSJed Brown ge = coef*be*(te - t0); 538*c4762a1bSJed Brown 539*c4762a1bSJed Brown tn = x[j+1][i]; 540*c4762a1bSJed Brown an = 0.5*(t0 + tn); 541*c4762a1bSJed Brown bn = PetscPowScalar(an,bm1); 542*c4762a1bSJed Brown /* dn = bn * an; */ 543*c4762a1bSJed Brown dn = PetscPowScalar(an,beta); 544*c4762a1bSJed Brown gn = coef*bn*(tn - t0); 545*c4762a1bSJed Brown 546*c4762a1bSJed Brown v[0] = -hydhx*(dw - gw); col[0].j = j; col[0].i = i-1; 547*c4762a1bSJed Brown v[1] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i; 548*c4762a1bSJed Brown v[2] = -hydhx*(de + ge); col[2].j = j; col[2].i = i+1; 549*c4762a1bSJed Brown v[3] = -hxdhy*(dn + gn); col[3].j = j+1; col[3].i = i; 550*c4762a1bSJed Brown ierr = MatSetValuesStencil(B,1,&row,4,col,v,INSERT_VALUES);CHKERRQ(ierr); 551*c4762a1bSJed Brown 552*c4762a1bSJed Brown /* top boundary,and i <> 0 or mx-1 */ 553*c4762a1bSJed Brown } else if (j == my-1) { 554*c4762a1bSJed Brown 555*c4762a1bSJed Brown tw = x[j][i-1]; 556*c4762a1bSJed Brown aw = 0.5*(t0 + tw); 557*c4762a1bSJed Brown bw = PetscPowScalar(aw,bm1); 558*c4762a1bSJed Brown /* dw = bw * aw */ 559*c4762a1bSJed Brown dw = PetscPowScalar(aw,beta); 560*c4762a1bSJed Brown gw = coef*bw*(t0 - tw); 561*c4762a1bSJed Brown 562*c4762a1bSJed Brown te = x[j][i+1]; 563*c4762a1bSJed Brown ae = 0.5*(t0 + te); 564*c4762a1bSJed Brown be = PetscPowScalar(ae,bm1); 565*c4762a1bSJed Brown /* de = be * ae; */ 566*c4762a1bSJed Brown de = PetscPowScalar(ae,beta); 567*c4762a1bSJed Brown ge = coef*be*(te - t0); 568*c4762a1bSJed Brown 569*c4762a1bSJed Brown ts = x[j-1][i]; 570*c4762a1bSJed Brown as = 0.5*(t0 + ts); 571*c4762a1bSJed Brown bs = PetscPowScalar(as,bm1); 572*c4762a1bSJed Brown /* ds = bs * as; */ 573*c4762a1bSJed Brown ds = PetscPowScalar(as,beta); 574*c4762a1bSJed Brown gs = coef*bs*(t0 - ts); 575*c4762a1bSJed Brown 576*c4762a1bSJed Brown v[0] = -hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i; 577*c4762a1bSJed Brown v[1] = -hydhx*(dw - gw); col[1].j = j; col[1].i = i-1; 578*c4762a1bSJed Brown v[2] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge); col[2].j = row.j = j; col[2].i = row.i = i; 579*c4762a1bSJed Brown v[3] = -hydhx*(de + ge); col[3].j = j; col[3].i = i+1; 580*c4762a1bSJed Brown ierr = MatSetValuesStencil(B,1,&row,4,col,v,INSERT_VALUES);CHKERRQ(ierr); 581*c4762a1bSJed Brown 582*c4762a1bSJed Brown } 583*c4762a1bSJed Brown } 584*c4762a1bSJed Brown } 585*c4762a1bSJed Brown ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 586*c4762a1bSJed Brown ierr = DMDAVecRestoreArray(da,localX,&x);CHKERRQ(ierr); 587*c4762a1bSJed Brown ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 588*c4762a1bSJed Brown ierr = DMRestoreLocalVector(da,&localX);CHKERRQ(ierr); 589*c4762a1bSJed Brown if (jac != B) { 590*c4762a1bSJed Brown ierr = MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 591*c4762a1bSJed Brown ierr = MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 592*c4762a1bSJed Brown } 593*c4762a1bSJed Brown 594*c4762a1bSJed Brown ierr = PetscLogFlops((41.0 + 8.0*POWFLOP)*xm*ym);CHKERRQ(ierr); 595*c4762a1bSJed Brown PetscFunctionReturn(0); 596*c4762a1bSJed Brown } 597*c4762a1bSJed Brown 598*c4762a1bSJed Brown 599*c4762a1bSJed Brown 600*c4762a1bSJed Brown /*TEST 601*c4762a1bSJed Brown 602*c4762a1bSJed Brown test: 603*c4762a1bSJed Brown args: -pc_type mg -ksp_type fgmres -da_refine 2 -pc_mg_galerkin pmat -snes_view 604*c4762a1bSJed Brown requires: !single 605*c4762a1bSJed Brown 606*c4762a1bSJed Brown TEST*/ 607