1 static char help[] = "Solves a time-dependent linear PDE with discontinuous right hand side.\n"; 2 3 /* ------------------------------------------------------------------------ 4 5 This program solves the one-dimensional quench front problem modeling a cooled 6 liquid rising on a hot metal rod 7 u_t = u_xx + g(u), 8 with 9 g(u) = -Au if u <= u_c, 10 = 0 if u > u_c 11 on the domain 0 <= x <= 1, with the boundary conditions 12 u(t,0) = 0, u_x(t,1) = 0, 13 and the initial condition 14 u(0,x) = 0 if 0 <= x <= 0.1, 15 = (x - 0.1)/0.15 if 0.1 < x < 0.25 16 = 1 if 0.25 <= x <= 1 17 We discretize the right-hand side using finite differences with 18 uniform grid spacing h: 19 u_xx = (u_{i+1} - 2u_{i} + u_{i-1})/(h^2) 20 21 Reference: L. Shampine and S. Thompson, "Event Location for Ordinary Differential Equations", 22 http://www.radford.edu/~thompson/webddes/eventsweb.pdf 23 ------------------------------------------------------------------------- */ 24 25 #include <petscdmda.h> 26 #include <petscts.h> 27 /* 28 User-defined application context - contains data needed by the 29 application-provided call-back routines. 30 */ 31 typedef struct { 32 PetscReal A; 33 PetscReal uc; 34 PetscInt *sw; 35 } AppCtx; 36 37 PetscErrorCode InitialConditions(Vec U, DM da, AppCtx *app) 38 { 39 Vec xcoord; 40 PetscScalar *x, *u; 41 PetscInt lsize, M, xs, xm, i; 42 43 PetscFunctionBeginUser; 44 PetscCall(DMGetCoordinates(da, &xcoord)); 45 PetscCall(DMDAVecGetArrayRead(da, xcoord, &x)); 46 47 PetscCall(VecGetLocalSize(U, &lsize)); 48 PetscCall(PetscMalloc1(lsize, &app->sw)); 49 50 PetscCall(DMDAVecGetArray(da, U, &u)); 51 52 PetscCall(DMDAGetInfo(da, 0, &M, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)); 53 PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0)); 54 55 for (i = xs; i < xs + xm; i++) { 56 if (x[i] <= 0.1) u[i] = 0.; 57 else if (x[i] > 0.1 && x[i] < 0.25) u[i] = (x[i] - 0.1) / 0.15; 58 else u[i] = 1.0; 59 60 app->sw[i - xs] = 1; 61 } 62 PetscCall(DMDAVecRestoreArray(da, U, &u)); 63 PetscCall(DMDAVecRestoreArrayRead(da, xcoord, &x)); 64 PetscFunctionReturn(PETSC_SUCCESS); 65 } 66 67 PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscReal *fvalue, void *ctx) 68 { 69 AppCtx *app = (AppCtx *)ctx; 70 const PetscScalar *u; 71 PetscInt i, lsize; 72 73 PetscFunctionBeginUser; 74 PetscCall(VecGetLocalSize(U, &lsize)); 75 PetscCall(VecGetArrayRead(U, &u)); 76 for (i = 0; i < lsize; i++) fvalue[i] = PetscRealPart(u[i]) - app->uc; 77 PetscCall(VecRestoreArrayRead(U, &u)); 78 PetscFunctionReturn(PETSC_SUCCESS); 79 } 80 81 PetscErrorCode PostEventFunction(TS ts, PetscInt nevents_zero, PetscInt events_zero[], PetscReal t, Vec U, PetscBool forwardsolve, void *ctx) 82 { 83 AppCtx *app = (AppCtx *)ctx; 84 PetscInt i, idx; 85 86 PetscFunctionBeginUser; 87 for (i = 0; i < nevents_zero; i++) { 88 idx = events_zero[i]; 89 app->sw[idx] = 0; 90 } 91 PetscFunctionReturn(PETSC_SUCCESS); 92 } 93 94 /* 95 Defines the ODE passed to the ODE solver 96 */ 97 static PetscErrorCode IFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, void *ctx) 98 { 99 AppCtx *app = (AppCtx *)ctx; 100 PetscScalar *f; 101 const PetscScalar *u, *udot; 102 DM da; 103 PetscInt M, xs, xm, i; 104 PetscReal h, h2; 105 Vec Ulocal; 106 107 PetscFunctionBeginUser; 108 PetscCall(TSGetDM(ts, &da)); 109 110 PetscCall(DMDAGetInfo(da, 0, &M, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)); 111 PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0)); 112 113 PetscCall(DMGetLocalVector(da, &Ulocal)); 114 PetscCall(DMGlobalToLocalBegin(da, U, INSERT_VALUES, Ulocal)); 115 PetscCall(DMGlobalToLocalEnd(da, U, INSERT_VALUES, Ulocal)); 116 117 h = 1.0 / (M - 1); 118 h2 = h * h; 119 PetscCall(DMDAVecGetArrayRead(da, Udot, &udot)); 120 PetscCall(DMDAVecGetArrayRead(da, Ulocal, &u)); 121 PetscCall(DMDAVecGetArray(da, F, &f)); 122 123 for (i = xs; i < xs + xm; i++) { 124 if (i == 0) { 125 f[i] = u[i]; 126 } else if (i == M - 1) { 127 f[i] = (u[i] - u[i - 1]) / h; 128 } else { 129 f[i] = (u[i + 1] - 2 * u[i] + u[i - 1]) / h2 + app->sw[i - xs] * (-app->A * u[i]) - udot[i]; 130 } 131 } 132 133 PetscCall(DMDAVecRestoreArrayRead(da, Udot, &udot)); 134 PetscCall(DMDAVecRestoreArrayRead(da, Ulocal, &u)); 135 PetscCall(DMDAVecRestoreArray(da, F, &f)); 136 PetscCall(DMRestoreLocalVector(da, &Ulocal)); 137 138 PetscFunctionReturn(PETSC_SUCCESS); 139 } 140 141 /* 142 Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian. 143 */ 144 static PetscErrorCode IJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat A, Mat B, void *ctx) 145 { 146 AppCtx *app = (AppCtx *)ctx; 147 DM da; 148 MatStencil row, col[3]; 149 PetscScalar v[3]; 150 PetscInt M, xs, xm, i; 151 PetscReal h, h2; 152 153 PetscFunctionBeginUser; 154 PetscCall(TSGetDM(ts, &da)); 155 156 PetscCall(DMDAGetInfo(da, 0, &M, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)); 157 PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0)); 158 159 h = 1.0 / (M - 1); 160 h2 = h * h; 161 for (i = xs; i < xs + xm; i++) { 162 row.i = i; 163 if (i == 0) { 164 v[0] = 1.0; 165 PetscCall(MatSetValuesStencil(A, 1, &row, 1, &row, v, INSERT_VALUES)); 166 } else if (i == M - 1) { 167 col[0].i = i; 168 v[0] = 1 / h; 169 col[1].i = i - 1; 170 v[1] = -1 / h; 171 PetscCall(MatSetValuesStencil(A, 1, &row, 2, col, v, INSERT_VALUES)); 172 } else { 173 col[0].i = i + 1; 174 v[0] = 1 / h2; 175 col[1].i = i; 176 v[1] = -2 / h2 + app->sw[i - xs] * (-app->A) - a; 177 col[2].i = i - 1; 178 v[2] = 1 / h2; 179 PetscCall(MatSetValuesStencil(A, 1, &row, 3, col, v, INSERT_VALUES)); 180 } 181 } 182 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 183 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 184 PetscFunctionReturn(PETSC_SUCCESS); 185 } 186 187 int main(int argc, char **argv) 188 { 189 TS ts; /* ODE integrator */ 190 Vec U; /* solution will be stored here */ 191 Mat J; /* Jacobian matrix */ 192 PetscInt n = 16; 193 AppCtx app; 194 DM da; 195 196 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 197 Initialize program 198 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 199 PetscFunctionBeginUser; 200 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 201 202 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ex22 options", ""); 203 { 204 app.A = 200000; 205 PetscCall(PetscOptionsReal("-A", "", "", app.A, &app.A, NULL)); 206 app.uc = 0.5; 207 PetscCall(PetscOptionsReal("-uc", "", "", app.uc, &app.uc, NULL)); 208 } 209 PetscOptionsEnd(); 210 211 PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, -n, 1, 1, 0, &da)); 212 PetscCall(DMSetFromOptions(da)); 213 PetscCall(DMSetUp(da)); 214 PetscCall(DMDASetUniformCoordinates(da, 0.0, 1.0, 0, 0, 0, 0)); 215 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 216 Create necessary matrix and vectors 217 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 218 PetscCall(DMCreateMatrix(da, &J)); 219 PetscCall(DMCreateGlobalVector(da, &U)); 220 221 PetscCall(InitialConditions(U, da, &app)); 222 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 223 Create timestepping solver context 224 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 225 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 226 PetscCall(TSSetProblemType(ts, TS_NONLINEAR)); 227 PetscCall(TSSetType(ts, TSROSW)); 228 PetscCall(TSSetIFunction(ts, NULL, (TSIFunctionFn *)IFunction, (void *)&app)); 229 PetscCall(TSSetIJacobian(ts, J, J, (TSIJacobianFn *)IJacobian, (void *)&app)); 230 231 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 232 Set initial conditions 233 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 234 PetscCall(TSSetSolution(ts, U)); 235 236 PetscCall(TSSetDM(ts, da)); 237 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 238 Set solver options 239 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 240 PetscCall(TSSetTimeStep(ts, 0.1)); 241 PetscCall(TSSetMaxTime(ts, 30.0)); 242 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); 243 PetscCall(TSSetFromOptions(ts)); 244 245 PetscInt lsize; 246 PetscCall(VecGetLocalSize(U, &lsize)); 247 PetscInt *direction; 248 PetscBool *terminate; 249 PetscInt i; 250 PetscCall(PetscMalloc1(lsize, &direction)); 251 PetscCall(PetscMalloc1(lsize, &terminate)); 252 for (i = 0; i < lsize; i++) { 253 direction[i] = -1; 254 terminate[i] = PETSC_FALSE; 255 } 256 PetscCall(TSSetEventHandler(ts, lsize, direction, terminate, EventFunction, PostEventFunction, (void *)&app)); 257 258 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 259 Run timestepping solver 260 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 261 PetscCall(TSSolve(ts, U)); 262 263 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 264 Free work space. All PETSc objects should be destroyed when they are no longer needed. 265 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 266 267 PetscCall(MatDestroy(&J)); 268 PetscCall(VecDestroy(&U)); 269 PetscCall(DMDestroy(&da)); 270 PetscCall(TSDestroy(&ts)); 271 PetscCall(PetscFree(direction)); 272 PetscCall(PetscFree(terminate)); 273 274 PetscCall(PetscFree(app.sw)); 275 PetscCall(PetscFinalize()); 276 return 0; 277 } 278