1 /* 2 This example implements the model described in 3 4 Rauenzahn, Mousseau, Knoll. "Temporal accuracy of the nonequilibrium radiation diffusion 5 equations employing a Saha ionization model" 2005. 6 7 The paper discusses three examples, the first two are nondimensional with a simple 8 ionization model. The third example is fully dimensional and uses the Saha ionization 9 model with realistic parameters. 10 */ 11 12 #include <petscts.h> 13 #include <petscdm.h> 14 #include <petscdmda.h> 15 16 typedef enum { 17 BC_DIRICHLET, 18 BC_NEUMANN, 19 BC_ROBIN 20 } BCType; 21 static const char *const BCTypes[] = {"DIRICHLET", "NEUMANN", "ROBIN", "BCType", "BC_", 0}; 22 typedef enum { 23 JACOBIAN_ANALYTIC, 24 JACOBIAN_MATRIXFREE, 25 JACOBIAN_FD_COLORING, 26 JACOBIAN_FD_FULL 27 } JacobianType; 28 static const char *const JacobianTypes[] = {"ANALYTIC", "MATRIXFREE", "FD_COLORING", "FD_FULL", "JacobianType", "FD_", 0}; 29 typedef enum { 30 DISCRETIZATION_FD, 31 DISCRETIZATION_FE 32 } DiscretizationType; 33 static const char *const DiscretizationTypes[] = {"FD", "FE", "DiscretizationType", "DISCRETIZATION_", 0}; 34 typedef enum { 35 QUADRATURE_GAUSS1, 36 QUADRATURE_GAUSS2, 37 QUADRATURE_GAUSS3, 38 QUADRATURE_GAUSS4, 39 QUADRATURE_LOBATTO2, 40 QUADRATURE_LOBATTO3 41 } QuadratureType; 42 static const char *const QuadratureTypes[] = {"GAUSS1", "GAUSS2", "GAUSS3", "GAUSS4", "LOBATTO2", "LOBATTO3", "QuadratureType", "QUADRATURE_", 0}; 43 44 typedef struct { 45 PetscScalar E; /* radiation energy */ 46 PetscScalar T; /* material temperature */ 47 } RDNode; 48 49 typedef struct { 50 PetscReal meter, kilogram, second, Kelvin; /* Fundamental units */ 51 PetscReal Joule, Watt; /* Derived units */ 52 } RDUnit; 53 54 typedef struct _n_RD *RD; 55 56 struct _n_RD { 57 void (*MaterialEnergy)(RD, const RDNode *, PetscScalar *, RDNode *); 58 DM da; 59 PetscBool monitor_residual; 60 DiscretizationType discretization; 61 QuadratureType quadrature; 62 JacobianType jacobian; 63 PetscInt initial; 64 BCType leftbc; 65 PetscBool view_draw; 66 char view_binary[PETSC_MAX_PATH_LEN]; 67 PetscBool test_diff; 68 PetscBool endpoint; 69 PetscBool bclimit; 70 PetscBool bcmidpoint; 71 RDUnit unit; 72 73 /* model constants, see Table 2 and RDCreate() */ 74 PetscReal rho, K_R, K_p, I_H, m_p, m_e, h, k, c, sigma_b, beta, gamma; 75 76 /* Domain and boundary conditions */ 77 PetscReal Eapplied; /* Radiation flux from the left */ 78 PetscReal L; /* Length of domain */ 79 PetscReal final_time; 80 }; 81 82 static PetscErrorCode RDDestroy(RD *rd) { 83 PetscFunctionBeginUser; 84 PetscCall(DMDestroy(&(*rd)->da)); 85 PetscCall(PetscFree(*rd)); 86 PetscFunctionReturn(0); 87 } 88 89 /* The paper has a time derivative for material energy (Eq 2) which is a dependent variable (computable from temperature 90 * and density through an uninvertible relation). Computing this derivative is trivial for trapezoid rule (used in the 91 * paper), but does not generalize nicely to higher order integrators. Here we use the implicit form which provides 92 * time derivatives of the independent variables (radiation energy and temperature), so we must compute the time 93 * derivative of material energy ourselves (could be done using AD). 94 * 95 * There are multiple ionization models, this interface dispatches to the one currently in use. 96 */ 97 static void RDMaterialEnergy(RD rd, const RDNode *n, PetscScalar *Em, RDNode *dEm) { 98 rd->MaterialEnergy(rd, n, Em, dEm); 99 } 100 101 /* Solves a quadratic equation while propagating tangents */ 102 static void QuadraticSolve(PetscScalar a, PetscScalar a_t, PetscScalar b, PetscScalar b_t, PetscScalar c, PetscScalar c_t, PetscScalar *x, PetscScalar *x_t) { 103 PetscScalar disc = b * b - 4. * a * c, disc_t = 2. * b * b_t - 4. * a_t * c - 4. * a * c_t, num = -b + PetscSqrtScalar(disc), /* choose positive sign */ 104 num_t = -b_t + 0.5 / PetscSqrtScalar(disc) * disc_t, den = 2. * a, den_t = 2. * a_t; 105 *x = num / den; 106 *x_t = (num_t * den - num * den_t) / PetscSqr(den); 107 } 108 109 /* The primary model presented in the paper */ 110 static void RDMaterialEnergy_Saha(RD rd, const RDNode *n, PetscScalar *inEm, RDNode *dEm) { 111 PetscScalar Em, alpha, alpha_t, T = n->T, T_t = 1., chi = rd->I_H / (rd->k * T), chi_t = -chi / T * T_t, a = 1., a_t = 0, b = 4. * rd->m_p / rd->rho * PetscPowScalarReal(2. * PETSC_PI * rd->m_e * rd->I_H / PetscSqr(rd->h), 1.5) * PetscExpScalar(-chi) * PetscPowScalarReal(chi, 1.5), /* Eq 7 */ 112 b_t = -b * chi_t + 1.5 * b / chi * chi_t, c = -b, c_t = -b_t; 113 QuadraticSolve(a, a_t, b, b_t, c, c_t, &alpha, &alpha_t); /* Solve Eq 7 for alpha */ 114 Em = rd->k * T / rd->m_p * (1.5 * (1. + alpha) + alpha * chi); /* Eq 6 */ 115 if (inEm) *inEm = Em; 116 if (dEm) { 117 dEm->E = 0; 118 dEm->T = Em / T * T_t + rd->k * T / rd->m_p * (1.5 * alpha_t + alpha_t * chi + alpha * chi_t); 119 } 120 } 121 /* Reduced ionization model, Eq 30 */ 122 static void RDMaterialEnergy_Reduced(RD rd, const RDNode *n, PetscScalar *Em, RDNode *dEm) { 123 PetscScalar alpha, alpha_t, T = n->T, T_t = 1., chi = -0.3 / T, chi_t = -chi / T * T_t, a = 1., a_t = 0., b = PetscExpScalar(chi), b_t = b * chi_t, c = -b, c_t = -b_t; 124 QuadraticSolve(a, a_t, b, b_t, c, c_t, &alpha, &alpha_t); 125 if (Em) *Em = (1. + alpha) * T + 0.3 * alpha; 126 if (dEm) { 127 dEm->E = 0; 128 dEm->T = alpha_t * T + (1. + alpha) * T_t + 0.3 * alpha_t; 129 } 130 } 131 132 /* Eq 5 */ 133 static void RDSigma_R(RD rd, RDNode *n, PetscScalar *sigma_R, RDNode *dsigma_R) { 134 *sigma_R = rd->K_R * rd->rho * PetscPowScalar(n->T, -rd->gamma); 135 dsigma_R->E = 0; 136 dsigma_R->T = -rd->gamma * (*sigma_R) / n->T; 137 } 138 139 /* Eq 4 */ 140 static void RDDiffusionCoefficient(RD rd, PetscBool limit, RDNode *n, RDNode *nx, PetscScalar *D_R, RDNode *dD_R, RDNode *dxD_R) { 141 PetscScalar sigma_R, denom; 142 RDNode dsigma_R, ddenom, dxdenom; 143 144 RDSigma_R(rd, n, &sigma_R, &dsigma_R); 145 denom = 3. * rd->rho * sigma_R + (int)limit * PetscAbsScalar(nx->E) / n->E; 146 ddenom.E = -(int)limit * PetscAbsScalar(nx->E) / PetscSqr(n->E); 147 ddenom.T = 3. * rd->rho * dsigma_R.T; 148 dxdenom.E = (int)limit * (PetscRealPart(nx->E) < 0 ? -1. : 1.) / n->E; 149 dxdenom.T = 0; 150 *D_R = rd->c / denom; 151 if (dD_R) { 152 dD_R->E = -rd->c / PetscSqr(denom) * ddenom.E; 153 dD_R->T = -rd->c / PetscSqr(denom) * ddenom.T; 154 } 155 if (dxD_R) { 156 dxD_R->E = -rd->c / PetscSqr(denom) * dxdenom.E; 157 dxD_R->T = -rd->c / PetscSqr(denom) * dxdenom.T; 158 } 159 } 160 161 static PetscErrorCode RDStateView(RD rd, Vec X, Vec Xdot, Vec F) { 162 DMDALocalInfo info; 163 PetscInt i; 164 const RDNode *x, *xdot, *f; 165 MPI_Comm comm; 166 167 PetscFunctionBeginUser; 168 PetscCall(PetscObjectGetComm((PetscObject)rd->da, &comm)); 169 PetscCall(DMDAGetLocalInfo(rd->da, &info)); 170 PetscCall(DMDAVecGetArrayRead(rd->da, X, (void *)&x)); 171 PetscCall(DMDAVecGetArrayRead(rd->da, Xdot, (void *)&xdot)); 172 PetscCall(DMDAVecGetArrayRead(rd->da, F, (void *)&f)); 173 for (i = info.xs; i < info.xs + info.xm; i++) { 174 PetscCall(PetscSynchronizedPrintf(comm, "x[%" PetscInt_FMT "] (%10.2G,%10.2G) (%10.2G,%10.2G) (%10.2G,%10.2G)\n", i, (double)PetscRealPart(x[i].E), (double)PetscRealPart(x[i].T), (double)PetscRealPart(xdot[i].E), (double)PetscRealPart(xdot[i].T), 175 (double)PetscRealPart(f[i].E), (double)PetscRealPart(f[i].T))); 176 } 177 PetscCall(DMDAVecRestoreArrayRead(rd->da, X, (void *)&x)); 178 PetscCall(DMDAVecRestoreArrayRead(rd->da, Xdot, (void *)&xdot)); 179 PetscCall(DMDAVecRestoreArrayRead(rd->da, F, (void *)&f)); 180 PetscCall(PetscSynchronizedFlush(comm, PETSC_STDOUT)); 181 PetscFunctionReturn(0); 182 } 183 184 static PetscScalar RDRadiation(RD rd, const RDNode *n, RDNode *dn) { 185 PetscScalar sigma_p = rd->K_p * rd->rho * PetscPowScalar(n->T, -rd->beta), sigma_p_T = -rd->beta * sigma_p / n->T, tmp = 4. * rd->sigma_b * PetscSqr(PetscSqr(n->T)) / rd->c - n->E, tmp_E = -1., tmp_T = 4. * rd->sigma_b * 4 * n->T * (PetscSqr(n->T)) / rd->c, rad = sigma_p * rd->c * rd->rho * tmp, rad_E = sigma_p * rd->c * rd->rho * tmp_E, rad_T = rd->c * rd->rho * (sigma_p_T * tmp + sigma_p * tmp_T); 186 if (dn) { 187 dn->E = rad_E; 188 dn->T = rad_T; 189 } 190 return rad; 191 } 192 193 static PetscScalar RDDiffusion(RD rd, PetscReal hx, const RDNode x[], PetscInt i, RDNode d[]) { 194 PetscReal ihx = 1. / hx; 195 RDNode n_L, nx_L, n_R, nx_R, dD_L, dxD_L, dD_R, dxD_R, dfluxL[2], dfluxR[2]; 196 PetscScalar D_L, D_R, fluxL, fluxR; 197 198 n_L.E = 0.5 * (x[i - 1].E + x[i].E); 199 n_L.T = 0.5 * (x[i - 1].T + x[i].T); 200 nx_L.E = (x[i].E - x[i - 1].E) / hx; 201 nx_L.T = (x[i].T - x[i - 1].T) / hx; 202 RDDiffusionCoefficient(rd, PETSC_TRUE, &n_L, &nx_L, &D_L, &dD_L, &dxD_L); 203 fluxL = D_L * nx_L.E; 204 dfluxL[0].E = -ihx * D_L + (0.5 * dD_L.E - ihx * dxD_L.E) * nx_L.E; 205 dfluxL[1].E = +ihx * D_L + (0.5 * dD_L.E + ihx * dxD_L.E) * nx_L.E; 206 dfluxL[0].T = (0.5 * dD_L.T - ihx * dxD_L.T) * nx_L.E; 207 dfluxL[1].T = (0.5 * dD_L.T + ihx * dxD_L.T) * nx_L.E; 208 209 n_R.E = 0.5 * (x[i].E + x[i + 1].E); 210 n_R.T = 0.5 * (x[i].T + x[i + 1].T); 211 nx_R.E = (x[i + 1].E - x[i].E) / hx; 212 nx_R.T = (x[i + 1].T - x[i].T) / hx; 213 RDDiffusionCoefficient(rd, PETSC_TRUE, &n_R, &nx_R, &D_R, &dD_R, &dxD_R); 214 fluxR = D_R * nx_R.E; 215 dfluxR[0].E = -ihx * D_R + (0.5 * dD_R.E - ihx * dxD_R.E) * nx_R.E; 216 dfluxR[1].E = +ihx * D_R + (0.5 * dD_R.E + ihx * dxD_R.E) * nx_R.E; 217 dfluxR[0].T = (0.5 * dD_R.T - ihx * dxD_R.T) * nx_R.E; 218 dfluxR[1].T = (0.5 * dD_R.T + ihx * dxD_R.T) * nx_R.E; 219 220 if (d) { 221 d[0].E = -ihx * dfluxL[0].E; 222 d[0].T = -ihx * dfluxL[0].T; 223 d[1].E = ihx * (dfluxR[0].E - dfluxL[1].E); 224 d[1].T = ihx * (dfluxR[0].T - dfluxL[1].T); 225 d[2].E = ihx * dfluxR[1].E; 226 d[2].T = ihx * dfluxR[1].T; 227 } 228 return ihx * (fluxR - fluxL); 229 } 230 231 static PetscErrorCode RDGetLocalArrays(RD rd, TS ts, Vec X, Vec Xdot, PetscReal *Theta, PetscReal *dt, Vec *X0loc, RDNode **x0, Vec *Xloc, RDNode **x, Vec *Xloc_t, RDNode **xdot) { 232 PetscBool istheta; 233 234 PetscFunctionBeginUser; 235 PetscCall(DMGetLocalVector(rd->da, X0loc)); 236 PetscCall(DMGetLocalVector(rd->da, Xloc)); 237 PetscCall(DMGetLocalVector(rd->da, Xloc_t)); 238 239 PetscCall(DMGlobalToLocalBegin(rd->da, X, INSERT_VALUES, *Xloc)); 240 PetscCall(DMGlobalToLocalEnd(rd->da, X, INSERT_VALUES, *Xloc)); 241 PetscCall(DMGlobalToLocalBegin(rd->da, Xdot, INSERT_VALUES, *Xloc_t)); 242 PetscCall(DMGlobalToLocalEnd(rd->da, Xdot, INSERT_VALUES, *Xloc_t)); 243 244 /* 245 The following is a hack to subvert TSTHETA which is like an implicit midpoint method to behave more like a trapezoid 246 rule. These methods have equivalent linear stability, but the nonlinear stability is somewhat different. The 247 radiation system is inconvenient to write in explicit form because the ionization model is "on the left". 248 */ 249 PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSTHETA, &istheta)); 250 if (istheta && rd->endpoint) { 251 PetscCall(TSThetaGetTheta(ts, Theta)); 252 } else *Theta = 1.; 253 254 PetscCall(TSGetTimeStep(ts, dt)); 255 PetscCall(VecWAXPY(*X0loc, -(*Theta) * (*dt), *Xloc_t, *Xloc)); /* back out the value at the start of this step */ 256 if (rd->endpoint) { PetscCall(VecWAXPY(*Xloc, *dt, *Xloc_t, *X0loc)); /* move the abscissa to the end of the step */ } 257 258 PetscCall(DMDAVecGetArray(rd->da, *X0loc, x0)); 259 PetscCall(DMDAVecGetArray(rd->da, *Xloc, x)); 260 PetscCall(DMDAVecGetArray(rd->da, *Xloc_t, xdot)); 261 PetscFunctionReturn(0); 262 } 263 264 static PetscErrorCode RDRestoreLocalArrays(RD rd, Vec *X0loc, RDNode **x0, Vec *Xloc, RDNode **x, Vec *Xloc_t, RDNode **xdot) { 265 PetscFunctionBeginUser; 266 PetscCall(DMDAVecRestoreArray(rd->da, *X0loc, x0)); 267 PetscCall(DMDAVecRestoreArray(rd->da, *Xloc, x)); 268 PetscCall(DMDAVecRestoreArray(rd->da, *Xloc_t, xdot)); 269 PetscCall(DMRestoreLocalVector(rd->da, X0loc)); 270 PetscCall(DMRestoreLocalVector(rd->da, Xloc)); 271 PetscCall(DMRestoreLocalVector(rd->da, Xloc_t)); 272 PetscFunctionReturn(0); 273 } 274 275 static PetscErrorCode PETSC_UNUSED RDCheckDomain_Private(RD rd, TS ts, Vec X, PetscBool *in) { 276 PetscInt minloc; 277 PetscReal min; 278 279 PetscFunctionBeginUser; 280 PetscCall(VecMin(X, &minloc, &min)); 281 if (min < 0) { 282 SNES snes; 283 *in = PETSC_FALSE; 284 PetscCall(TSGetSNES(ts, &snes)); 285 PetscCall(SNESSetFunctionDomainError(snes)); 286 PetscCall(PetscInfo(ts, "Domain violation at %" PetscInt_FMT " field %" PetscInt_FMT " value %g\n", minloc / 2, minloc % 2, (double)min)); 287 } else *in = PETSC_TRUE; 288 PetscFunctionReturn(0); 289 } 290 291 /* Energy and temperature must remain positive */ 292 #define RDCheckDomain(rd, ts, X) \ 293 do { \ 294 PetscBool _in; \ 295 PetscCall(RDCheckDomain_Private(rd, ts, X, &_in)); \ 296 if (!_in) PetscFunctionReturn(0); \ 297 } while (0) 298 299 static PetscErrorCode RDIFunction_FD(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ctx) { 300 RD rd = (RD)ctx; 301 RDNode *x, *x0, *xdot, *f; 302 Vec X0loc, Xloc, Xloc_t; 303 PetscReal hx, Theta, dt; 304 DMDALocalInfo info; 305 PetscInt i; 306 307 PetscFunctionBeginUser; 308 PetscCall(RDGetLocalArrays(rd, ts, X, Xdot, &Theta, &dt, &X0loc, &x0, &Xloc, &x, &Xloc_t, &xdot)); 309 PetscCall(DMDAVecGetArray(rd->da, F, &f)); 310 PetscCall(DMDAGetLocalInfo(rd->da, &info)); 311 VecZeroEntries(F); 312 313 hx = rd->L / (info.mx - 1); 314 315 for (i = info.xs; i < info.xs + info.xm; i++) { 316 PetscReal rho = rd->rho; 317 PetscScalar Em_t, rad; 318 319 rad = (1. - Theta) * RDRadiation(rd, &x0[i], 0) + Theta * RDRadiation(rd, &x[i], 0); 320 if (rd->endpoint) { 321 PetscScalar Em0, Em1; 322 RDMaterialEnergy(rd, &x0[i], &Em0, NULL); 323 RDMaterialEnergy(rd, &x[i], &Em1, NULL); 324 Em_t = (Em1 - Em0) / dt; 325 } else { 326 RDNode dEm; 327 RDMaterialEnergy(rd, &x[i], NULL, &dEm); 328 Em_t = dEm.E * xdot[i].E + dEm.T * xdot[i].T; 329 } 330 /* Residuals are multiplied by the volume element (hx). */ 331 /* The temperature equation does not have boundary conditions */ 332 f[i].T = hx * (rho * Em_t + rad); 333 334 if (i == 0) { /* Left boundary condition */ 335 PetscScalar D_R, bcTheta = rd->bcmidpoint ? Theta : 1.; 336 RDNode n, nx; 337 338 n.E = (1. - bcTheta) * x0[0].E + bcTheta * x[0].E; 339 n.T = (1. - bcTheta) * x0[0].T + bcTheta * x[0].T; 340 nx.E = ((1. - bcTheta) * (x0[1].E - x0[0].E) + bcTheta * (x[1].E - x[0].E)) / hx; 341 nx.T = ((1. - bcTheta) * (x0[1].T - x0[0].T) + bcTheta * (x[1].T - x[0].T)) / hx; 342 switch (rd->leftbc) { 343 case BC_ROBIN: 344 RDDiffusionCoefficient(rd, rd->bclimit, &n, &nx, &D_R, 0, 0); 345 f[0].E = hx * (n.E - 2. * D_R * nx.E - rd->Eapplied); 346 break; 347 case BC_NEUMANN: f[0].E = x[1].E - x[0].E; break; 348 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Case %" PetscInt_FMT, rd->initial); 349 } 350 } else if (i == info.mx - 1) { /* Right boundary */ 351 f[i].E = x[i].E - x[i - 1].E; /* Homogeneous Neumann */ 352 } else { 353 PetscScalar diff = (1. - Theta) * RDDiffusion(rd, hx, x0, i, 0) + Theta * RDDiffusion(rd, hx, x, i, 0); 354 f[i].E = hx * (xdot[i].E - diff - rad); 355 } 356 } 357 PetscCall(RDRestoreLocalArrays(rd, &X0loc, &x0, &Xloc, &x, &Xloc_t, &xdot)); 358 PetscCall(DMDAVecRestoreArray(rd->da, F, &f)); 359 if (rd->monitor_residual) PetscCall(RDStateView(rd, X, Xdot, F)); 360 PetscFunctionReturn(0); 361 } 362 363 static PetscErrorCode RDIJacobian_FD(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat A, Mat B, void *ctx) { 364 RD rd = (RD)ctx; 365 RDNode *x, *x0, *xdot; 366 Vec X0loc, Xloc, Xloc_t; 367 PetscReal hx, Theta, dt; 368 DMDALocalInfo info; 369 PetscInt i; 370 371 PetscFunctionBeginUser; 372 PetscCall(RDGetLocalArrays(rd, ts, X, Xdot, &Theta, &dt, &X0loc, &x0, &Xloc, &x, &Xloc_t, &xdot)); 373 PetscCall(DMDAGetLocalInfo(rd->da, &info)); 374 hx = rd->L / (info.mx - 1); 375 PetscCall(MatZeroEntries(B)); 376 377 for (i = info.xs; i < info.xs + info.xm; i++) { 378 PetscInt col[3]; 379 PetscReal rho = rd->rho; 380 PetscScalar /*Em_t,rad,*/ K[2][6]; 381 RDNode dEm_t, drad; 382 383 /*rad = (1.-Theta)* */ RDRadiation(rd, &x0[i], 0); /* + Theta* */ 384 RDRadiation(rd, &x[i], &drad); 385 386 if (rd->endpoint) { 387 PetscScalar Em0, Em1; 388 RDNode dEm1; 389 RDMaterialEnergy(rd, &x0[i], &Em0, NULL); 390 RDMaterialEnergy(rd, &x[i], &Em1, &dEm1); 391 /*Em_t = (Em1 - Em0) / (Theta*dt);*/ 392 dEm_t.E = dEm1.E / (Theta * dt); 393 dEm_t.T = dEm1.T / (Theta * dt); 394 } else { 395 const PetscScalar epsilon = x[i].T * PETSC_SQRT_MACHINE_EPSILON; 396 RDNode n1; 397 RDNode dEm, dEm1; 398 PetscScalar Em_TT; 399 400 n1.E = x[i].E; 401 n1.T = x[i].T + epsilon; 402 RDMaterialEnergy(rd, &x[i], NULL, &dEm); 403 RDMaterialEnergy(rd, &n1, NULL, &dEm1); 404 /* The Jacobian needs another derivative. We finite difference here instead of 405 * propagating second derivatives through the ionization model. */ 406 Em_TT = (dEm1.T - dEm.T) / epsilon; 407 /*Em_t = dEm.E * xdot[i].E + dEm.T * xdot[i].T;*/ 408 dEm_t.E = dEm.E * a; 409 dEm_t.T = dEm.T * a + Em_TT * xdot[i].T; 410 } 411 412 PetscCall(PetscMemzero(K, sizeof(K))); 413 /* Residuals are multiplied by the volume element (hx). */ 414 if (i == 0) { 415 PetscScalar D, bcTheta = rd->bcmidpoint ? Theta : 1.; 416 RDNode n, nx; 417 RDNode dD, dxD; 418 419 n.E = (1. - bcTheta) * x0[0].E + bcTheta * x[0].E; 420 n.T = (1. - bcTheta) * x0[0].T + bcTheta * x[0].T; 421 nx.E = ((1. - bcTheta) * (x0[1].E - x0[0].E) + bcTheta * (x[1].E - x[0].E)) / hx; 422 nx.T = ((1. - bcTheta) * (x0[1].T - x0[0].T) + bcTheta * (x[1].T - x[0].T)) / hx; 423 switch (rd->leftbc) { 424 case BC_ROBIN: 425 RDDiffusionCoefficient(rd, rd->bclimit, &n, &nx, &D, &dD, &dxD); 426 K[0][1 * 2 + 0] = (bcTheta / Theta) * hx * (1. - 2. * D * (-1. / hx) - 2. * nx.E * dD.E + 2. * nx.E * dxD.E / hx); 427 K[0][1 * 2 + 1] = (bcTheta / Theta) * hx * (-2. * nx.E * dD.T); 428 K[0][2 * 2 + 0] = (bcTheta / Theta) * hx * (-2. * D * (1. / hx) - 2. * nx.E * dD.E - 2. * nx.E * dxD.E / hx); 429 break; 430 case BC_NEUMANN: 431 K[0][1 * 2 + 0] = -1. / Theta; 432 K[0][2 * 2 + 0] = 1. / Theta; 433 break; 434 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Case %" PetscInt_FMT, rd->initial); 435 } 436 } else if (i == info.mx - 1) { 437 K[0][0 * 2 + 0] = -1. / Theta; 438 K[0][1 * 2 + 0] = 1. / Theta; 439 } else { 440 /*PetscScalar diff;*/ 441 RDNode ddiff[3]; 442 /*diff = (1.-Theta)*RDDiffusion(rd,hx,x0,i,0) + Theta* */ RDDiffusion(rd, hx, x, i, ddiff); 443 K[0][0 * 2 + 0] = -hx * ddiff[0].E; 444 K[0][0 * 2 + 1] = -hx * ddiff[0].T; 445 K[0][1 * 2 + 0] = hx * (a - ddiff[1].E - drad.E); 446 K[0][1 * 2 + 1] = hx * (-ddiff[1].T - drad.T); 447 K[0][2 * 2 + 0] = -hx * ddiff[2].E; 448 K[0][2 * 2 + 1] = -hx * ddiff[2].T; 449 } 450 451 K[1][1 * 2 + 0] = hx * (rho * dEm_t.E + drad.E); 452 K[1][1 * 2 + 1] = hx * (rho * dEm_t.T + drad.T); 453 454 col[0] = i - 1; 455 col[1] = i; 456 col[2] = i + 1 < info.mx ? i + 1 : -1; 457 PetscCall(MatSetValuesBlocked(B, 1, &i, 3, col, &K[0][0], INSERT_VALUES)); 458 } 459 PetscCall(RDRestoreLocalArrays(rd, &X0loc, &x0, &Xloc, &x, &Xloc_t, &xdot)); 460 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 461 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 462 if (A != B) { 463 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 464 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 465 } 466 PetscFunctionReturn(0); 467 } 468 469 /* Evaluate interpolants and derivatives at a select quadrature point */ 470 static void RDEvaluate(PetscReal interp[][2], PetscReal deriv[][2], PetscInt q, const RDNode x[], PetscInt i, RDNode *n, RDNode *nx) { 471 PetscInt j; 472 n->E = 0; 473 n->T = 0; 474 nx->E = 0; 475 nx->T = 0; 476 for (j = 0; j < 2; j++) { 477 n->E += interp[q][j] * x[i + j].E; 478 n->T += interp[q][j] * x[i + j].T; 479 nx->E += deriv[q][j] * x[i + j].E; 480 nx->T += deriv[q][j] * x[i + j].T; 481 } 482 } 483 484 /* 485 Various quadrature rules. The nonlinear terms are non-polynomial so no standard quadrature will be exact. 486 */ 487 static PetscErrorCode RDGetQuadrature(RD rd, PetscReal hx, PetscInt *nq, PetscReal weight[], PetscReal interp[][2], PetscReal deriv[][2]) { 488 PetscInt q, j; 489 const PetscReal *refweight, (*refinterp)[2], (*refderiv)[2]; 490 491 PetscFunctionBeginUser; 492 switch (rd->quadrature) { 493 case QUADRATURE_GAUSS1: { 494 static const PetscReal ww[1] = {1.}; 495 static const PetscReal ii[1][2] = { 496 {0.5, 0.5} 497 }; 498 static const PetscReal dd[1][2] = { 499 {-1., 1.} 500 }; 501 *nq = 1; 502 refweight = ww; 503 refinterp = ii; 504 refderiv = dd; 505 } break; 506 case QUADRATURE_GAUSS2: { 507 static const PetscReal ii[2][2] = { 508 {0.78867513459481287, 0.21132486540518713}, 509 {0.21132486540518713, 0.78867513459481287} 510 }; 511 static const PetscReal dd[2][2] = { 512 {-1., 1.}, 513 {-1., 1.} 514 }; 515 static const PetscReal ww[2] = {0.5, 0.5}; 516 *nq = 2; 517 refweight = ww; 518 refinterp = ii; 519 refderiv = dd; 520 } break; 521 case QUADRATURE_GAUSS3: { 522 static const PetscReal ii[3][2] = { 523 {0.8872983346207417, 0.1127016653792583}, 524 {0.5, 0.5 }, 525 {0.1127016653792583, 0.8872983346207417} 526 }; 527 static const PetscReal dd[3][2] = { 528 {-1, 1}, 529 {-1, 1}, 530 {-1, 1} 531 }; 532 static const PetscReal ww[3] = {5. / 18, 8. / 18, 5. / 18}; 533 *nq = 3; 534 refweight = ww; 535 refinterp = ii; 536 refderiv = dd; 537 } break; 538 case QUADRATURE_GAUSS4: { 539 static const PetscReal ii[][2] = { 540 {0.93056815579702623, 0.069431844202973658}, 541 {0.66999052179242813, 0.33000947820757187 }, 542 {0.33000947820757187, 0.66999052179242813 }, 543 {0.069431844202973658, 0.93056815579702623 } 544 }; 545 static const PetscReal dd[][2] = { 546 {-1, 1}, 547 {-1, 1}, 548 {-1, 1}, 549 {-1, 1} 550 }; 551 static const PetscReal ww[] = {0.17392742256872692, 0.3260725774312731, 0.3260725774312731, 0.17392742256872692}; 552 553 *nq = 4; 554 refweight = ww; 555 refinterp = ii; 556 refderiv = dd; 557 } break; 558 case QUADRATURE_LOBATTO2: { 559 static const PetscReal ii[2][2] = { 560 {1., 0.}, 561 {0., 1.} 562 }; 563 static const PetscReal dd[2][2] = { 564 {-1., 1.}, 565 {-1., 1.} 566 }; 567 static const PetscReal ww[2] = {0.5, 0.5}; 568 *nq = 2; 569 refweight = ww; 570 refinterp = ii; 571 refderiv = dd; 572 } break; 573 case QUADRATURE_LOBATTO3: { 574 static const PetscReal ii[3][2] = { 575 {1, 0 }, 576 {0.5, 0.5}, 577 {0, 1 } 578 }; 579 static const PetscReal dd[3][2] = { 580 {-1, 1}, 581 {-1, 1}, 582 {-1, 1} 583 }; 584 static const PetscReal ww[3] = {1. / 6, 4. / 6, 1. / 6}; 585 *nq = 3; 586 refweight = ww; 587 refinterp = ii; 588 refderiv = dd; 589 } break; 590 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unknown quadrature %d", (int)rd->quadrature); 591 } 592 593 for (q = 0; q < *nq; q++) { 594 weight[q] = refweight[q] * hx; 595 for (j = 0; j < 2; j++) { 596 interp[q][j] = refinterp[q][j]; 597 deriv[q][j] = refderiv[q][j] / hx; 598 } 599 } 600 PetscFunctionReturn(0); 601 } 602 603 /* 604 Finite element version 605 */ 606 static PetscErrorCode RDIFunction_FE(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ctx) { 607 RD rd = (RD)ctx; 608 RDNode *x, *x0, *xdot, *f; 609 Vec X0loc, Xloc, Xloc_t, Floc; 610 PetscReal hx, Theta, dt, weight[5], interp[5][2], deriv[5][2]; 611 DMDALocalInfo info; 612 PetscInt i, j, q, nq; 613 614 PetscFunctionBeginUser; 615 PetscCall(RDGetLocalArrays(rd, ts, X, Xdot, &Theta, &dt, &X0loc, &x0, &Xloc, &x, &Xloc_t, &xdot)); 616 617 PetscCall(DMGetLocalVector(rd->da, &Floc)); 618 PetscCall(VecZeroEntries(Floc)); 619 PetscCall(DMDAVecGetArray(rd->da, Floc, &f)); 620 PetscCall(DMDAGetLocalInfo(rd->da, &info)); 621 622 /* Set up shape functions and quadrature for elements (assumes a uniform grid) */ 623 hx = rd->L / (info.mx - 1); 624 PetscCall(RDGetQuadrature(rd, hx, &nq, weight, interp, deriv)); 625 626 for (i = info.xs; i < PetscMin(info.xs + info.xm, info.mx - 1); i++) { 627 for (q = 0; q < nq; q++) { 628 PetscReal rho = rd->rho; 629 PetscScalar Em_t, rad, D_R, D0_R; 630 RDNode n, n0, nx, n0x, nt, ntx; 631 RDEvaluate(interp, deriv, q, x, i, &n, &nx); 632 RDEvaluate(interp, deriv, q, x0, i, &n0, &n0x); 633 RDEvaluate(interp, deriv, q, xdot, i, &nt, &ntx); 634 635 rad = (1. - Theta) * RDRadiation(rd, &n0, 0) + Theta * RDRadiation(rd, &n, 0); 636 if (rd->endpoint) { 637 PetscScalar Em0, Em1; 638 RDMaterialEnergy(rd, &n0, &Em0, NULL); 639 RDMaterialEnergy(rd, &n, &Em1, NULL); 640 Em_t = (Em1 - Em0) / dt; 641 } else { 642 RDNode dEm; 643 RDMaterialEnergy(rd, &n, NULL, &dEm); 644 Em_t = dEm.E * nt.E + dEm.T * nt.T; 645 } 646 RDDiffusionCoefficient(rd, PETSC_TRUE, &n0, &n0x, &D0_R, 0, 0); 647 RDDiffusionCoefficient(rd, PETSC_TRUE, &n, &nx, &D_R, 0, 0); 648 for (j = 0; j < 2; j++) { 649 f[i + j].E += (deriv[q][j] * weight[q] * ((1. - Theta) * D0_R * n0x.E + Theta * D_R * nx.E) + interp[q][j] * weight[q] * (nt.E - rad)); 650 f[i + j].T += interp[q][j] * weight[q] * (rho * Em_t + rad); 651 } 652 } 653 } 654 if (info.xs == 0) { 655 switch (rd->leftbc) { 656 case BC_ROBIN: { 657 PetscScalar D_R, D_R_bc; 658 PetscReal ratio, bcTheta = rd->bcmidpoint ? Theta : 1.; 659 RDNode n, nx; 660 661 n.E = (1 - bcTheta) * x0[0].E + bcTheta * x[0].E; 662 n.T = (1 - bcTheta) * x0[0].T + bcTheta * x[0].T; 663 nx.E = (x[1].E - x[0].E) / hx; 664 nx.T = (x[1].T - x[0].T) / hx; 665 RDDiffusionCoefficient(rd, PETSC_TRUE, &n, &nx, &D_R, 0, 0); 666 RDDiffusionCoefficient(rd, rd->bclimit, &n, &nx, &D_R_bc, 0, 0); 667 ratio = PetscRealPart(D_R / D_R_bc); 668 PetscCheck(ratio <= 1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Limited diffusivity is greater than unlimited"); 669 PetscCheck(ratio >= 1e-3, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Heavily limited diffusivity"); 670 f[0].E += -ratio * 0.5 * (rd->Eapplied - n.E); 671 } break; 672 case BC_NEUMANN: 673 /* homogeneous Neumann is the natural condition */ 674 break; 675 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Case %" PetscInt_FMT, rd->initial); 676 } 677 } 678 679 PetscCall(RDRestoreLocalArrays(rd, &X0loc, &x0, &Xloc, &x, &Xloc_t, &xdot)); 680 PetscCall(DMDAVecRestoreArray(rd->da, Floc, &f)); 681 PetscCall(VecZeroEntries(F)); 682 PetscCall(DMLocalToGlobalBegin(rd->da, Floc, ADD_VALUES, F)); 683 PetscCall(DMLocalToGlobalEnd(rd->da, Floc, ADD_VALUES, F)); 684 PetscCall(DMRestoreLocalVector(rd->da, &Floc)); 685 686 if (rd->monitor_residual) PetscCall(RDStateView(rd, X, Xdot, F)); 687 PetscFunctionReturn(0); 688 } 689 690 static PetscErrorCode RDIJacobian_FE(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat A, Mat B, void *ctx) { 691 RD rd = (RD)ctx; 692 RDNode *x, *x0, *xdot; 693 Vec X0loc, Xloc, Xloc_t; 694 PetscReal hx, Theta, dt, weight[5], interp[5][2], deriv[5][2]; 695 DMDALocalInfo info; 696 PetscInt i, j, k, q, nq; 697 PetscScalar K[4][4]; 698 699 PetscFunctionBeginUser; 700 PetscCall(RDGetLocalArrays(rd, ts, X, Xdot, &Theta, &dt, &X0loc, &x0, &Xloc, &x, &Xloc_t, &xdot)); 701 PetscCall(DMDAGetLocalInfo(rd->da, &info)); 702 hx = rd->L / (info.mx - 1); 703 PetscCall(RDGetQuadrature(rd, hx, &nq, weight, interp, deriv)); 704 PetscCall(MatZeroEntries(B)); 705 for (i = info.xs; i < PetscMin(info.xs + info.xm, info.mx - 1); i++) { 706 PetscInt rc[2]; 707 708 rc[0] = i; 709 rc[1] = i + 1; 710 PetscCall(PetscMemzero(K, sizeof(K))); 711 for (q = 0; q < nq; q++) { 712 PetscScalar D_R; 713 PETSC_UNUSED PetscScalar rad; 714 RDNode n, nx, nt, ntx, drad, dD_R, dxD_R, dEm; 715 RDEvaluate(interp, deriv, q, x, i, &n, &nx); 716 RDEvaluate(interp, deriv, q, xdot, i, &nt, &ntx); 717 rad = RDRadiation(rd, &n, &drad); 718 RDDiffusionCoefficient(rd, PETSC_TRUE, &n, &nx, &D_R, &dD_R, &dxD_R); 719 RDMaterialEnergy(rd, &n, NULL, &dEm); 720 for (j = 0; j < 2; j++) { 721 for (k = 0; k < 2; k++) { 722 K[j * 2 + 0][k * 2 + 0] += (+interp[q][j] * weight[q] * (a - drad.E) * interp[q][k] + deriv[q][j] * weight[q] * ((D_R + dxD_R.E * nx.E) * deriv[q][k] + dD_R.E * nx.E * interp[q][k])); 723 K[j * 2 + 0][k * 2 + 1] += (+interp[q][j] * weight[q] * (-drad.T * interp[q][k]) + deriv[q][j] * weight[q] * (dxD_R.T * deriv[q][k] + dD_R.T * interp[q][k]) * nx.E); 724 K[j * 2 + 1][k * 2 + 0] += interp[q][j] * weight[q] * drad.E * interp[q][k]; 725 K[j * 2 + 1][k * 2 + 1] += interp[q][j] * weight[q] * (a * rd->rho * dEm.T + drad.T) * interp[q][k]; 726 } 727 } 728 } 729 PetscCall(MatSetValuesBlocked(B, 2, rc, 2, rc, &K[0][0], ADD_VALUES)); 730 } 731 if (info.xs == 0) { 732 switch (rd->leftbc) { 733 case BC_ROBIN: { 734 PetscScalar D_R, D_R_bc; 735 PetscReal ratio; 736 RDNode n, nx; 737 738 n.E = (1 - Theta) * x0[0].E + Theta * x[0].E; 739 n.T = (1 - Theta) * x0[0].T + Theta * x[0].T; 740 nx.E = (x[1].E - x[0].E) / hx; 741 nx.T = (x[1].T - x[0].T) / hx; 742 RDDiffusionCoefficient(rd, PETSC_TRUE, &n, &nx, &D_R, 0, 0); 743 RDDiffusionCoefficient(rd, rd->bclimit, &n, &nx, &D_R_bc, 0, 0); 744 ratio = PetscRealPart(D_R / D_R_bc); 745 PetscCall(MatSetValue(B, 0, 0, ratio * 0.5, ADD_VALUES)); 746 } break; 747 case BC_NEUMANN: 748 /* homogeneous Neumann is the natural condition */ 749 break; 750 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Case %" PetscInt_FMT, rd->initial); 751 } 752 } 753 754 PetscCall(RDRestoreLocalArrays(rd, &X0loc, &x0, &Xloc, &x, &Xloc_t, &xdot)); 755 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 756 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 757 if (A != B) { 758 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 759 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 760 } 761 PetscFunctionReturn(0); 762 } 763 764 /* Temperature that is in equilibrium with the radiation density */ 765 static PetscScalar RDRadiationTemperature(RD rd, PetscScalar E) { 766 return PetscPowScalar(E * rd->c / (4. * rd->sigma_b), 0.25); 767 } 768 769 static PetscErrorCode RDInitialState(RD rd, Vec X) { 770 DMDALocalInfo info; 771 PetscInt i; 772 RDNode *x; 773 774 PetscFunctionBeginUser; 775 PetscCall(DMDAGetLocalInfo(rd->da, &info)); 776 PetscCall(DMDAVecGetArray(rd->da, X, &x)); 777 for (i = info.xs; i < info.xs + info.xm; i++) { 778 PetscReal coord = i * rd->L / (info.mx - 1); 779 switch (rd->initial) { 780 case 1: 781 x[i].E = 0.001; 782 x[i].T = RDRadiationTemperature(rd, x[i].E); 783 break; 784 case 2: 785 x[i].E = 0.001 + 100. * PetscExpReal(-PetscSqr(coord / 0.1)); 786 x[i].T = RDRadiationTemperature(rd, x[i].E); 787 break; 788 case 3: 789 x[i].E = 7.56e-2 * rd->unit.Joule / PetscPowScalarInt(rd->unit.meter, 3); 790 x[i].T = RDRadiationTemperature(rd, x[i].E); 791 break; 792 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No initial state %" PetscInt_FMT, rd->initial); 793 } 794 } 795 PetscCall(DMDAVecRestoreArray(rd->da, X, &x)); 796 PetscFunctionReturn(0); 797 } 798 799 static PetscErrorCode RDView(RD rd, Vec X, PetscViewer viewer) { 800 Vec Y; 801 const RDNode *x; 802 PetscScalar *y; 803 PetscInt i, m, M; 804 const PetscInt *lx; 805 DM da; 806 MPI_Comm comm; 807 808 PetscFunctionBeginUser; 809 /* 810 Create a DMDA (one dof per node, zero stencil width, same layout) to hold Trad 811 (radiation temperature). It is not necessary to create a DMDA for this, but this way 812 output and visualization will have meaningful variable names and correct scales. 813 */ 814 PetscCall(DMDAGetInfo(rd->da, 0, &M, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)); 815 PetscCall(DMDAGetOwnershipRanges(rd->da, &lx, 0, 0)); 816 PetscCall(PetscObjectGetComm((PetscObject)rd->da, &comm)); 817 PetscCall(DMDACreate1d(comm, DM_BOUNDARY_NONE, M, 1, 0, lx, &da)); 818 PetscCall(DMSetFromOptions(da)); 819 PetscCall(DMSetUp(da)); 820 PetscCall(DMDASetUniformCoordinates(da, 0., rd->L, 0., 0., 0., 0.)); 821 PetscCall(DMDASetFieldName(da, 0, "T_rad")); 822 PetscCall(DMCreateGlobalVector(da, &Y)); 823 824 /* Compute the radiation temperature from the solution at each node */ 825 PetscCall(VecGetLocalSize(Y, &m)); 826 PetscCall(VecGetArrayRead(X, (const PetscScalar **)&x)); 827 PetscCall(VecGetArray(Y, &y)); 828 for (i = 0; i < m; i++) y[i] = RDRadiationTemperature(rd, x[i].E); 829 PetscCall(VecRestoreArrayRead(X, (const PetscScalar **)&x)); 830 PetscCall(VecRestoreArray(Y, &y)); 831 832 PetscCall(VecView(Y, viewer)); 833 PetscCall(VecDestroy(&Y)); 834 PetscCall(DMDestroy(&da)); 835 PetscFunctionReturn(0); 836 } 837 838 static PetscErrorCode RDTestDifferentiation(RD rd) { 839 MPI_Comm comm; 840 RDNode n, nx; 841 PetscScalar epsilon; 842 843 PetscFunctionBeginUser; 844 PetscCall(PetscObjectGetComm((PetscObject)rd->da, &comm)); 845 epsilon = 1e-8; 846 { 847 RDNode dEm, fdEm; 848 PetscScalar T0 = 1000., T1 = T0 * (1. + epsilon), Em0, Em1; 849 n.E = 1.; 850 n.T = T0; 851 rd->MaterialEnergy(rd, &n, &Em0, &dEm); 852 n.E = 1. + epsilon; 853 n.T = T0; 854 rd->MaterialEnergy(rd, &n, &Em1, 0); 855 fdEm.E = (Em1 - Em0) / epsilon; 856 n.E = 1.; 857 n.T = T1; 858 rd->MaterialEnergy(rd, &n, &Em1, 0); 859 fdEm.T = (Em1 - Em0) / (T0 * epsilon); 860 PetscCall(PetscPrintf(comm, "dEm {%g,%g}, fdEm {%g,%g}, diff {%g,%g}\n", (double)PetscRealPart(dEm.E), (double)PetscRealPart(dEm.T), (double)PetscRealPart(fdEm.E), (double)PetscRealPart(fdEm.T), (double)PetscRealPart(dEm.E - fdEm.E), 861 (double)PetscRealPart(dEm.T - fdEm.T))); 862 } 863 { 864 PetscScalar D0, D; 865 RDNode dD, dxD, fdD, fdxD; 866 n.E = 1.; 867 n.T = 1.; 868 nx.E = 1.; 869 n.T = 1.; 870 RDDiffusionCoefficient(rd, rd->bclimit, &n, &nx, &D0, &dD, &dxD); 871 n.E = 1. + epsilon; 872 n.T = 1.; 873 nx.E = 1.; 874 n.T = 1.; 875 RDDiffusionCoefficient(rd, rd->bclimit, &n, &nx, &D, 0, 0); 876 fdD.E = (D - D0) / epsilon; 877 n.E = 1; 878 n.T = 1. + epsilon; 879 nx.E = 1.; 880 n.T = 1.; 881 RDDiffusionCoefficient(rd, rd->bclimit, &n, &nx, &D, 0, 0); 882 fdD.T = (D - D0) / epsilon; 883 n.E = 1; 884 n.T = 1.; 885 nx.E = 1. + epsilon; 886 n.T = 1.; 887 RDDiffusionCoefficient(rd, rd->bclimit, &n, &nx, &D, 0, 0); 888 fdxD.E = (D - D0) / epsilon; 889 n.E = 1; 890 n.T = 1.; 891 nx.E = 1.; 892 n.T = 1. + epsilon; 893 RDDiffusionCoefficient(rd, rd->bclimit, &n, &nx, &D, 0, 0); 894 fdxD.T = (D - D0) / epsilon; 895 PetscCall(PetscPrintf(comm, "dD {%g,%g}, fdD {%g,%g}, diff {%g,%g}\n", (double)PetscRealPart(dD.E), (double)PetscRealPart(dD.T), (double)PetscRealPart(fdD.E), (double)PetscRealPart(fdD.T), (double)PetscRealPart(dD.E - fdD.E), 896 (double)PetscRealPart(dD.T - fdD.T))); 897 PetscCall(PetscPrintf(comm, "dxD {%g,%g}, fdxD {%g,%g}, diffx {%g,%g}\n", (double)PetscRealPart(dxD.E), (double)PetscRealPart(dxD.T), (double)PetscRealPart(fdxD.E), (double)PetscRealPart(fdxD.T), (double)PetscRealPart(dxD.E - fdxD.E), 898 (double)PetscRealPart(dxD.T - fdxD.T))); 899 } 900 { 901 PetscInt i; 902 PetscReal hx = 1.; 903 PetscScalar a0; 904 RDNode n0[3], n1[3], d[3], fd[3]; 905 906 n0[0].E = 1.; 907 n0[0].T = 1.; 908 n0[1].E = 5.; 909 n0[1].T = 3.; 910 n0[2].E = 4.; 911 n0[2].T = 2.; 912 a0 = RDDiffusion(rd, hx, n0, 1, d); 913 for (i = 0; i < 3; i++) { 914 PetscCall(PetscMemcpy(n1, n0, sizeof(n0))); 915 n1[i].E += epsilon; 916 fd[i].E = (RDDiffusion(rd, hx, n1, 1, 0) - a0) / epsilon; 917 PetscCall(PetscMemcpy(n1, n0, sizeof(n0))); 918 n1[i].T += epsilon; 919 fd[i].T = (RDDiffusion(rd, hx, n1, 1, 0) - a0) / epsilon; 920 PetscCall(PetscPrintf(comm, "ddiff[%" PetscInt_FMT "] {%g,%g}, fd {%g %g}, diff {%g,%g}\n", i, (double)PetscRealPart(d[i].E), (double)PetscRealPart(d[i].T), (double)PetscRealPart(fd[i].E), (double)PetscRealPart(fd[i].T), 921 (double)PetscRealPart(d[i].E - fd[i].E), (double)PetscRealPart(d[i].T - fd[i].T))); 922 } 923 } 924 { 925 PetscScalar rad0, rad; 926 RDNode drad, fdrad; 927 n.E = 1.; 928 n.T = 1.; 929 rad0 = RDRadiation(rd, &n, &drad); 930 n.E = 1. + epsilon; 931 n.T = 1.; 932 rad = RDRadiation(rd, &n, 0); 933 fdrad.E = (rad - rad0) / epsilon; 934 n.E = 1.; 935 n.T = 1. + epsilon; 936 rad = RDRadiation(rd, &n, 0); 937 fdrad.T = (rad - rad0) / epsilon; 938 PetscCall(PetscPrintf(comm, "drad {%g,%g}, fdrad {%g,%g}, diff {%g,%g}\n", (double)PetscRealPart(drad.E), (double)PetscRealPart(drad.T), (double)PetscRealPart(fdrad.E), (double)PetscRealPart(fdrad.T), (double)PetscRealPart(drad.E - drad.E), 939 (double)PetscRealPart(drad.T - fdrad.T))); 940 } 941 PetscFunctionReturn(0); 942 } 943 944 static PetscErrorCode RDCreate(MPI_Comm comm, RD *inrd) { 945 RD rd; 946 PetscReal meter = 0, kilogram = 0, second = 0, Kelvin = 0, Joule = 0, Watt = 0; 947 948 PetscFunctionBeginUser; 949 *inrd = 0; 950 PetscCall(PetscNew(&rd)); 951 952 PetscOptionsBegin(comm, NULL, "Options for nonequilibrium radiation-diffusion with RD ionization", NULL); 953 { 954 rd->initial = 1; 955 PetscCall(PetscOptionsInt("-rd_initial", "Initial condition (1=Marshak, 2=Blast, 3=Marshak+)", "", rd->initial, &rd->initial, 0)); 956 switch (rd->initial) { 957 case 1: 958 case 2: 959 rd->unit.kilogram = 1.; 960 rd->unit.meter = 1.; 961 rd->unit.second = 1.; 962 rd->unit.Kelvin = 1.; 963 break; 964 case 3: 965 rd->unit.kilogram = 1.e12; 966 rd->unit.meter = 1.; 967 rd->unit.second = 1.e9; 968 rd->unit.Kelvin = 1.; 969 break; 970 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unknown initial condition %" PetscInt_FMT, rd->initial); 971 } 972 /* Fundamental units */ 973 PetscCall(PetscOptionsReal("-rd_unit_meter", "Length of 1 meter in nondimensional units", "", rd->unit.meter, &rd->unit.meter, 0)); 974 PetscCall(PetscOptionsReal("-rd_unit_kilogram", "Mass of 1 kilogram in nondimensional units", "", rd->unit.kilogram, &rd->unit.kilogram, 0)); 975 PetscCall(PetscOptionsReal("-rd_unit_second", "Time of a second in nondimensional units", "", rd->unit.second, &rd->unit.second, 0)); 976 PetscCall(PetscOptionsReal("-rd_unit_Kelvin", "Temperature of a Kelvin in nondimensional units", "", rd->unit.Kelvin, &rd->unit.Kelvin, 0)); 977 /* Derived units */ 978 rd->unit.Joule = rd->unit.kilogram * PetscSqr(rd->unit.meter / rd->unit.second); 979 rd->unit.Watt = rd->unit.Joule / rd->unit.second; 980 /* Local aliases */ 981 meter = rd->unit.meter; 982 kilogram = rd->unit.kilogram; 983 second = rd->unit.second; 984 Kelvin = rd->unit.Kelvin; 985 Joule = rd->unit.Joule; 986 Watt = rd->unit.Watt; 987 988 PetscCall(PetscOptionsBool("-rd_monitor_residual", "Display residuals every time they are evaluated", "", rd->monitor_residual, &rd->monitor_residual, NULL)); 989 PetscCall(PetscOptionsEnum("-rd_discretization", "Discretization type", "", DiscretizationTypes, (PetscEnum)rd->discretization, (PetscEnum *)&rd->discretization, NULL)); 990 if (rd->discretization == DISCRETIZATION_FE) { 991 rd->quadrature = QUADRATURE_GAUSS2; 992 PetscCall(PetscOptionsEnum("-rd_quadrature", "Finite element quadrature", "", QuadratureTypes, (PetscEnum)rd->quadrature, (PetscEnum *)&rd->quadrature, NULL)); 993 } 994 PetscCall(PetscOptionsEnum("-rd_jacobian", "Type of finite difference Jacobian", "", JacobianTypes, (PetscEnum)rd->jacobian, (PetscEnum *)&rd->jacobian, NULL)); 995 switch (rd->initial) { 996 case 1: 997 rd->leftbc = BC_ROBIN; 998 rd->Eapplied = 4 * rd->unit.Joule / PetscPowRealInt(rd->unit.meter, 3); 999 rd->L = 1. * rd->unit.meter; 1000 rd->beta = 3.0; 1001 rd->gamma = 3.0; 1002 rd->final_time = 3 * second; 1003 break; 1004 case 2: 1005 rd->leftbc = BC_NEUMANN; 1006 rd->Eapplied = 0.; 1007 rd->L = 1. * rd->unit.meter; 1008 rd->beta = 3.0; 1009 rd->gamma = 3.0; 1010 rd->final_time = 1 * second; 1011 break; 1012 case 3: 1013 rd->leftbc = BC_ROBIN; 1014 rd->Eapplied = 7.503e6 * rd->unit.Joule / PetscPowRealInt(rd->unit.meter, 3); 1015 rd->L = 5. * rd->unit.meter; 1016 rd->beta = 3.5; 1017 rd->gamma = 3.5; 1018 rd->final_time = 20e-9 * second; 1019 break; 1020 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Initial %" PetscInt_FMT, rd->initial); 1021 } 1022 PetscCall(PetscOptionsEnum("-rd_leftbc", "Left boundary condition", "", BCTypes, (PetscEnum)rd->leftbc, (PetscEnum *)&rd->leftbc, NULL)); 1023 PetscCall(PetscOptionsReal("-rd_E_applied", "Radiation flux at left end of domain", "", rd->Eapplied, &rd->Eapplied, NULL)); 1024 PetscCall(PetscOptionsReal("-rd_beta", "Thermal exponent for photon absorption", "", rd->beta, &rd->beta, NULL)); 1025 PetscCall(PetscOptionsReal("-rd_gamma", "Thermal exponent for diffusion coefficient", "", rd->gamma, &rd->gamma, NULL)); 1026 PetscCall(PetscOptionsBool("-rd_view_draw", "Draw final solution", "", rd->view_draw, &rd->view_draw, NULL)); 1027 PetscCall(PetscOptionsBool("-rd_endpoint", "Discretize using endpoints (like trapezoid rule) instead of midpoint", "", rd->endpoint, &rd->endpoint, NULL)); 1028 PetscCall(PetscOptionsBool("-rd_bcmidpoint", "Impose the boundary condition at the midpoint (Theta) of the interval", "", rd->bcmidpoint, &rd->bcmidpoint, NULL)); 1029 PetscCall(PetscOptionsBool("-rd_bclimit", "Limit diffusion coefficient in definition of Robin boundary condition", "", rd->bclimit, &rd->bclimit, NULL)); 1030 PetscCall(PetscOptionsBool("-rd_test_diff", "Test differentiation in constitutive relations", "", rd->test_diff, &rd->test_diff, NULL)); 1031 PetscCall(PetscOptionsString("-rd_view_binary", "File name to hold final solution", "", rd->view_binary, rd->view_binary, sizeof(rd->view_binary), NULL)); 1032 } 1033 PetscOptionsEnd(); 1034 1035 switch (rd->initial) { 1036 case 1: 1037 case 2: 1038 rd->rho = 1.; 1039 rd->c = 1.; 1040 rd->K_R = 1.; 1041 rd->K_p = 1.; 1042 rd->sigma_b = 0.25; 1043 rd->MaterialEnergy = RDMaterialEnergy_Reduced; 1044 break; 1045 case 3: 1046 /* Table 2 */ 1047 rd->rho = 1.17e-3 * kilogram / (meter * meter * meter); /* density */ 1048 rd->K_R = 7.44e18 * PetscPowRealInt(meter, 5) * PetscPowReal(Kelvin, 3.5) * PetscPowRealInt(kilogram, -2); /* */ 1049 rd->K_p = 2.33e20 * PetscPowRealInt(meter, 5) * PetscPowReal(Kelvin, 3.5) * PetscPowRealInt(kilogram, -2); /* */ 1050 rd->I_H = 2.179e-18 * Joule; /* Hydrogen ionization potential */ 1051 rd->m_p = 1.673e-27 * kilogram; /* proton mass */ 1052 rd->m_e = 9.109e-31 * kilogram; /* electron mass */ 1053 rd->h = 6.626e-34 * Joule * second; /* Planck's constant */ 1054 rd->k = 1.381e-23 * Joule / Kelvin; /* Boltzman constant */ 1055 rd->c = 3.00e8 * meter / second; /* speed of light */ 1056 rd->sigma_b = 5.67e-8 * Watt * PetscPowRealInt(meter, -2) * PetscPowRealInt(Kelvin, -4); /* Stefan-Boltzman constant */ 1057 rd->MaterialEnergy = RDMaterialEnergy_Saha; 1058 break; 1059 } 1060 1061 PetscCall(DMDACreate1d(comm, DM_BOUNDARY_NONE, 20, sizeof(RDNode) / sizeof(PetscScalar), 1, NULL, &rd->da)); 1062 PetscCall(DMSetFromOptions(rd->da)); 1063 PetscCall(DMSetUp(rd->da)); 1064 PetscCall(DMDASetFieldName(rd->da, 0, "E")); 1065 PetscCall(DMDASetFieldName(rd->da, 1, "T")); 1066 PetscCall(DMDASetUniformCoordinates(rd->da, 0., 1., 0., 0., 0., 0.)); 1067 1068 *inrd = rd; 1069 PetscFunctionReturn(0); 1070 } 1071 1072 int main(int argc, char *argv[]) { 1073 RD rd; 1074 TS ts; 1075 SNES snes; 1076 Vec X; 1077 Mat A, B; 1078 PetscInt steps; 1079 PetscReal ftime; 1080 1081 PetscFunctionBeginUser; 1082 PetscCall(PetscInitialize(&argc, &argv, 0, NULL)); 1083 PetscCall(RDCreate(PETSC_COMM_WORLD, &rd)); 1084 PetscCall(DMCreateGlobalVector(rd->da, &X)); 1085 PetscCall(DMSetMatType(rd->da, MATAIJ)); 1086 PetscCall(DMCreateMatrix(rd->da, &B)); 1087 PetscCall(RDInitialState(rd, X)); 1088 1089 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 1090 PetscCall(TSSetProblemType(ts, TS_NONLINEAR)); 1091 PetscCall(TSSetType(ts, TSTHETA)); 1092 PetscCall(TSSetDM(ts, rd->da)); 1093 switch (rd->discretization) { 1094 case DISCRETIZATION_FD: 1095 PetscCall(TSSetIFunction(ts, NULL, RDIFunction_FD, rd)); 1096 if (rd->jacobian == JACOBIAN_ANALYTIC) PetscCall(TSSetIJacobian(ts, B, B, RDIJacobian_FD, rd)); 1097 break; 1098 case DISCRETIZATION_FE: 1099 PetscCall(TSSetIFunction(ts, NULL, RDIFunction_FE, rd)); 1100 if (rd->jacobian == JACOBIAN_ANALYTIC) PetscCall(TSSetIJacobian(ts, B, B, RDIJacobian_FE, rd)); 1101 break; 1102 } 1103 PetscCall(TSSetMaxTime(ts, rd->final_time)); 1104 PetscCall(TSSetTimeStep(ts, 1e-3)); 1105 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); 1106 PetscCall(TSSetFromOptions(ts)); 1107 1108 A = B; 1109 PetscCall(TSGetSNES(ts, &snes)); 1110 switch (rd->jacobian) { 1111 case JACOBIAN_ANALYTIC: break; 1112 case JACOBIAN_MATRIXFREE: break; 1113 case JACOBIAN_FD_COLORING: { 1114 PetscCall(SNESSetJacobian(snes, A, B, SNESComputeJacobianDefaultColor, 0)); 1115 } break; 1116 case JACOBIAN_FD_FULL: PetscCall(SNESSetJacobian(snes, A, B, SNESComputeJacobianDefault, ts)); break; 1117 } 1118 1119 if (rd->test_diff) PetscCall(RDTestDifferentiation(rd)); 1120 PetscCall(TSSolve(ts, X)); 1121 PetscCall(TSGetSolveTime(ts, &ftime)); 1122 PetscCall(TSGetStepNumber(ts, &steps)); 1123 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Steps %" PetscInt_FMT " final time %g\n", steps, (double)ftime)); 1124 if (rd->view_draw) PetscCall(RDView(rd, X, PETSC_VIEWER_DRAW_WORLD)); 1125 if (rd->view_binary[0]) { 1126 PetscViewer viewer; 1127 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, rd->view_binary, FILE_MODE_WRITE, &viewer)); 1128 PetscCall(RDView(rd, X, viewer)); 1129 PetscCall(PetscViewerDestroy(&viewer)); 1130 } 1131 PetscCall(VecDestroy(&X)); 1132 PetscCall(MatDestroy(&B)); 1133 PetscCall(RDDestroy(&rd)); 1134 PetscCall(TSDestroy(&ts)); 1135 PetscCall(PetscFinalize()); 1136 return 0; 1137 } 1138 /*TEST 1139 1140 test: 1141 args: -da_grid_x 20 -rd_initial 1 -rd_discretization fd -rd_jacobian fd_coloring -rd_endpoint -ts_max_time 1 -ts_dt 2e-1 -ts_theta_initial_guess_extrapolate 0 -ts_monitor -snes_monitor_short -ksp_monitor_short 1142 requires: !single 1143 1144 test: 1145 suffix: 2 1146 args: -da_grid_x 20 -rd_initial 1 -rd_discretization fe -rd_quadrature lobatto2 -rd_jacobian fd_coloring -rd_endpoint -ts_max_time 1 -ts_dt 2e-1 -ts_theta_initial_guess_extrapolate 0 -ts_monitor -snes_monitor_short -ksp_monitor_short 1147 requires: !single 1148 1149 test: 1150 suffix: 3 1151 nsize: 2 1152 args: -da_grid_x 20 -rd_initial 1 -rd_discretization fd -rd_jacobian analytic -rd_endpoint -ts_max_time 3 -ts_dt 1e-1 -ts_theta_initial_guess_extrapolate 0 -ts_monitor -snes_monitor_short -ksp_monitor_short 1153 requires: !single 1154 1155 TEST*/ 1156