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