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; 74c4762a1bSJed Brown ierr = PetscOptionsGetReal(NULL,NULL,"-tleft",&user.tleft,NULL);CHKERRQ(ierr); 75c4762a1bSJed Brown ierr = PetscOptionsGetReal(NULL,NULL,"-tright",&user.tright,NULL);CHKERRQ(ierr); 76c4762a1bSJed Brown ierr = PetscOptionsGetReal(NULL,NULL,"-beta",&user.beta,NULL);CHKERRQ(ierr); 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 */ 83c4762a1bSJed Brown ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr); 84c4762a1bSJed Brown 85c4762a1bSJed Brown /* 86c4762a1bSJed Brown Set the DMDA (grid structure) for the grids. 87c4762a1bSJed Brown */ 88c4762a1bSJed 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); 89c4762a1bSJed Brown ierr = DMSetFromOptions(da);CHKERRQ(ierr); 90c4762a1bSJed Brown ierr = DMSetUp(da);CHKERRQ(ierr); 91c4762a1bSJed Brown ierr = DMSetApplicationContext(da,&user);CHKERRQ(ierr); 92c4762a1bSJed Brown ierr = SNESSetDM(snes,(DM)da);CHKERRQ(ierr); 93c4762a1bSJed Brown 94c4762a1bSJed Brown /* 95c4762a1bSJed Brown Create the nonlinear solver, and tell it the functions to use 96c4762a1bSJed Brown */ 97c4762a1bSJed Brown ierr = SNESSetFunction(snes,NULL,FormFunction,&user);CHKERRQ(ierr); 98c4762a1bSJed Brown ierr = SNESSetJacobian(snes,NULL,NULL,FormJacobian,&user);CHKERRQ(ierr); 99c4762a1bSJed Brown ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 100c4762a1bSJed Brown ierr = SNESSetComputeInitialGuess(snes,FormInitialGuess,NULL);CHKERRQ(ierr); 101c4762a1bSJed Brown 102c4762a1bSJed Brown ierr = SNESSolve(snes,NULL,NULL);CHKERRQ(ierr); 103c4762a1bSJed Brown ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr); 104c4762a1bSJed Brown ierr = SNESGetLinearSolveIterations(snes,&lits);CHKERRQ(ierr); 105c4762a1bSJed Brown litspit = ((PetscReal)lits)/((PetscReal)its); 106c4762a1bSJed Brown ierr = PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n",its);CHKERRQ(ierr); 107c4762a1bSJed Brown ierr = PetscPrintf(PETSC_COMM_WORLD,"Number of Linear iterations = %D\n",lits);CHKERRQ(ierr); 108c4762a1bSJed Brown ierr = PetscPrintf(PETSC_COMM_WORLD,"Average Linear its / SNES = %e\n",(double)litspit);CHKERRQ(ierr); 109c4762a1bSJed Brown 110c4762a1bSJed Brown ierr = DMDestroy(&da);CHKERRQ(ierr); 111c4762a1bSJed Brown ierr = SNESDestroy(&snes);CHKERRQ(ierr); 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 PetscErrorCode ierr; 121c4762a1bSJed Brown PetscReal tleft; 122c4762a1bSJed Brown PetscScalar **x; 123c4762a1bSJed Brown DM da; 124c4762a1bSJed Brown 125c4762a1bSJed Brown PetscFunctionBeginUser; 126c4762a1bSJed Brown ierr = SNESGetDM(snes,&da);CHKERRQ(ierr); 127c4762a1bSJed Brown ierr = DMGetApplicationContext(da,&user);CHKERRQ(ierr); 128c4762a1bSJed Brown tleft = user->tleft; 129c4762a1bSJed Brown /* Get ghost points */ 130c4762a1bSJed Brown ierr = DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);CHKERRQ(ierr); 131c4762a1bSJed Brown ierr = DMDAVecGetArray(da,X,&x);CHKERRQ(ierr); 132c4762a1bSJed Brown 133c4762a1bSJed Brown /* Compute initial guess */ 134c4762a1bSJed Brown for (j=ys; j<ys+ym; j++) { 135c4762a1bSJed Brown for (i=xs; i<xs+xm; i++) { 136c4762a1bSJed Brown x[j][i] = tleft; 137c4762a1bSJed Brown } 138c4762a1bSJed Brown } 139c4762a1bSJed Brown ierr = DMDAVecRestoreArray(da,X,&x);CHKERRQ(ierr); 140c4762a1bSJed Brown PetscFunctionReturn(0); 141c4762a1bSJed Brown } 142c4762a1bSJed Brown /* -------------------- Evaluate Function F(x) --------------------- */ 143c4762a1bSJed Brown PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ptr) 144c4762a1bSJed Brown { 145c4762a1bSJed Brown AppCtx *user = (AppCtx*)ptr; 146c4762a1bSJed Brown PetscErrorCode ierr; 147c4762a1bSJed Brown PetscInt i,j,mx,my,xs,ys,xm,ym; 148c4762a1bSJed Brown PetscScalar zero = 0.0,one = 1.0; 149c4762a1bSJed Brown PetscScalar hx,hy,hxdhy,hydhx; 150c4762a1bSJed 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; 151c4762a1bSJed Brown PetscScalar tleft,tright,beta; 152c4762a1bSJed Brown PetscScalar **x,**f; 153c4762a1bSJed Brown Vec localX; 154c4762a1bSJed Brown DM da; 155c4762a1bSJed Brown 156c4762a1bSJed Brown PetscFunctionBeginUser; 157c4762a1bSJed Brown ierr = SNESGetDM(snes,&da);CHKERRQ(ierr); 158c4762a1bSJed Brown ierr = DMGetLocalVector(da,&localX);CHKERRQ(ierr); 159c4762a1bSJed Brown ierr = DMDAGetInfo(da,NULL,&mx,&my,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 160c4762a1bSJed Brown hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1); 161c4762a1bSJed Brown hxdhy = hx/hy; hydhx = hy/hx; 162c4762a1bSJed Brown tleft = user->tleft; tright = user->tright; 163c4762a1bSJed Brown beta = user->beta; 164c4762a1bSJed Brown 165c4762a1bSJed Brown /* Get ghost points */ 166c4762a1bSJed Brown ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);CHKERRQ(ierr); 167c4762a1bSJed Brown ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);CHKERRQ(ierr); 168c4762a1bSJed Brown ierr = DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);CHKERRQ(ierr); 169c4762a1bSJed Brown ierr = DMDAVecGetArray(da,localX,&x);CHKERRQ(ierr); 170c4762a1bSJed Brown ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr); 171c4762a1bSJed Brown 172c4762a1bSJed Brown /* Evaluate function */ 173c4762a1bSJed Brown for (j=ys; j<ys+ym; j++) { 174c4762a1bSJed Brown for (i=xs; i<xs+xm; i++) { 175c4762a1bSJed Brown t0 = x[j][i]; 176c4762a1bSJed Brown 177c4762a1bSJed Brown if (i > 0 && i < mx-1 && j > 0 && j < my-1) { 178c4762a1bSJed Brown 179c4762a1bSJed Brown /* general interior volume */ 180c4762a1bSJed Brown 181c4762a1bSJed Brown tw = x[j][i-1]; 182c4762a1bSJed Brown aw = 0.5*(t0 + tw); 183c4762a1bSJed Brown dw = PetscPowScalar(aw,beta); 184c4762a1bSJed Brown fw = dw*(t0 - tw); 185c4762a1bSJed Brown 186c4762a1bSJed Brown te = x[j][i+1]; 187c4762a1bSJed Brown ae = 0.5*(t0 + te); 188c4762a1bSJed Brown de = PetscPowScalar(ae,beta); 189c4762a1bSJed Brown fe = de*(te - t0); 190c4762a1bSJed Brown 191c4762a1bSJed Brown ts = x[j-1][i]; 192c4762a1bSJed Brown as = 0.5*(t0 + ts); 193c4762a1bSJed Brown ds = PetscPowScalar(as,beta); 194c4762a1bSJed Brown fs = ds*(t0 - ts); 195c4762a1bSJed Brown 196c4762a1bSJed Brown tn = x[j+1][i]; 197c4762a1bSJed Brown an = 0.5*(t0 + tn); 198c4762a1bSJed Brown dn = PetscPowScalar(an,beta); 199c4762a1bSJed Brown fn = dn*(tn - t0); 200c4762a1bSJed Brown 201c4762a1bSJed Brown } else if (i == 0) { 202c4762a1bSJed Brown 203c4762a1bSJed Brown /* left-hand boundary */ 204c4762a1bSJed Brown tw = tleft; 205c4762a1bSJed Brown aw = 0.5*(t0 + tw); 206c4762a1bSJed Brown dw = PetscPowScalar(aw,beta); 207c4762a1bSJed Brown fw = dw*(t0 - tw); 208c4762a1bSJed Brown 209c4762a1bSJed Brown te = x[j][i+1]; 210c4762a1bSJed Brown ae = 0.5*(t0 + te); 211c4762a1bSJed Brown de = PetscPowScalar(ae,beta); 212c4762a1bSJed Brown fe = de*(te - t0); 213c4762a1bSJed Brown 214c4762a1bSJed Brown if (j > 0) { 215c4762a1bSJed Brown ts = x[j-1][i]; 216c4762a1bSJed Brown as = 0.5*(t0 + ts); 217c4762a1bSJed Brown ds = PetscPowScalar(as,beta); 218c4762a1bSJed Brown fs = ds*(t0 - ts); 219c4762a1bSJed Brown } else fs = zero; 220c4762a1bSJed Brown 221c4762a1bSJed Brown if (j < my-1) { 222c4762a1bSJed Brown tn = x[j+1][i]; 223c4762a1bSJed Brown an = 0.5*(t0 + tn); 224c4762a1bSJed Brown dn = PetscPowScalar(an,beta); 225c4762a1bSJed Brown fn = dn*(tn - t0); 226c4762a1bSJed Brown } else fn = zero; 227c4762a1bSJed Brown 228c4762a1bSJed Brown } else if (i == mx-1) { 229c4762a1bSJed Brown 230c4762a1bSJed Brown /* right-hand boundary */ 231c4762a1bSJed Brown tw = x[j][i-1]; 232c4762a1bSJed Brown aw = 0.5*(t0 + tw); 233c4762a1bSJed Brown dw = PetscPowScalar(aw,beta); 234c4762a1bSJed Brown fw = dw*(t0 - tw); 235c4762a1bSJed Brown 236c4762a1bSJed Brown te = tright; 237c4762a1bSJed Brown ae = 0.5*(t0 + te); 238c4762a1bSJed Brown de = PetscPowScalar(ae,beta); 239c4762a1bSJed Brown fe = de*(te - t0); 240c4762a1bSJed Brown 241c4762a1bSJed Brown if (j > 0) { 242c4762a1bSJed Brown ts = x[j-1][i]; 243c4762a1bSJed Brown as = 0.5*(t0 + ts); 244c4762a1bSJed Brown ds = PetscPowScalar(as,beta); 245c4762a1bSJed Brown fs = ds*(t0 - ts); 246c4762a1bSJed Brown } else fs = zero; 247c4762a1bSJed Brown 248c4762a1bSJed Brown if (j < my-1) { 249c4762a1bSJed Brown tn = x[j+1][i]; 250c4762a1bSJed Brown an = 0.5*(t0 + tn); 251c4762a1bSJed Brown dn = PetscPowScalar(an,beta); 252c4762a1bSJed Brown fn = dn*(tn - t0); 253c4762a1bSJed Brown } else fn = zero; 254c4762a1bSJed Brown 255c4762a1bSJed Brown } else if (j == 0) { 256c4762a1bSJed Brown 257c4762a1bSJed Brown /* bottom boundary,and i <> 0 or mx-1 */ 258c4762a1bSJed Brown tw = x[j][i-1]; 259c4762a1bSJed Brown aw = 0.5*(t0 + tw); 260c4762a1bSJed Brown dw = PetscPowScalar(aw,beta); 261c4762a1bSJed Brown fw = dw*(t0 - tw); 262c4762a1bSJed Brown 263c4762a1bSJed Brown te = x[j][i+1]; 264c4762a1bSJed Brown ae = 0.5*(t0 + te); 265c4762a1bSJed Brown de = PetscPowScalar(ae,beta); 266c4762a1bSJed Brown fe = de*(te - t0); 267c4762a1bSJed Brown 268c4762a1bSJed Brown fs = zero; 269c4762a1bSJed Brown 270c4762a1bSJed Brown tn = x[j+1][i]; 271c4762a1bSJed Brown an = 0.5*(t0 + tn); 272c4762a1bSJed Brown dn = PetscPowScalar(an,beta); 273c4762a1bSJed Brown fn = dn*(tn - t0); 274c4762a1bSJed Brown 275c4762a1bSJed Brown } else if (j == my-1) { 276c4762a1bSJed Brown 277c4762a1bSJed Brown /* top boundary,and i <> 0 or mx-1 */ 278c4762a1bSJed Brown tw = x[j][i-1]; 279c4762a1bSJed Brown aw = 0.5*(t0 + tw); 280c4762a1bSJed Brown dw = PetscPowScalar(aw,beta); 281c4762a1bSJed Brown fw = dw*(t0 - tw); 282c4762a1bSJed Brown 283c4762a1bSJed Brown te = x[j][i+1]; 284c4762a1bSJed Brown ae = 0.5*(t0 + te); 285c4762a1bSJed Brown de = PetscPowScalar(ae,beta); 286c4762a1bSJed Brown fe = de*(te - t0); 287c4762a1bSJed Brown 288c4762a1bSJed Brown ts = x[j-1][i]; 289c4762a1bSJed Brown as = 0.5*(t0 + ts); 290c4762a1bSJed Brown ds = PetscPowScalar(as,beta); 291c4762a1bSJed Brown fs = ds*(t0 - ts); 292c4762a1bSJed Brown 293c4762a1bSJed Brown fn = zero; 294c4762a1bSJed Brown 295c4762a1bSJed Brown } 296c4762a1bSJed Brown 297c4762a1bSJed Brown f[j][i] = -hydhx*(fe-fw) - hxdhy*(fn-fs); 298c4762a1bSJed Brown 299c4762a1bSJed Brown } 300c4762a1bSJed Brown } 301c4762a1bSJed Brown ierr = DMDAVecRestoreArray(da,localX,&x);CHKERRQ(ierr); 302c4762a1bSJed Brown ierr = DMDAVecRestoreArray(da,F,&f);CHKERRQ(ierr); 303c4762a1bSJed Brown ierr = DMRestoreLocalVector(da,&localX);CHKERRQ(ierr); 304c4762a1bSJed Brown ierr = PetscLogFlops((22.0 + 4.0*POWFLOP)*ym*xm);CHKERRQ(ierr); 305c4762a1bSJed Brown PetscFunctionReturn(0); 306c4762a1bSJed Brown } 307c4762a1bSJed Brown /* -------------------- Evaluate Jacobian F(x) --------------------- */ 308c4762a1bSJed Brown PetscErrorCode FormJacobian(SNES snes,Vec X,Mat jac,Mat B,void *ptr) 309c4762a1bSJed Brown { 310c4762a1bSJed Brown AppCtx *user = (AppCtx*)ptr; 311c4762a1bSJed Brown PetscErrorCode ierr; 312c4762a1bSJed Brown PetscInt i,j,mx,my,xs,ys,xm,ym; 313c4762a1bSJed Brown PetscScalar one = 1.0,hx,hy,hxdhy,hydhx,t0,tn,ts,te,tw; 314c4762a1bSJed Brown PetscScalar dn,ds,de,dw,an,as,ae,aw,bn,bs,be,bw,gn,gs,ge,gw; 315c4762a1bSJed Brown PetscScalar tleft,tright,beta,bm1,coef; 316c4762a1bSJed Brown PetscScalar v[5],**x; 317c4762a1bSJed Brown Vec localX; 318c4762a1bSJed Brown MatStencil col[5],row; 319c4762a1bSJed Brown DM da; 320c4762a1bSJed Brown 321c4762a1bSJed Brown PetscFunctionBeginUser; 322c4762a1bSJed Brown ierr = SNESGetDM(snes,&da);CHKERRQ(ierr); 323c4762a1bSJed Brown ierr = DMGetLocalVector(da,&localX);CHKERRQ(ierr); 324c4762a1bSJed Brown ierr = DMDAGetInfo(da,NULL,&mx,&my,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 325c4762a1bSJed Brown hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1); 326c4762a1bSJed Brown hxdhy = hx/hy; hydhx = hy/hx; 327c4762a1bSJed Brown tleft = user->tleft; tright = user->tright; 328c4762a1bSJed Brown beta = user->beta; bm1 = user->bm1; coef = user->coef; 329c4762a1bSJed Brown 330c4762a1bSJed Brown /* Get ghost points */ 331c4762a1bSJed Brown ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);CHKERRQ(ierr); 332c4762a1bSJed Brown ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);CHKERRQ(ierr); 333c4762a1bSJed Brown ierr = DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);CHKERRQ(ierr); 334c4762a1bSJed Brown ierr = DMDAVecGetArray(da,localX,&x);CHKERRQ(ierr); 335c4762a1bSJed Brown 336c4762a1bSJed Brown /* Evaluate Jacobian of function */ 337c4762a1bSJed Brown for (j=ys; j<ys+ym; j++) { 338c4762a1bSJed Brown for (i=xs; i<xs+xm; i++) { 339c4762a1bSJed Brown t0 = x[j][i]; 340c4762a1bSJed Brown 341c4762a1bSJed Brown if (i > 0 && i < mx-1 && j > 0 && j < my-1) { 342c4762a1bSJed Brown 343c4762a1bSJed Brown /* general interior volume */ 344c4762a1bSJed Brown 345c4762a1bSJed Brown tw = x[j][i-1]; 346c4762a1bSJed Brown aw = 0.5*(t0 + tw); 347c4762a1bSJed Brown bw = PetscPowScalar(aw,bm1); 348c4762a1bSJed Brown /* dw = bw * aw */ 349c4762a1bSJed Brown dw = PetscPowScalar(aw,beta); 350c4762a1bSJed Brown gw = coef*bw*(t0 - tw); 351c4762a1bSJed Brown 352c4762a1bSJed Brown te = x[j][i+1]; 353c4762a1bSJed Brown ae = 0.5*(t0 + te); 354c4762a1bSJed Brown be = PetscPowScalar(ae,bm1); 355c4762a1bSJed Brown /* de = be * ae; */ 356c4762a1bSJed Brown de = PetscPowScalar(ae,beta); 357c4762a1bSJed Brown ge = coef*be*(te - t0); 358c4762a1bSJed Brown 359c4762a1bSJed Brown ts = x[j-1][i]; 360c4762a1bSJed Brown as = 0.5*(t0 + ts); 361c4762a1bSJed Brown bs = PetscPowScalar(as,bm1); 362c4762a1bSJed Brown /* ds = bs * as; */ 363c4762a1bSJed Brown ds = PetscPowScalar(as,beta); 364c4762a1bSJed Brown gs = coef*bs*(t0 - ts); 365c4762a1bSJed Brown 366c4762a1bSJed Brown tn = x[j+1][i]; 367c4762a1bSJed Brown an = 0.5*(t0 + tn); 368c4762a1bSJed Brown bn = PetscPowScalar(an,bm1); 369c4762a1bSJed Brown /* dn = bn * an; */ 370c4762a1bSJed Brown dn = PetscPowScalar(an,beta); 371c4762a1bSJed Brown gn = coef*bn*(tn - t0); 372c4762a1bSJed Brown 373c4762a1bSJed Brown v[0] = -hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i; 374c4762a1bSJed Brown v[1] = -hydhx*(dw - gw); col[1].j = j; col[1].i = i-1; 375c4762a1bSJed 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; 376c4762a1bSJed Brown v[3] = -hydhx*(de + ge); col[3].j = j; col[3].i = i+1; 377c4762a1bSJed Brown v[4] = -hxdhy*(dn + gn); col[4].j = j+1; col[4].i = i; 378c4762a1bSJed Brown ierr = MatSetValuesStencil(B,1,&row,5,col,v,INSERT_VALUES);CHKERRQ(ierr); 379c4762a1bSJed Brown 380c4762a1bSJed Brown } else if (i == 0) { 381c4762a1bSJed Brown 382c4762a1bSJed Brown /* left-hand boundary */ 383c4762a1bSJed Brown tw = tleft; 384c4762a1bSJed Brown aw = 0.5*(t0 + tw); 385c4762a1bSJed Brown bw = PetscPowScalar(aw,bm1); 386c4762a1bSJed Brown /* dw = bw * aw */ 387c4762a1bSJed Brown dw = PetscPowScalar(aw,beta); 388c4762a1bSJed Brown gw = coef*bw*(t0 - tw); 389c4762a1bSJed Brown 390c4762a1bSJed Brown te = x[j][i + 1]; 391c4762a1bSJed Brown ae = 0.5*(t0 + te); 392c4762a1bSJed Brown be = PetscPowScalar(ae,bm1); 393c4762a1bSJed Brown /* de = be * ae; */ 394c4762a1bSJed Brown de = PetscPowScalar(ae,beta); 395c4762a1bSJed Brown ge = coef*be*(te - t0); 396c4762a1bSJed Brown 397c4762a1bSJed Brown /* left-hand bottom boundary */ 398c4762a1bSJed Brown if (j == 0) { 399c4762a1bSJed Brown 400c4762a1bSJed Brown tn = x[j+1][i]; 401c4762a1bSJed Brown an = 0.5*(t0 + tn); 402c4762a1bSJed Brown bn = PetscPowScalar(an,bm1); 403c4762a1bSJed Brown /* dn = bn * an; */ 404c4762a1bSJed Brown dn = PetscPowScalar(an,beta); 405c4762a1bSJed Brown gn = coef*bn*(tn - t0); 406c4762a1bSJed Brown 407c4762a1bSJed Brown v[0] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[0].j = row.j = j; col[0].i = row.i = i; 408c4762a1bSJed Brown v[1] = -hydhx*(de + ge); col[1].j = j; col[1].i = i+1; 409c4762a1bSJed Brown v[2] = -hxdhy*(dn + gn); col[2].j = j+1; col[2].i = i; 410c4762a1bSJed Brown ierr = MatSetValuesStencil(B,1,&row,3,col,v,INSERT_VALUES);CHKERRQ(ierr); 411c4762a1bSJed Brown 412c4762a1bSJed Brown /* left-hand interior boundary */ 413c4762a1bSJed Brown } else if (j < my-1) { 414c4762a1bSJed Brown 415c4762a1bSJed Brown ts = x[j-1][i]; 416c4762a1bSJed Brown as = 0.5*(t0 + ts); 417c4762a1bSJed Brown bs = PetscPowScalar(as,bm1); 418c4762a1bSJed Brown /* ds = bs * as; */ 419c4762a1bSJed Brown ds = PetscPowScalar(as,beta); 420c4762a1bSJed Brown gs = coef*bs*(t0 - ts); 421c4762a1bSJed Brown 422c4762a1bSJed Brown tn = x[j+1][i]; 423c4762a1bSJed Brown an = 0.5*(t0 + tn); 424c4762a1bSJed Brown bn = PetscPowScalar(an,bm1); 425c4762a1bSJed Brown /* dn = bn * an; */ 426c4762a1bSJed Brown dn = PetscPowScalar(an,beta); 427c4762a1bSJed Brown gn = coef*bn*(tn - t0); 428c4762a1bSJed Brown 429c4762a1bSJed Brown v[0] = -hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i; 430c4762a1bSJed 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; 431c4762a1bSJed Brown v[2] = -hydhx*(de + ge); col[2].j = j; col[2].i = i+1; 432c4762a1bSJed Brown v[3] = -hxdhy*(dn + gn); col[3].j = j+1; col[3].i = i; 433c4762a1bSJed Brown ierr = MatSetValuesStencil(B,1,&row,4,col,v,INSERT_VALUES);CHKERRQ(ierr); 434c4762a1bSJed Brown /* left-hand top boundary */ 435c4762a1bSJed Brown } else { 436c4762a1bSJed Brown 437c4762a1bSJed Brown ts = x[j-1][i]; 438c4762a1bSJed Brown as = 0.5*(t0 + ts); 439c4762a1bSJed Brown bs = PetscPowScalar(as,bm1); 440c4762a1bSJed Brown /* ds = bs * as; */ 441c4762a1bSJed Brown ds = PetscPowScalar(as,beta); 442c4762a1bSJed Brown gs = coef*bs*(t0 - ts); 443c4762a1bSJed Brown 444c4762a1bSJed Brown v[0] = -hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i; 445c4762a1bSJed Brown v[1] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i; 446c4762a1bSJed Brown v[2] = -hydhx*(de + ge); col[2].j = j; col[2].i = i+1; 447c4762a1bSJed Brown ierr = MatSetValuesStencil(B,1,&row,3,col,v,INSERT_VALUES);CHKERRQ(ierr); 448c4762a1bSJed Brown } 449c4762a1bSJed Brown 450c4762a1bSJed Brown } else if (i == mx-1) { 451c4762a1bSJed Brown 452c4762a1bSJed Brown /* right-hand boundary */ 453c4762a1bSJed Brown tw = x[j][i-1]; 454c4762a1bSJed Brown aw = 0.5*(t0 + tw); 455c4762a1bSJed Brown bw = PetscPowScalar(aw,bm1); 456c4762a1bSJed Brown /* dw = bw * aw */ 457c4762a1bSJed Brown dw = PetscPowScalar(aw,beta); 458c4762a1bSJed Brown gw = coef*bw*(t0 - tw); 459c4762a1bSJed Brown 460c4762a1bSJed Brown te = tright; 461c4762a1bSJed Brown ae = 0.5*(t0 + te); 462c4762a1bSJed Brown be = PetscPowScalar(ae,bm1); 463c4762a1bSJed Brown /* de = be * ae; */ 464c4762a1bSJed Brown de = PetscPowScalar(ae,beta); 465c4762a1bSJed Brown ge = coef*be*(te - t0); 466c4762a1bSJed Brown 467c4762a1bSJed Brown /* right-hand bottom boundary */ 468c4762a1bSJed Brown if (j == 0) { 469c4762a1bSJed Brown 470c4762a1bSJed Brown tn = x[j+1][i]; 471c4762a1bSJed Brown an = 0.5*(t0 + tn); 472c4762a1bSJed Brown bn = PetscPowScalar(an,bm1); 473c4762a1bSJed Brown /* dn = bn * an; */ 474c4762a1bSJed Brown dn = PetscPowScalar(an,beta); 475c4762a1bSJed Brown gn = coef*bn*(tn - t0); 476c4762a1bSJed Brown 477c4762a1bSJed Brown v[0] = -hydhx*(dw - gw); col[0].j = j; col[0].i = i-1; 478c4762a1bSJed Brown v[1] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i; 479c4762a1bSJed Brown v[2] = -hxdhy*(dn + gn); col[2].j = j+1; col[2].i = i; 480c4762a1bSJed Brown ierr = MatSetValuesStencil(B,1,&row,3,col,v,INSERT_VALUES);CHKERRQ(ierr); 481c4762a1bSJed Brown 482c4762a1bSJed Brown /* right-hand interior boundary */ 483c4762a1bSJed Brown } else if (j < my-1) { 484c4762a1bSJed Brown 485c4762a1bSJed Brown ts = x[j-1][i]; 486c4762a1bSJed Brown as = 0.5*(t0 + ts); 487c4762a1bSJed Brown bs = PetscPowScalar(as,bm1); 488c4762a1bSJed Brown /* ds = bs * as; */ 489c4762a1bSJed Brown ds = PetscPowScalar(as,beta); 490c4762a1bSJed Brown gs = coef*bs*(t0 - ts); 491c4762a1bSJed Brown 492c4762a1bSJed Brown tn = x[j+1][i]; 493c4762a1bSJed Brown an = 0.5*(t0 + tn); 494c4762a1bSJed Brown bn = PetscPowScalar(an,bm1); 495c4762a1bSJed Brown /* dn = bn * an; */ 496c4762a1bSJed Brown dn = PetscPowScalar(an,beta); 497c4762a1bSJed Brown gn = coef*bn*(tn - t0); 498c4762a1bSJed Brown 499c4762a1bSJed Brown v[0] = -hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i; 500c4762a1bSJed Brown v[1] = -hydhx*(dw - gw); col[1].j = j; col[1].i = i-1; 501c4762a1bSJed 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; 502c4762a1bSJed Brown v[3] = -hxdhy*(dn + gn); col[3].j = j+1; col[3].i = i; 503c4762a1bSJed Brown ierr = MatSetValuesStencil(B,1,&row,4,col,v,INSERT_VALUES);CHKERRQ(ierr); 504c4762a1bSJed Brown /* right-hand top boundary */ 505c4762a1bSJed Brown } else { 506c4762a1bSJed Brown 507c4762a1bSJed Brown ts = x[j-1][i]; 508c4762a1bSJed Brown as = 0.5*(t0 + ts); 509c4762a1bSJed Brown bs = PetscPowScalar(as,bm1); 510c4762a1bSJed Brown /* ds = bs * as; */ 511c4762a1bSJed Brown ds = PetscPowScalar(as,beta); 512c4762a1bSJed Brown gs = coef*bs*(t0 - ts); 513c4762a1bSJed Brown 514c4762a1bSJed Brown v[0] = -hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i; 515c4762a1bSJed Brown v[1] = -hydhx*(dw - gw); col[1].j = j; col[1].i = i-1; 516c4762a1bSJed Brown v[2] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge); col[2].j = row.j = j; col[2].i = row.i = i; 517c4762a1bSJed Brown ierr = MatSetValuesStencil(B,1,&row,3,col,v,INSERT_VALUES);CHKERRQ(ierr); 518c4762a1bSJed Brown } 519c4762a1bSJed Brown 520c4762a1bSJed Brown /* bottom boundary,and i <> 0 or mx-1 */ 521c4762a1bSJed Brown } else if (j == 0) { 522c4762a1bSJed Brown 523c4762a1bSJed Brown tw = x[j][i-1]; 524c4762a1bSJed Brown aw = 0.5*(t0 + tw); 525c4762a1bSJed Brown bw = PetscPowScalar(aw,bm1); 526c4762a1bSJed Brown /* dw = bw * aw */ 527c4762a1bSJed Brown dw = PetscPowScalar(aw,beta); 528c4762a1bSJed Brown gw = coef*bw*(t0 - tw); 529c4762a1bSJed Brown 530c4762a1bSJed Brown te = x[j][i+1]; 531c4762a1bSJed Brown ae = 0.5*(t0 + te); 532c4762a1bSJed Brown be = PetscPowScalar(ae,bm1); 533c4762a1bSJed Brown /* de = be * ae; */ 534c4762a1bSJed Brown de = PetscPowScalar(ae,beta); 535c4762a1bSJed Brown ge = coef*be*(te - t0); 536c4762a1bSJed Brown 537c4762a1bSJed Brown tn = x[j+1][i]; 538c4762a1bSJed Brown an = 0.5*(t0 + tn); 539c4762a1bSJed Brown bn = PetscPowScalar(an,bm1); 540c4762a1bSJed Brown /* dn = bn * an; */ 541c4762a1bSJed Brown dn = PetscPowScalar(an,beta); 542c4762a1bSJed Brown gn = coef*bn*(tn - t0); 543c4762a1bSJed Brown 544c4762a1bSJed Brown v[0] = -hydhx*(dw - gw); col[0].j = j; col[0].i = i-1; 545c4762a1bSJed Brown v[1] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i; 546c4762a1bSJed Brown v[2] = -hydhx*(de + ge); col[2].j = j; col[2].i = i+1; 547c4762a1bSJed Brown v[3] = -hxdhy*(dn + gn); col[3].j = j+1; col[3].i = i; 548c4762a1bSJed Brown ierr = MatSetValuesStencil(B,1,&row,4,col,v,INSERT_VALUES);CHKERRQ(ierr); 549c4762a1bSJed Brown 550c4762a1bSJed Brown /* top boundary,and i <> 0 or mx-1 */ 551c4762a1bSJed Brown } else if (j == my-1) { 552c4762a1bSJed Brown 553c4762a1bSJed Brown tw = x[j][i-1]; 554c4762a1bSJed Brown aw = 0.5*(t0 + tw); 555c4762a1bSJed Brown bw = PetscPowScalar(aw,bm1); 556c4762a1bSJed Brown /* dw = bw * aw */ 557c4762a1bSJed Brown dw = PetscPowScalar(aw,beta); 558c4762a1bSJed Brown gw = coef*bw*(t0 - tw); 559c4762a1bSJed Brown 560c4762a1bSJed Brown te = x[j][i+1]; 561c4762a1bSJed Brown ae = 0.5*(t0 + te); 562c4762a1bSJed Brown be = PetscPowScalar(ae,bm1); 563c4762a1bSJed Brown /* de = be * ae; */ 564c4762a1bSJed Brown de = PetscPowScalar(ae,beta); 565c4762a1bSJed Brown ge = coef*be*(te - t0); 566c4762a1bSJed Brown 567c4762a1bSJed Brown ts = x[j-1][i]; 568c4762a1bSJed Brown as = 0.5*(t0 + ts); 569c4762a1bSJed Brown bs = PetscPowScalar(as,bm1); 570c4762a1bSJed Brown /* ds = bs * as; */ 571c4762a1bSJed Brown ds = PetscPowScalar(as,beta); 572c4762a1bSJed Brown gs = coef*bs*(t0 - ts); 573c4762a1bSJed Brown 574c4762a1bSJed Brown v[0] = -hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i; 575c4762a1bSJed Brown v[1] = -hydhx*(dw - gw); col[1].j = j; col[1].i = i-1; 576c4762a1bSJed Brown v[2] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge); col[2].j = row.j = j; col[2].i = row.i = i; 577c4762a1bSJed Brown v[3] = -hydhx*(de + ge); col[3].j = j; col[3].i = i+1; 578c4762a1bSJed Brown ierr = MatSetValuesStencil(B,1,&row,4,col,v,INSERT_VALUES);CHKERRQ(ierr); 579c4762a1bSJed Brown 580c4762a1bSJed Brown } 581c4762a1bSJed Brown } 582c4762a1bSJed Brown } 583c4762a1bSJed Brown ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 584c4762a1bSJed Brown ierr = DMDAVecRestoreArray(da,localX,&x);CHKERRQ(ierr); 585c4762a1bSJed Brown ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 586c4762a1bSJed Brown ierr = DMRestoreLocalVector(da,&localX);CHKERRQ(ierr); 587c4762a1bSJed Brown if (jac != B) { 588c4762a1bSJed Brown ierr = MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 589c4762a1bSJed Brown ierr = MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 590c4762a1bSJed Brown } 591c4762a1bSJed Brown 592c4762a1bSJed Brown ierr = PetscLogFlops((41.0 + 8.0*POWFLOP)*xm*ym);CHKERRQ(ierr); 593c4762a1bSJed Brown PetscFunctionReturn(0); 594c4762a1bSJed Brown } 595c4762a1bSJed Brown 596c4762a1bSJed Brown /*TEST 597c4762a1bSJed Brown 598c4762a1bSJed Brown test: 599*41ba4c6cSHeeho Park suffix: 1 600c4762a1bSJed Brown args: -pc_type mg -ksp_type fgmres -da_refine 2 -pc_mg_galerkin pmat -snes_view 601c4762a1bSJed Brown requires: !single 602c4762a1bSJed Brown 603*41ba4c6cSHeeho Park test: 604*41ba4c6cSHeeho Park suffix: 2 605*41ba4c6cSHeeho Park args: -pc_type mg -ksp_type fgmres -da_refine 2 -pc_mg_galerkin pmat -snes_view -snes_type newtontrdc 606*41ba4c6cSHeeho Park requires: !single 607*41ba4c6cSHeeho Park 608c4762a1bSJed Brown TEST*/ 609