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, PetscCtx 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, PetscCtx 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, PetscCtx 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 PetscFunctionReturn(PETSC_SUCCESS); 138 } 139 140 /* 141 Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian. 142 */ 143 static PetscErrorCode IJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat A, Mat B, PetscCtx ctx) 144 { 145 AppCtx *app = (AppCtx *)ctx; 146 DM da; 147 MatStencil row, col[3]; 148 PetscScalar v[3]; 149 PetscInt M, xs, xm, i; 150 PetscReal h, h2; 151 152 PetscFunctionBeginUser; 153 PetscCall(TSGetDM(ts, &da)); 154 155 PetscCall(DMDAGetInfo(da, 0, &M, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)); 156 PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0)); 157 158 h = 1.0 / (M - 1); 159 h2 = h * h; 160 for (i = xs; i < xs + xm; i++) { 161 row.i = i; 162 if (i == 0) { 163 v[0] = 1.0; 164 PetscCall(MatSetValuesStencil(A, 1, &row, 1, &row, v, INSERT_VALUES)); 165 } else if (i == M - 1) { 166 col[0].i = i; 167 v[0] = 1 / h; 168 col[1].i = i - 1; 169 v[1] = -1 / h; 170 PetscCall(MatSetValuesStencil(A, 1, &row, 2, col, v, INSERT_VALUES)); 171 } else { 172 col[0].i = i + 1; 173 v[0] = 1 / h2; 174 col[1].i = i; 175 v[1] = -2 / h2 + app->sw[i - xs] * (-app->A) - a; 176 col[2].i = i - 1; 177 v[2] = 1 / h2; 178 PetscCall(MatSetValuesStencil(A, 1, &row, 3, col, v, INSERT_VALUES)); 179 } 180 } 181 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 182 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 183 PetscFunctionReturn(PETSC_SUCCESS); 184 } 185 186 int main(int argc, char **argv) 187 { 188 TS ts; /* ODE integrator */ 189 Vec U; /* solution will be stored here */ 190 Mat J; /* Jacobian matrix */ 191 PetscInt n = 16; 192 AppCtx app; 193 DM da; 194 195 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 196 Initialize program 197 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 198 PetscFunctionBeginUser; 199 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 200 201 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ex22 options", ""); 202 { 203 app.A = 200000; 204 PetscCall(PetscOptionsReal("-A", "", "", app.A, &app.A, NULL)); 205 app.uc = 0.5; 206 PetscCall(PetscOptionsReal("-uc", "", "", app.uc, &app.uc, NULL)); 207 } 208 PetscOptionsEnd(); 209 210 PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, -n, 1, 1, 0, &da)); 211 PetscCall(DMSetFromOptions(da)); 212 PetscCall(DMSetUp(da)); 213 PetscCall(DMDASetUniformCoordinates(da, 0.0, 1.0, 0, 0, 0, 0)); 214 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 215 Create necessary matrix and vectors 216 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 217 PetscCall(DMCreateMatrix(da, &J)); 218 PetscCall(DMCreateGlobalVector(da, &U)); 219 220 PetscCall(InitialConditions(U, da, &app)); 221 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 222 Create timestepping solver context 223 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 224 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 225 PetscCall(TSSetProblemType(ts, TS_NONLINEAR)); 226 PetscCall(TSSetType(ts, TSROSW)); 227 PetscCall(TSSetIFunction(ts, NULL, (TSIFunctionFn *)IFunction, (void *)&app)); 228 PetscCall(TSSetIJacobian(ts, J, J, (TSIJacobianFn *)IJacobian, (void *)&app)); 229 230 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 231 Set initial conditions 232 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 233 PetscCall(TSSetSolution(ts, U)); 234 235 PetscCall(TSSetDM(ts, da)); 236 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 237 Set solver options 238 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 239 PetscCall(TSSetTimeStep(ts, 0.1)); 240 PetscCall(TSSetMaxTime(ts, 30.0)); 241 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); 242 PetscCall(TSSetFromOptions(ts)); 243 244 PetscInt lsize; 245 PetscCall(VecGetLocalSize(U, &lsize)); 246 PetscInt *direction; 247 PetscBool *terminate; 248 PetscInt i; 249 PetscCall(PetscMalloc1(lsize, &direction)); 250 PetscCall(PetscMalloc1(lsize, &terminate)); 251 for (i = 0; i < lsize; i++) { 252 direction[i] = -1; 253 terminate[i] = PETSC_FALSE; 254 } 255 PetscCall(TSSetEventHandler(ts, lsize, direction, terminate, EventFunction, PostEventFunction, (void *)&app)); 256 257 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 258 Run timestepping solver 259 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 260 PetscCall(TSSolve(ts, U)); 261 262 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 263 Free work space. All PETSc objects should be destroyed when they are no longer needed. 264 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 265 266 PetscCall(MatDestroy(&J)); 267 PetscCall(VecDestroy(&U)); 268 PetscCall(DMDestroy(&da)); 269 PetscCall(TSDestroy(&ts)); 270 PetscCall(PetscFree(direction)); 271 PetscCall(PetscFree(terminate)); 272 273 PetscCall(PetscFree(app.sw)); 274 PetscCall(PetscFinalize()); 275 return 0; 276 } 277