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 /* 15c4762a1bSJed Brown 16c4762a1bSJed Brown This example models the partial differential equation 17c4762a1bSJed Brown 18c4762a1bSJed Brown - Div(alpha* T^beta (GRAD T)) = 0. 19c4762a1bSJed Brown 20c4762a1bSJed Brown where beta = 2.5 and alpha = 1.0 21c4762a1bSJed Brown 22c4762a1bSJed Brown BC: T_left = 1.0, T_right = 0.1, dT/dn_top = dTdn_bottom = 0. 23c4762a1bSJed Brown 24c4762a1bSJed Brown in the unit square, which is uniformly discretized in each of x and 25c4762a1bSJed Brown y in this simple encoding. The degrees of freedom are cell centered. 26c4762a1bSJed Brown 27c4762a1bSJed Brown A finite volume approximation with the usual 5-point stencil 28c4762a1bSJed Brown is used to discretize the boundary value problem to obtain a 29c4762a1bSJed Brown nonlinear system of equations. 30c4762a1bSJed Brown 31c4762a1bSJed Brown This code was contributed by David Keyes 32c4762a1bSJed Brown 33c4762a1bSJed Brown */ 34c4762a1bSJed Brown 35c4762a1bSJed Brown #include <petscsnes.h> 36c4762a1bSJed Brown #include <petscdm.h> 37c4762a1bSJed Brown #include <petscdmda.h> 38c4762a1bSJed Brown 39c4762a1bSJed Brown /* User-defined application context */ 40c4762a1bSJed Brown 41c4762a1bSJed Brown typedef struct { 42c4762a1bSJed Brown PetscReal tleft, tright; /* Dirichlet boundary conditions */ 43c4762a1bSJed Brown PetscReal beta, bm1, coef; /* nonlinear diffusivity parameterizations */ 44c4762a1bSJed Brown } AppCtx; 45c4762a1bSJed Brown 46c4762a1bSJed Brown #define POWFLOP 5 /* assume a pow() takes five flops */ 47c4762a1bSJed Brown 48c4762a1bSJed Brown extern PetscErrorCode FormInitialGuess(SNES, Vec, void *); 49c4762a1bSJed Brown extern PetscErrorCode FormFunction(SNES, Vec, Vec, void *); 50c4762a1bSJed Brown extern PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *); 51c4762a1bSJed Brown 52d71ae5a4SJacob Faibussowitsch int main(int argc, char **argv) 53d71ae5a4SJacob Faibussowitsch { 54c4762a1bSJed Brown SNES snes; 55c4762a1bSJed Brown AppCtx user; 56c4762a1bSJed Brown PetscInt its, lits; 57c4762a1bSJed Brown PetscReal litspit; 58c4762a1bSJed Brown DM da; 59c4762a1bSJed Brown 60327415f7SBarry Smith PetscFunctionBeginUser; 619566063dSJacob Faibussowitsch PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 62c4762a1bSJed Brown 63c4762a1bSJed Brown /* set problem parameters */ 64c4762a1bSJed Brown user.tleft = 1.0; 65c4762a1bSJed Brown user.tright = 0.1; 66c4762a1bSJed Brown user.beta = 2.5; 679566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetReal(NULL, NULL, "-tleft", &user.tleft, NULL)); 689566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetReal(NULL, NULL, "-tright", &user.tright, NULL)); 699566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetReal(NULL, NULL, "-beta", &user.beta, NULL)); 70c4762a1bSJed Brown user.bm1 = user.beta - 1.0; 71c4762a1bSJed Brown user.coef = user.beta / 2.0; 72c4762a1bSJed Brown 73c4762a1bSJed Brown /* 74c4762a1bSJed Brown Create the multilevel DM data structure 75c4762a1bSJed Brown */ 769566063dSJacob Faibussowitsch PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes)); 77c4762a1bSJed Brown 78c4762a1bSJed Brown /* 79c4762a1bSJed Brown Set the DMDA (grid structure) for the grids. 80c4762a1bSJed Brown */ 819566063dSJacob Faibussowitsch PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 5, 5, PETSC_DECIDE, PETSC_DECIDE, 1, 1, 0, 0, &da)); 829566063dSJacob Faibussowitsch PetscCall(DMSetFromOptions(da)); 839566063dSJacob Faibussowitsch PetscCall(DMSetUp(da)); 849566063dSJacob Faibussowitsch PetscCall(DMSetApplicationContext(da, &user)); 859566063dSJacob Faibussowitsch PetscCall(SNESSetDM(snes, (DM)da)); 86c4762a1bSJed Brown 87c4762a1bSJed Brown /* 88c4762a1bSJed Brown Create the nonlinear solver, and tell it the functions to use 89c4762a1bSJed Brown */ 909566063dSJacob Faibussowitsch PetscCall(SNESSetFunction(snes, NULL, FormFunction, &user)); 919566063dSJacob Faibussowitsch PetscCall(SNESSetJacobian(snes, NULL, NULL, FormJacobian, &user)); 929566063dSJacob Faibussowitsch PetscCall(SNESSetFromOptions(snes)); 939566063dSJacob Faibussowitsch PetscCall(SNESSetComputeInitialGuess(snes, FormInitialGuess, NULL)); 94c4762a1bSJed Brown 959566063dSJacob Faibussowitsch PetscCall(SNESSolve(snes, NULL, NULL)); 969566063dSJacob Faibussowitsch PetscCall(SNESGetIterationNumber(snes, &its)); 979566063dSJacob Faibussowitsch PetscCall(SNESGetLinearSolveIterations(snes, &lits)); 98c4762a1bSJed Brown litspit = ((PetscReal)lits) / ((PetscReal)its); 9963a3b9bcSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %" PetscInt_FMT "\n", its)); 10063a3b9bcSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of Linear iterations = %" PetscInt_FMT "\n", lits)); 1019566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Average Linear its / SNES = %e\n", (double)litspit)); 102c4762a1bSJed Brown 1039566063dSJacob Faibussowitsch PetscCall(DMDestroy(&da)); 1049566063dSJacob Faibussowitsch PetscCall(SNESDestroy(&snes)); 1059566063dSJacob Faibussowitsch PetscCall(PetscFinalize()); 106b122ec5aSJacob Faibussowitsch return 0; 107c4762a1bSJed Brown } 108c4762a1bSJed Brown /* -------------------- Form initial approximation ----------------- */ 109d71ae5a4SJacob Faibussowitsch PetscErrorCode FormInitialGuess(SNES snes, Vec X, void *ctx) 110d71ae5a4SJacob Faibussowitsch { 111c4762a1bSJed Brown AppCtx *user; 112c4762a1bSJed Brown PetscInt i, j, xs, ys, xm, ym; 113c4762a1bSJed Brown PetscReal tleft; 114c4762a1bSJed Brown PetscScalar **x; 115c4762a1bSJed Brown DM da; 116c4762a1bSJed Brown 117c4762a1bSJed Brown PetscFunctionBeginUser; 1189566063dSJacob Faibussowitsch PetscCall(SNESGetDM(snes, &da)); 1199566063dSJacob Faibussowitsch PetscCall(DMGetApplicationContext(da, &user)); 120c4762a1bSJed Brown tleft = user->tleft; 121c4762a1bSJed Brown /* Get ghost points */ 1229566063dSJacob Faibussowitsch PetscCall(DMDAGetCorners(da, &xs, &ys, 0, &xm, &ym, 0)); 1239566063dSJacob Faibussowitsch PetscCall(DMDAVecGetArray(da, X, &x)); 124c4762a1bSJed Brown 125c4762a1bSJed Brown /* Compute initial guess */ 126c4762a1bSJed Brown for (j = ys; j < ys + ym; j++) { 127ad540459SPierre Jolivet for (i = xs; i < xs + xm; i++) x[j][i] = tleft; 128c4762a1bSJed Brown } 1299566063dSJacob Faibussowitsch PetscCall(DMDAVecRestoreArray(da, X, &x)); 130*3ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 131c4762a1bSJed Brown } 132c4762a1bSJed Brown /* -------------------- Evaluate Function F(x) --------------------- */ 133d71ae5a4SJacob Faibussowitsch PetscErrorCode FormFunction(SNES snes, Vec X, Vec F, void *ptr) 134d71ae5a4SJacob Faibussowitsch { 135c4762a1bSJed Brown AppCtx *user = (AppCtx *)ptr; 136c4762a1bSJed Brown PetscInt i, j, mx, my, xs, ys, xm, ym; 137c4762a1bSJed Brown PetscScalar zero = 0.0, one = 1.0; 138c4762a1bSJed Brown PetscScalar hx, hy, hxdhy, hydhx; 139c4762a1bSJed 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; 140c4762a1bSJed Brown PetscScalar tleft, tright, beta; 141c4762a1bSJed Brown PetscScalar **x, **f; 142c4762a1bSJed Brown Vec localX; 143c4762a1bSJed Brown DM da; 144c4762a1bSJed Brown 145c4762a1bSJed Brown PetscFunctionBeginUser; 1469566063dSJacob Faibussowitsch PetscCall(SNESGetDM(snes, &da)); 1479566063dSJacob Faibussowitsch PetscCall(DMGetLocalVector(da, &localX)); 1489566063dSJacob Faibussowitsch PetscCall(DMDAGetInfo(da, NULL, &mx, &my, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)); 1499371c9d4SSatish Balay hx = one / (PetscReal)(mx - 1); 1509371c9d4SSatish Balay hy = one / (PetscReal)(my - 1); 1519371c9d4SSatish Balay hxdhy = hx / hy; 1529371c9d4SSatish Balay hydhx = hy / hx; 1539371c9d4SSatish Balay tleft = user->tleft; 1549371c9d4SSatish Balay tright = user->tright; 155c4762a1bSJed Brown beta = user->beta; 156c4762a1bSJed Brown 157c4762a1bSJed Brown /* Get ghost points */ 1589566063dSJacob Faibussowitsch PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX)); 1599566063dSJacob Faibussowitsch PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX)); 1609566063dSJacob Faibussowitsch PetscCall(DMDAGetCorners(da, &xs, &ys, 0, &xm, &ym, 0)); 1619566063dSJacob Faibussowitsch PetscCall(DMDAVecGetArray(da, localX, &x)); 1629566063dSJacob Faibussowitsch PetscCall(DMDAVecGetArray(da, F, &f)); 163c4762a1bSJed Brown 164c4762a1bSJed Brown /* Evaluate function */ 165c4762a1bSJed Brown for (j = ys; j < ys + ym; j++) { 166c4762a1bSJed Brown for (i = xs; i < xs + xm; i++) { 167c4762a1bSJed Brown t0 = x[j][i]; 168c4762a1bSJed Brown 169c4762a1bSJed Brown if (i > 0 && i < mx - 1 && j > 0 && j < my - 1) { 170c4762a1bSJed Brown /* general interior volume */ 171c4762a1bSJed Brown 172c4762a1bSJed Brown tw = x[j][i - 1]; 173c4762a1bSJed Brown aw = 0.5 * (t0 + tw); 174c4762a1bSJed Brown dw = PetscPowScalar(aw, beta); 175c4762a1bSJed Brown fw = dw * (t0 - tw); 176c4762a1bSJed Brown 177c4762a1bSJed Brown te = x[j][i + 1]; 178c4762a1bSJed Brown ae = 0.5 * (t0 + te); 179c4762a1bSJed Brown de = PetscPowScalar(ae, beta); 180c4762a1bSJed Brown fe = de * (te - t0); 181c4762a1bSJed Brown 182c4762a1bSJed Brown ts = x[j - 1][i]; 183c4762a1bSJed Brown as = 0.5 * (t0 + ts); 184c4762a1bSJed Brown ds = PetscPowScalar(as, beta); 185c4762a1bSJed Brown fs = ds * (t0 - ts); 186c4762a1bSJed Brown 187c4762a1bSJed Brown tn = x[j + 1][i]; 188c4762a1bSJed Brown an = 0.5 * (t0 + tn); 189c4762a1bSJed Brown dn = PetscPowScalar(an, beta); 190c4762a1bSJed Brown fn = dn * (tn - t0); 191c4762a1bSJed Brown 192c4762a1bSJed Brown } else if (i == 0) { 193c4762a1bSJed Brown /* left-hand boundary */ 194c4762a1bSJed Brown tw = tleft; 195c4762a1bSJed Brown aw = 0.5 * (t0 + tw); 196c4762a1bSJed Brown dw = PetscPowScalar(aw, beta); 197c4762a1bSJed Brown fw = dw * (t0 - tw); 198c4762a1bSJed Brown 199c4762a1bSJed Brown te = x[j][i + 1]; 200c4762a1bSJed Brown ae = 0.5 * (t0 + te); 201c4762a1bSJed Brown de = PetscPowScalar(ae, beta); 202c4762a1bSJed Brown fe = de * (te - t0); 203c4762a1bSJed Brown 204c4762a1bSJed Brown if (j > 0) { 205c4762a1bSJed Brown ts = x[j - 1][i]; 206c4762a1bSJed Brown as = 0.5 * (t0 + ts); 207c4762a1bSJed Brown ds = PetscPowScalar(as, beta); 208c4762a1bSJed Brown fs = ds * (t0 - ts); 209c4762a1bSJed Brown } else fs = zero; 210c4762a1bSJed Brown 211c4762a1bSJed Brown if (j < my - 1) { 212c4762a1bSJed Brown tn = x[j + 1][i]; 213c4762a1bSJed Brown an = 0.5 * (t0 + tn); 214c4762a1bSJed Brown dn = PetscPowScalar(an, beta); 215c4762a1bSJed Brown fn = dn * (tn - t0); 216c4762a1bSJed Brown } else fn = zero; 217c4762a1bSJed Brown 218c4762a1bSJed Brown } else if (i == mx - 1) { 219c4762a1bSJed Brown /* right-hand boundary */ 220c4762a1bSJed Brown tw = x[j][i - 1]; 221c4762a1bSJed Brown aw = 0.5 * (t0 + tw); 222c4762a1bSJed Brown dw = PetscPowScalar(aw, beta); 223c4762a1bSJed Brown fw = dw * (t0 - tw); 224c4762a1bSJed Brown 225c4762a1bSJed Brown te = tright; 226c4762a1bSJed Brown ae = 0.5 * (t0 + te); 227c4762a1bSJed Brown de = PetscPowScalar(ae, beta); 228c4762a1bSJed Brown fe = de * (te - t0); 229c4762a1bSJed Brown 230c4762a1bSJed Brown if (j > 0) { 231c4762a1bSJed Brown ts = x[j - 1][i]; 232c4762a1bSJed Brown as = 0.5 * (t0 + ts); 233c4762a1bSJed Brown ds = PetscPowScalar(as, beta); 234c4762a1bSJed Brown fs = ds * (t0 - ts); 235c4762a1bSJed Brown } else fs = zero; 236c4762a1bSJed Brown 237c4762a1bSJed Brown if (j < my - 1) { 238c4762a1bSJed Brown tn = x[j + 1][i]; 239c4762a1bSJed Brown an = 0.5 * (t0 + tn); 240c4762a1bSJed Brown dn = PetscPowScalar(an, beta); 241c4762a1bSJed Brown fn = dn * (tn - t0); 242c4762a1bSJed Brown } else fn = zero; 243c4762a1bSJed Brown 244c4762a1bSJed Brown } else if (j == 0) { 245c4762a1bSJed Brown /* bottom boundary,and i <> 0 or mx-1 */ 246c4762a1bSJed Brown tw = x[j][i - 1]; 247c4762a1bSJed Brown aw = 0.5 * (t0 + tw); 248c4762a1bSJed Brown dw = PetscPowScalar(aw, beta); 249c4762a1bSJed Brown fw = dw * (t0 - tw); 250c4762a1bSJed Brown 251c4762a1bSJed Brown te = x[j][i + 1]; 252c4762a1bSJed Brown ae = 0.5 * (t0 + te); 253c4762a1bSJed Brown de = PetscPowScalar(ae, beta); 254c4762a1bSJed Brown fe = de * (te - t0); 255c4762a1bSJed Brown 256c4762a1bSJed Brown fs = zero; 257c4762a1bSJed Brown 258c4762a1bSJed Brown tn = x[j + 1][i]; 259c4762a1bSJed Brown an = 0.5 * (t0 + tn); 260c4762a1bSJed Brown dn = PetscPowScalar(an, beta); 261c4762a1bSJed Brown fn = dn * (tn - t0); 262c4762a1bSJed Brown 263c4762a1bSJed Brown } else if (j == my - 1) { 264c4762a1bSJed Brown /* top boundary,and i <> 0 or mx-1 */ 265c4762a1bSJed Brown tw = x[j][i - 1]; 266c4762a1bSJed Brown aw = 0.5 * (t0 + tw); 267c4762a1bSJed Brown dw = PetscPowScalar(aw, beta); 268c4762a1bSJed Brown fw = dw * (t0 - tw); 269c4762a1bSJed Brown 270c4762a1bSJed Brown te = x[j][i + 1]; 271c4762a1bSJed Brown ae = 0.5 * (t0 + te); 272c4762a1bSJed Brown de = PetscPowScalar(ae, beta); 273c4762a1bSJed Brown fe = de * (te - t0); 274c4762a1bSJed Brown 275c4762a1bSJed Brown ts = x[j - 1][i]; 276c4762a1bSJed Brown as = 0.5 * (t0 + ts); 277c4762a1bSJed Brown ds = PetscPowScalar(as, beta); 278c4762a1bSJed Brown fs = ds * (t0 - ts); 279c4762a1bSJed Brown 280c4762a1bSJed Brown fn = zero; 281c4762a1bSJed Brown } 282c4762a1bSJed Brown 283c4762a1bSJed Brown f[j][i] = -hydhx * (fe - fw) - hxdhy * (fn - fs); 284c4762a1bSJed Brown } 285c4762a1bSJed Brown } 2869566063dSJacob Faibussowitsch PetscCall(DMDAVecRestoreArray(da, localX, &x)); 2879566063dSJacob Faibussowitsch PetscCall(DMDAVecRestoreArray(da, F, &f)); 2889566063dSJacob Faibussowitsch PetscCall(DMRestoreLocalVector(da, &localX)); 2899566063dSJacob Faibussowitsch PetscCall(PetscLogFlops((22.0 + 4.0 * POWFLOP) * ym * xm)); 290*3ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 291c4762a1bSJed Brown } 292c4762a1bSJed Brown /* -------------------- Evaluate Jacobian F(x) --------------------- */ 293d71ae5a4SJacob Faibussowitsch PetscErrorCode FormJacobian(SNES snes, Vec X, Mat jac, Mat B, void *ptr) 294d71ae5a4SJacob Faibussowitsch { 295c4762a1bSJed Brown AppCtx *user = (AppCtx *)ptr; 296c4762a1bSJed Brown PetscInt i, j, mx, my, xs, ys, xm, ym; 297c4762a1bSJed Brown PetscScalar one = 1.0, hx, hy, hxdhy, hydhx, t0, tn, ts, te, tw; 298c4762a1bSJed Brown PetscScalar dn, ds, de, dw, an, as, ae, aw, bn, bs, be, bw, gn, gs, ge, gw; 299c4762a1bSJed Brown PetscScalar tleft, tright, beta, bm1, coef; 300c4762a1bSJed Brown PetscScalar v[5], **x; 301c4762a1bSJed Brown Vec localX; 302c4762a1bSJed Brown MatStencil col[5], row; 303c4762a1bSJed Brown DM da; 304c4762a1bSJed Brown 305c4762a1bSJed Brown PetscFunctionBeginUser; 3069566063dSJacob Faibussowitsch PetscCall(SNESGetDM(snes, &da)); 3079566063dSJacob Faibussowitsch PetscCall(DMGetLocalVector(da, &localX)); 3089566063dSJacob Faibussowitsch PetscCall(DMDAGetInfo(da, NULL, &mx, &my, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)); 3099371c9d4SSatish Balay hx = one / (PetscReal)(mx - 1); 3109371c9d4SSatish Balay hy = one / (PetscReal)(my - 1); 3119371c9d4SSatish Balay hxdhy = hx / hy; 3129371c9d4SSatish Balay hydhx = hy / hx; 3139371c9d4SSatish Balay tleft = user->tleft; 3149371c9d4SSatish Balay tright = user->tright; 3159371c9d4SSatish Balay beta = user->beta; 3169371c9d4SSatish Balay bm1 = user->bm1; 3179371c9d4SSatish Balay coef = user->coef; 318c4762a1bSJed Brown 319c4762a1bSJed Brown /* Get ghost points */ 3209566063dSJacob Faibussowitsch PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX)); 3219566063dSJacob Faibussowitsch PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX)); 3229566063dSJacob Faibussowitsch PetscCall(DMDAGetCorners(da, &xs, &ys, 0, &xm, &ym, 0)); 3239566063dSJacob Faibussowitsch PetscCall(DMDAVecGetArray(da, localX, &x)); 324c4762a1bSJed Brown 325c4762a1bSJed Brown /* Evaluate Jacobian of function */ 326c4762a1bSJed Brown for (j = ys; j < ys + ym; j++) { 327c4762a1bSJed Brown for (i = xs; i < xs + xm; i++) { 328c4762a1bSJed Brown t0 = x[j][i]; 329c4762a1bSJed Brown 330c4762a1bSJed Brown if (i > 0 && i < mx - 1 && j > 0 && j < my - 1) { 331c4762a1bSJed Brown /* general interior volume */ 332c4762a1bSJed Brown 333c4762a1bSJed Brown tw = x[j][i - 1]; 334c4762a1bSJed Brown aw = 0.5 * (t0 + tw); 335c4762a1bSJed Brown bw = PetscPowScalar(aw, bm1); 336c4762a1bSJed Brown /* dw = bw * aw */ 337c4762a1bSJed Brown dw = PetscPowScalar(aw, beta); 338c4762a1bSJed Brown gw = coef * bw * (t0 - tw); 339c4762a1bSJed Brown 340c4762a1bSJed Brown te = x[j][i + 1]; 341c4762a1bSJed Brown ae = 0.5 * (t0 + te); 342c4762a1bSJed Brown be = PetscPowScalar(ae, bm1); 343c4762a1bSJed Brown /* de = be * ae; */ 344c4762a1bSJed Brown de = PetscPowScalar(ae, beta); 345c4762a1bSJed Brown ge = coef * be * (te - t0); 346c4762a1bSJed Brown 347c4762a1bSJed Brown ts = x[j - 1][i]; 348c4762a1bSJed Brown as = 0.5 * (t0 + ts); 349c4762a1bSJed Brown bs = PetscPowScalar(as, bm1); 350c4762a1bSJed Brown /* ds = bs * as; */ 351c4762a1bSJed Brown ds = PetscPowScalar(as, beta); 352c4762a1bSJed Brown gs = coef * bs * (t0 - ts); 353c4762a1bSJed Brown 354c4762a1bSJed Brown tn = x[j + 1][i]; 355c4762a1bSJed Brown an = 0.5 * (t0 + tn); 356c4762a1bSJed Brown bn = PetscPowScalar(an, bm1); 357c4762a1bSJed Brown /* dn = bn * an; */ 358c4762a1bSJed Brown dn = PetscPowScalar(an, beta); 359c4762a1bSJed Brown gn = coef * bn * (tn - t0); 360c4762a1bSJed Brown 3619371c9d4SSatish Balay v[0] = -hxdhy * (ds - gs); 3629371c9d4SSatish Balay col[0].j = j - 1; 3639371c9d4SSatish Balay col[0].i = i; 3649371c9d4SSatish Balay v[1] = -hydhx * (dw - gw); 3659371c9d4SSatish Balay col[1].j = j; 3669371c9d4SSatish Balay col[1].i = i - 1; 3679371c9d4SSatish Balay v[2] = hxdhy * (ds + dn + gs - gn) + hydhx * (dw + de + gw - ge); 3689371c9d4SSatish Balay col[2].j = row.j = j; 3699371c9d4SSatish Balay col[2].i = row.i = i; 3709371c9d4SSatish Balay v[3] = -hydhx * (de + ge); 3719371c9d4SSatish Balay col[3].j = j; 3729371c9d4SSatish Balay col[3].i = i + 1; 3739371c9d4SSatish Balay v[4] = -hxdhy * (dn + gn); 3749371c9d4SSatish Balay col[4].j = j + 1; 3759371c9d4SSatish Balay col[4].i = i; 3769566063dSJacob Faibussowitsch PetscCall(MatSetValuesStencil(B, 1, &row, 5, col, v, INSERT_VALUES)); 377c4762a1bSJed Brown 378c4762a1bSJed Brown } else if (i == 0) { 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 tn = x[j + 1][i]; 397c4762a1bSJed Brown an = 0.5 * (t0 + tn); 398c4762a1bSJed Brown bn = PetscPowScalar(an, bm1); 399c4762a1bSJed Brown /* dn = bn * an; */ 400c4762a1bSJed Brown dn = PetscPowScalar(an, beta); 401c4762a1bSJed Brown gn = coef * bn * (tn - t0); 402c4762a1bSJed Brown 4039371c9d4SSatish Balay v[0] = hxdhy * (dn - gn) + hydhx * (dw + de + gw - ge); 4049371c9d4SSatish Balay col[0].j = row.j = j; 4059371c9d4SSatish Balay col[0].i = row.i = i; 4069371c9d4SSatish Balay v[1] = -hydhx * (de + ge); 4079371c9d4SSatish Balay col[1].j = j; 4089371c9d4SSatish Balay col[1].i = i + 1; 4099371c9d4SSatish Balay v[2] = -hxdhy * (dn + gn); 4109371c9d4SSatish Balay col[2].j = j + 1; 4119371c9d4SSatish Balay col[2].i = i; 4129566063dSJacob Faibussowitsch PetscCall(MatSetValuesStencil(B, 1, &row, 3, col, v, INSERT_VALUES)); 413c4762a1bSJed Brown 414c4762a1bSJed Brown /* left-hand interior boundary */ 415c4762a1bSJed Brown } else if (j < my - 1) { 416c4762a1bSJed Brown ts = x[j - 1][i]; 417c4762a1bSJed Brown as = 0.5 * (t0 + ts); 418c4762a1bSJed Brown bs = PetscPowScalar(as, bm1); 419c4762a1bSJed Brown /* ds = bs * as; */ 420c4762a1bSJed Brown ds = PetscPowScalar(as, beta); 421c4762a1bSJed Brown gs = coef * bs * (t0 - ts); 422c4762a1bSJed Brown 423c4762a1bSJed Brown tn = x[j + 1][i]; 424c4762a1bSJed Brown an = 0.5 * (t0 + tn); 425c4762a1bSJed Brown bn = PetscPowScalar(an, bm1); 426c4762a1bSJed Brown /* dn = bn * an; */ 427c4762a1bSJed Brown dn = PetscPowScalar(an, beta); 428c4762a1bSJed Brown gn = coef * bn * (tn - t0); 429c4762a1bSJed Brown 4309371c9d4SSatish Balay v[0] = -hxdhy * (ds - gs); 4319371c9d4SSatish Balay col[0].j = j - 1; 4329371c9d4SSatish Balay col[0].i = i; 4339371c9d4SSatish Balay v[1] = hxdhy * (ds + dn + gs - gn) + hydhx * (dw + de + gw - ge); 4349371c9d4SSatish Balay col[1].j = row.j = j; 4359371c9d4SSatish Balay col[1].i = row.i = i; 4369371c9d4SSatish Balay v[2] = -hydhx * (de + ge); 4379371c9d4SSatish Balay col[2].j = j; 4389371c9d4SSatish Balay col[2].i = i + 1; 4399371c9d4SSatish Balay v[3] = -hxdhy * (dn + gn); 4409371c9d4SSatish Balay col[3].j = j + 1; 4419371c9d4SSatish Balay col[3].i = i; 4429566063dSJacob Faibussowitsch PetscCall(MatSetValuesStencil(B, 1, &row, 4, col, v, INSERT_VALUES)); 443c4762a1bSJed Brown /* left-hand top boundary */ 444c4762a1bSJed Brown } else { 445c4762a1bSJed Brown ts = x[j - 1][i]; 446c4762a1bSJed Brown as = 0.5 * (t0 + ts); 447c4762a1bSJed Brown bs = PetscPowScalar(as, bm1); 448c4762a1bSJed Brown /* ds = bs * as; */ 449c4762a1bSJed Brown ds = PetscPowScalar(as, beta); 450c4762a1bSJed Brown gs = coef * bs * (t0 - ts); 451c4762a1bSJed Brown 4529371c9d4SSatish Balay v[0] = -hxdhy * (ds - gs); 4539371c9d4SSatish Balay col[0].j = j - 1; 4549371c9d4SSatish Balay col[0].i = i; 4559371c9d4SSatish Balay v[1] = hxdhy * (ds + gs) + hydhx * (dw + de + gw - ge); 4569371c9d4SSatish Balay col[1].j = row.j = j; 4579371c9d4SSatish Balay col[1].i = row.i = i; 4589371c9d4SSatish Balay v[2] = -hydhx * (de + ge); 4599371c9d4SSatish Balay col[2].j = j; 4609371c9d4SSatish Balay col[2].i = i + 1; 4619566063dSJacob Faibussowitsch PetscCall(MatSetValuesStencil(B, 1, &row, 3, col, v, INSERT_VALUES)); 462c4762a1bSJed Brown } 463c4762a1bSJed Brown 464c4762a1bSJed Brown } else if (i == mx - 1) { 465c4762a1bSJed Brown /* right-hand boundary */ 466c4762a1bSJed Brown tw = x[j][i - 1]; 467c4762a1bSJed Brown aw = 0.5 * (t0 + tw); 468c4762a1bSJed Brown bw = PetscPowScalar(aw, bm1); 469c4762a1bSJed Brown /* dw = bw * aw */ 470c4762a1bSJed Brown dw = PetscPowScalar(aw, beta); 471c4762a1bSJed Brown gw = coef * bw * (t0 - tw); 472c4762a1bSJed Brown 473c4762a1bSJed Brown te = tright; 474c4762a1bSJed Brown ae = 0.5 * (t0 + te); 475c4762a1bSJed Brown be = PetscPowScalar(ae, bm1); 476c4762a1bSJed Brown /* de = be * ae; */ 477c4762a1bSJed Brown de = PetscPowScalar(ae, beta); 478c4762a1bSJed Brown ge = coef * be * (te - t0); 479c4762a1bSJed Brown 480c4762a1bSJed Brown /* right-hand bottom boundary */ 481c4762a1bSJed Brown if (j == 0) { 482c4762a1bSJed Brown tn = x[j + 1][i]; 483c4762a1bSJed Brown an = 0.5 * (t0 + tn); 484c4762a1bSJed Brown bn = PetscPowScalar(an, bm1); 485c4762a1bSJed Brown /* dn = bn * an; */ 486c4762a1bSJed Brown dn = PetscPowScalar(an, beta); 487c4762a1bSJed Brown gn = coef * bn * (tn - t0); 488c4762a1bSJed Brown 4899371c9d4SSatish Balay v[0] = -hydhx * (dw - gw); 4909371c9d4SSatish Balay col[0].j = j; 4919371c9d4SSatish Balay col[0].i = i - 1; 4929371c9d4SSatish Balay v[1] = hxdhy * (dn - gn) + hydhx * (dw + de + gw - ge); 4939371c9d4SSatish Balay col[1].j = row.j = j; 4949371c9d4SSatish Balay col[1].i = row.i = i; 4959371c9d4SSatish Balay v[2] = -hxdhy * (dn + gn); 4969371c9d4SSatish Balay col[2].j = j + 1; 4979371c9d4SSatish Balay col[2].i = i; 4989566063dSJacob Faibussowitsch PetscCall(MatSetValuesStencil(B, 1, &row, 3, col, v, INSERT_VALUES)); 499c4762a1bSJed Brown 500c4762a1bSJed Brown /* right-hand interior boundary */ 501c4762a1bSJed Brown } else if (j < my - 1) { 502c4762a1bSJed Brown ts = x[j - 1][i]; 503c4762a1bSJed Brown as = 0.5 * (t0 + ts); 504c4762a1bSJed Brown bs = PetscPowScalar(as, bm1); 505c4762a1bSJed Brown /* ds = bs * as; */ 506c4762a1bSJed Brown ds = PetscPowScalar(as, beta); 507c4762a1bSJed Brown gs = coef * bs * (t0 - ts); 508c4762a1bSJed Brown 509c4762a1bSJed Brown tn = x[j + 1][i]; 510c4762a1bSJed Brown an = 0.5 * (t0 + tn); 511c4762a1bSJed Brown bn = PetscPowScalar(an, bm1); 512c4762a1bSJed Brown /* dn = bn * an; */ 513c4762a1bSJed Brown dn = PetscPowScalar(an, beta); 514c4762a1bSJed Brown gn = coef * bn * (tn - t0); 515c4762a1bSJed Brown 5169371c9d4SSatish Balay v[0] = -hxdhy * (ds - gs); 5179371c9d4SSatish Balay col[0].j = j - 1; 5189371c9d4SSatish Balay col[0].i = i; 5199371c9d4SSatish Balay v[1] = -hydhx * (dw - gw); 5209371c9d4SSatish Balay col[1].j = j; 5219371c9d4SSatish Balay col[1].i = i - 1; 5229371c9d4SSatish Balay v[2] = hxdhy * (ds + dn + gs - gn) + hydhx * (dw + de + gw - ge); 5239371c9d4SSatish Balay col[2].j = row.j = j; 5249371c9d4SSatish Balay col[2].i = row.i = i; 5259371c9d4SSatish Balay v[3] = -hxdhy * (dn + gn); 5269371c9d4SSatish Balay col[3].j = j + 1; 5279371c9d4SSatish Balay col[3].i = i; 5289566063dSJacob Faibussowitsch PetscCall(MatSetValuesStencil(B, 1, &row, 4, col, v, INSERT_VALUES)); 529c4762a1bSJed Brown /* right-hand top boundary */ 530c4762a1bSJed Brown } else { 531c4762a1bSJed Brown ts = x[j - 1][i]; 532c4762a1bSJed Brown as = 0.5 * (t0 + ts); 533c4762a1bSJed Brown bs = PetscPowScalar(as, bm1); 534c4762a1bSJed Brown /* ds = bs * as; */ 535c4762a1bSJed Brown ds = PetscPowScalar(as, beta); 536c4762a1bSJed Brown gs = coef * bs * (t0 - ts); 537c4762a1bSJed Brown 5389371c9d4SSatish Balay v[0] = -hxdhy * (ds - gs); 5399371c9d4SSatish Balay col[0].j = j - 1; 5409371c9d4SSatish Balay col[0].i = i; 5419371c9d4SSatish Balay v[1] = -hydhx * (dw - gw); 5429371c9d4SSatish Balay col[1].j = j; 5439371c9d4SSatish Balay col[1].i = i - 1; 5449371c9d4SSatish Balay v[2] = hxdhy * (ds + gs) + hydhx * (dw + de + gw - ge); 5459371c9d4SSatish Balay col[2].j = row.j = j; 5469371c9d4SSatish Balay col[2].i = row.i = i; 5479566063dSJacob Faibussowitsch PetscCall(MatSetValuesStencil(B, 1, &row, 3, col, v, INSERT_VALUES)); 548c4762a1bSJed Brown } 549c4762a1bSJed Brown 550c4762a1bSJed Brown /* bottom boundary,and i <> 0 or mx-1 */ 551c4762a1bSJed Brown } else if (j == 0) { 552c4762a1bSJed Brown tw = x[j][i - 1]; 553c4762a1bSJed Brown aw = 0.5 * (t0 + tw); 554c4762a1bSJed Brown bw = PetscPowScalar(aw, bm1); 555c4762a1bSJed Brown /* dw = bw * aw */ 556c4762a1bSJed Brown dw = PetscPowScalar(aw, beta); 557c4762a1bSJed Brown gw = coef * bw * (t0 - tw); 558c4762a1bSJed Brown 559c4762a1bSJed Brown te = x[j][i + 1]; 560c4762a1bSJed Brown ae = 0.5 * (t0 + te); 561c4762a1bSJed Brown be = PetscPowScalar(ae, bm1); 562c4762a1bSJed Brown /* de = be * ae; */ 563c4762a1bSJed Brown de = PetscPowScalar(ae, beta); 564c4762a1bSJed Brown ge = coef * be * (te - t0); 565c4762a1bSJed Brown 566c4762a1bSJed Brown tn = x[j + 1][i]; 567c4762a1bSJed Brown an = 0.5 * (t0 + tn); 568c4762a1bSJed Brown bn = PetscPowScalar(an, bm1); 569c4762a1bSJed Brown /* dn = bn * an; */ 570c4762a1bSJed Brown dn = PetscPowScalar(an, beta); 571c4762a1bSJed Brown gn = coef * bn * (tn - t0); 572c4762a1bSJed Brown 5739371c9d4SSatish Balay v[0] = -hydhx * (dw - gw); 5749371c9d4SSatish Balay col[0].j = j; 5759371c9d4SSatish Balay col[0].i = i - 1; 5769371c9d4SSatish Balay v[1] = hxdhy * (dn - gn) + hydhx * (dw + de + gw - ge); 5779371c9d4SSatish Balay col[1].j = row.j = j; 5789371c9d4SSatish Balay col[1].i = row.i = i; 5799371c9d4SSatish Balay v[2] = -hydhx * (de + ge); 5809371c9d4SSatish Balay col[2].j = j; 5819371c9d4SSatish Balay col[2].i = i + 1; 5829371c9d4SSatish Balay v[3] = -hxdhy * (dn + gn); 5839371c9d4SSatish Balay col[3].j = j + 1; 5849371c9d4SSatish Balay col[3].i = i; 5859566063dSJacob Faibussowitsch PetscCall(MatSetValuesStencil(B, 1, &row, 4, col, v, INSERT_VALUES)); 586c4762a1bSJed Brown 587c4762a1bSJed Brown /* top boundary,and i <> 0 or mx-1 */ 588c4762a1bSJed Brown } else if (j == my - 1) { 589c4762a1bSJed Brown tw = x[j][i - 1]; 590c4762a1bSJed Brown aw = 0.5 * (t0 + tw); 591c4762a1bSJed Brown bw = PetscPowScalar(aw, bm1); 592c4762a1bSJed Brown /* dw = bw * aw */ 593c4762a1bSJed Brown dw = PetscPowScalar(aw, beta); 594c4762a1bSJed Brown gw = coef * bw * (t0 - tw); 595c4762a1bSJed Brown 596c4762a1bSJed Brown te = x[j][i + 1]; 597c4762a1bSJed Brown ae = 0.5 * (t0 + te); 598c4762a1bSJed Brown be = PetscPowScalar(ae, bm1); 599c4762a1bSJed Brown /* de = be * ae; */ 600c4762a1bSJed Brown de = PetscPowScalar(ae, beta); 601c4762a1bSJed Brown ge = coef * be * (te - t0); 602c4762a1bSJed Brown 603c4762a1bSJed Brown ts = x[j - 1][i]; 604c4762a1bSJed Brown as = 0.5 * (t0 + ts); 605c4762a1bSJed Brown bs = PetscPowScalar(as, bm1); 606c4762a1bSJed Brown /* ds = bs * as; */ 607c4762a1bSJed Brown ds = PetscPowScalar(as, beta); 608c4762a1bSJed Brown gs = coef * bs * (t0 - ts); 609c4762a1bSJed Brown 6109371c9d4SSatish Balay v[0] = -hxdhy * (ds - gs); 6119371c9d4SSatish Balay col[0].j = j - 1; 6129371c9d4SSatish Balay col[0].i = i; 6139371c9d4SSatish Balay v[1] = -hydhx * (dw - gw); 6149371c9d4SSatish Balay col[1].j = j; 6159371c9d4SSatish Balay col[1].i = i - 1; 6169371c9d4SSatish Balay v[2] = hxdhy * (ds + gs) + hydhx * (dw + de + gw - ge); 6179371c9d4SSatish Balay col[2].j = row.j = j; 6189371c9d4SSatish Balay col[2].i = row.i = i; 6199371c9d4SSatish Balay v[3] = -hydhx * (de + ge); 6209371c9d4SSatish Balay col[3].j = j; 6219371c9d4SSatish Balay col[3].i = i + 1; 6229566063dSJacob Faibussowitsch PetscCall(MatSetValuesStencil(B, 1, &row, 4, col, v, INSERT_VALUES)); 623c4762a1bSJed Brown } 624c4762a1bSJed Brown } 625c4762a1bSJed Brown } 6269566063dSJacob Faibussowitsch PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 6279566063dSJacob Faibussowitsch PetscCall(DMDAVecRestoreArray(da, localX, &x)); 6289566063dSJacob Faibussowitsch PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 6299566063dSJacob Faibussowitsch PetscCall(DMRestoreLocalVector(da, &localX)); 630c4762a1bSJed Brown if (jac != B) { 6319566063dSJacob Faibussowitsch PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY)); 6329566063dSJacob Faibussowitsch PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY)); 633c4762a1bSJed Brown } 634c4762a1bSJed Brown 6359566063dSJacob Faibussowitsch PetscCall(PetscLogFlops((41.0 + 8.0 * POWFLOP) * xm * ym)); 636*3ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 637c4762a1bSJed Brown } 638c4762a1bSJed Brown 639c4762a1bSJed Brown /*TEST 640c4762a1bSJed Brown 641c4762a1bSJed Brown test: 64241ba4c6cSHeeho Park suffix: 1 643c4762a1bSJed Brown args: -pc_type mg -ksp_type fgmres -da_refine 2 -pc_mg_galerkin pmat -snes_view 644c4762a1bSJed Brown requires: !single 645c4762a1bSJed Brown 64641ba4c6cSHeeho Park test: 64741ba4c6cSHeeho Park suffix: 2 64841ba4c6cSHeeho Park args: -pc_type mg -ksp_type fgmres -da_refine 2 -pc_mg_galerkin pmat -snes_view -snes_type newtontrdc 64941ba4c6cSHeeho Park requires: !single 65041ba4c6cSHeeho Park 651c4762a1bSJed Brown TEST*/ 652