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] = 495 { 496 1. 497 }, 498 ii[1][2] = {{0.5, 0.5}}, dd[1][2] = {{-1., 1.}}; 499 *nq = 1; 500 refweight = ww; 501 refinterp = ii; 502 refderiv = dd; 503 } break; 504 case QUADRATURE_GAUSS2: { 505 static const PetscReal ii[2][2] = 506 { 507 {0.78867513459481287, 0.21132486540518713}, 508 {0.21132486540518713, 0.78867513459481287} 509 }, 510 dd[2][2] = {{-1., 1.}, {-1., 1.}}, ww[2] = {0.5, 0.5}; 511 *nq = 2; 512 refweight = ww; 513 refinterp = ii; 514 refderiv = dd; 515 } break; 516 case QUADRATURE_GAUSS3: { 517 static const PetscReal ii[3][2] = 518 { 519 {0.8872983346207417, 0.1127016653792583}, 520 {0.5, 0.5 }, 521 {0.1127016653792583, 0.8872983346207417} 522 }, 523 dd[3][2] = {{-1, 1}, {-1, 1}, {-1, 1}}, ww[3] = {5. / 18, 8. / 18, 5. / 18}; 524 *nq = 3; 525 refweight = ww; 526 refinterp = ii; 527 refderiv = dd; 528 } break; 529 case QUADRATURE_GAUSS4: { 530 static const PetscReal ii[][2] = 531 { 532 {0.93056815579702623, 0.069431844202973658}, 533 {0.66999052179242813, 0.33000947820757187 }, 534 {0.33000947820757187, 0.66999052179242813 }, 535 {0.069431844202973658, 0.93056815579702623 } 536 }, 537 dd[][2] = {{-1, 1}, {-1, 1}, {-1, 1}, {-1, 1}}, ww[] = {0.17392742256872692, 0.3260725774312731, 0.3260725774312731, 0.17392742256872692}; 538 539 *nq = 4; 540 refweight = ww; 541 refinterp = ii; 542 refderiv = dd; 543 } break; 544 case QUADRATURE_LOBATTO2: { 545 static const PetscReal ii[2][2] = 546 { 547 {1., 0.}, 548 {0., 1.} 549 }, 550 dd[2][2] = {{-1., 1.}, {-1., 1.}}, ww[2] = {0.5, 0.5}; 551 *nq = 2; 552 refweight = ww; 553 refinterp = ii; 554 refderiv = dd; 555 } break; 556 case QUADRATURE_LOBATTO3: { 557 static const PetscReal ii[3][2] = 558 { 559 {1, 0 }, 560 {0.5, 0.5}, 561 {0, 1 } 562 }, 563 dd[3][2] = {{-1, 1}, {-1, 1}, {-1, 1}}, ww[3] = {1. / 6, 4. / 6, 1. / 6}; 564 *nq = 3; 565 refweight = ww; 566 refinterp = ii; 567 refderiv = dd; 568 } break; 569 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unknown quadrature %d", (int)rd->quadrature); 570 } 571 572 for (q = 0; q < *nq; q++) { 573 weight[q] = refweight[q] * hx; 574 for (j = 0; j < 2; j++) { 575 interp[q][j] = refinterp[q][j]; 576 deriv[q][j] = refderiv[q][j] / hx; 577 } 578 } 579 PetscFunctionReturn(0); 580 } 581 582 /* 583 Finite element version 584 */ 585 static PetscErrorCode RDIFunction_FE(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ctx) { 586 RD rd = (RD)ctx; 587 RDNode *x, *x0, *xdot, *f; 588 Vec X0loc, Xloc, Xloc_t, Floc; 589 PetscReal hx, Theta, dt, weight[5], interp[5][2], deriv[5][2]; 590 DMDALocalInfo info; 591 PetscInt i, j, q, nq; 592 593 PetscFunctionBeginUser; 594 PetscCall(RDGetLocalArrays(rd, ts, X, Xdot, &Theta, &dt, &X0loc, &x0, &Xloc, &x, &Xloc_t, &xdot)); 595 596 PetscCall(DMGetLocalVector(rd->da, &Floc)); 597 PetscCall(VecZeroEntries(Floc)); 598 PetscCall(DMDAVecGetArray(rd->da, Floc, &f)); 599 PetscCall(DMDAGetLocalInfo(rd->da, &info)); 600 601 /* Set up shape functions and quadrature for elements (assumes a uniform grid) */ 602 hx = rd->L / (info.mx - 1); 603 PetscCall(RDGetQuadrature(rd, hx, &nq, weight, interp, deriv)); 604 605 for (i = info.xs; i < PetscMin(info.xs + info.xm, info.mx - 1); i++) { 606 for (q = 0; q < nq; q++) { 607 PetscReal rho = rd->rho; 608 PetscScalar Em_t, rad, D_R, D0_R; 609 RDNode n, n0, nx, n0x, nt, ntx; 610 RDEvaluate(interp, deriv, q, x, i, &n, &nx); 611 RDEvaluate(interp, deriv, q, x0, i, &n0, &n0x); 612 RDEvaluate(interp, deriv, q, xdot, i, &nt, &ntx); 613 614 rad = (1. - Theta) * RDRadiation(rd, &n0, 0) + Theta * RDRadiation(rd, &n, 0); 615 if (rd->endpoint) { 616 PetscScalar Em0, Em1; 617 RDMaterialEnergy(rd, &n0, &Em0, NULL); 618 RDMaterialEnergy(rd, &n, &Em1, NULL); 619 Em_t = (Em1 - Em0) / dt; 620 } else { 621 RDNode dEm; 622 RDMaterialEnergy(rd, &n, NULL, &dEm); 623 Em_t = dEm.E * nt.E + dEm.T * nt.T; 624 } 625 RDDiffusionCoefficient(rd, PETSC_TRUE, &n0, &n0x, &D0_R, 0, 0); 626 RDDiffusionCoefficient(rd, PETSC_TRUE, &n, &nx, &D_R, 0, 0); 627 for (j = 0; j < 2; j++) { 628 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)); 629 f[i + j].T += interp[q][j] * weight[q] * (rho * Em_t + rad); 630 } 631 } 632 } 633 if (info.xs == 0) { 634 switch (rd->leftbc) { 635 case BC_ROBIN: { 636 PetscScalar D_R, D_R_bc; 637 PetscReal ratio, bcTheta = rd->bcmidpoint ? Theta : 1.; 638 RDNode n, nx; 639 640 n.E = (1 - bcTheta) * x0[0].E + bcTheta * x[0].E; 641 n.T = (1 - bcTheta) * x0[0].T + bcTheta * x[0].T; 642 nx.E = (x[1].E - x[0].E) / hx; 643 nx.T = (x[1].T - x[0].T) / hx; 644 RDDiffusionCoefficient(rd, PETSC_TRUE, &n, &nx, &D_R, 0, 0); 645 RDDiffusionCoefficient(rd, rd->bclimit, &n, &nx, &D_R_bc, 0, 0); 646 ratio = PetscRealPart(D_R / D_R_bc); 647 PetscCheck(ratio <= 1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Limited diffusivity is greater than unlimited"); 648 PetscCheck(ratio >= 1e-3, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Heavily limited diffusivity"); 649 f[0].E += -ratio * 0.5 * (rd->Eapplied - n.E); 650 } break; 651 case BC_NEUMANN: 652 /* homogeneous Neumann is the natural condition */ 653 break; 654 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Case %" PetscInt_FMT, rd->initial); 655 } 656 } 657 658 PetscCall(RDRestoreLocalArrays(rd, &X0loc, &x0, &Xloc, &x, &Xloc_t, &xdot)); 659 PetscCall(DMDAVecRestoreArray(rd->da, Floc, &f)); 660 PetscCall(VecZeroEntries(F)); 661 PetscCall(DMLocalToGlobalBegin(rd->da, Floc, ADD_VALUES, F)); 662 PetscCall(DMLocalToGlobalEnd(rd->da, Floc, ADD_VALUES, F)); 663 PetscCall(DMRestoreLocalVector(rd->da, &Floc)); 664 665 if (rd->monitor_residual) PetscCall(RDStateView(rd, X, Xdot, F)); 666 PetscFunctionReturn(0); 667 } 668 669 static PetscErrorCode RDIJacobian_FE(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat A, Mat B, void *ctx) { 670 RD rd = (RD)ctx; 671 RDNode *x, *x0, *xdot; 672 Vec X0loc, Xloc, Xloc_t; 673 PetscReal hx, Theta, dt, weight[5], interp[5][2], deriv[5][2]; 674 DMDALocalInfo info; 675 PetscInt i, j, k, q, nq; 676 PetscScalar K[4][4]; 677 678 PetscFunctionBeginUser; 679 PetscCall(RDGetLocalArrays(rd, ts, X, Xdot, &Theta, &dt, &X0loc, &x0, &Xloc, &x, &Xloc_t, &xdot)); 680 PetscCall(DMDAGetLocalInfo(rd->da, &info)); 681 hx = rd->L / (info.mx - 1); 682 PetscCall(RDGetQuadrature(rd, hx, &nq, weight, interp, deriv)); 683 PetscCall(MatZeroEntries(B)); 684 for (i = info.xs; i < PetscMin(info.xs + info.xm, info.mx - 1); i++) { 685 PetscInt rc[2]; 686 687 rc[0] = i; 688 rc[1] = i + 1; 689 PetscCall(PetscMemzero(K, sizeof(K))); 690 for (q = 0; q < nq; q++) { 691 PetscScalar D_R; 692 PETSC_UNUSED PetscScalar rad; 693 RDNode n, nx, nt, ntx, drad, dD_R, dxD_R, dEm; 694 RDEvaluate(interp, deriv, q, x, i, &n, &nx); 695 RDEvaluate(interp, deriv, q, xdot, i, &nt, &ntx); 696 rad = RDRadiation(rd, &n, &drad); 697 RDDiffusionCoefficient(rd, PETSC_TRUE, &n, &nx, &D_R, &dD_R, &dxD_R); 698 RDMaterialEnergy(rd, &n, NULL, &dEm); 699 for (j = 0; j < 2; j++) { 700 for (k = 0; k < 2; k++) { 701 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])); 702 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); 703 K[j * 2 + 1][k * 2 + 0] += interp[q][j] * weight[q] * drad.E * interp[q][k]; 704 K[j * 2 + 1][k * 2 + 1] += interp[q][j] * weight[q] * (a * rd->rho * dEm.T + drad.T) * interp[q][k]; 705 } 706 } 707 } 708 PetscCall(MatSetValuesBlocked(B, 2, rc, 2, rc, &K[0][0], ADD_VALUES)); 709 } 710 if (info.xs == 0) { 711 switch (rd->leftbc) { 712 case BC_ROBIN: { 713 PetscScalar D_R, D_R_bc; 714 PetscReal ratio; 715 RDNode n, nx; 716 717 n.E = (1 - Theta) * x0[0].E + Theta * x[0].E; 718 n.T = (1 - Theta) * x0[0].T + Theta * x[0].T; 719 nx.E = (x[1].E - x[0].E) / hx; 720 nx.T = (x[1].T - x[0].T) / hx; 721 RDDiffusionCoefficient(rd, PETSC_TRUE, &n, &nx, &D_R, 0, 0); 722 RDDiffusionCoefficient(rd, rd->bclimit, &n, &nx, &D_R_bc, 0, 0); 723 ratio = PetscRealPart(D_R / D_R_bc); 724 PetscCall(MatSetValue(B, 0, 0, ratio * 0.5, ADD_VALUES)); 725 } break; 726 case BC_NEUMANN: 727 /* homogeneous Neumann is the natural condition */ 728 break; 729 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Case %" PetscInt_FMT, rd->initial); 730 } 731 } 732 733 PetscCall(RDRestoreLocalArrays(rd, &X0loc, &x0, &Xloc, &x, &Xloc_t, &xdot)); 734 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 735 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 736 if (A != B) { 737 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 738 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 739 } 740 PetscFunctionReturn(0); 741 } 742 743 /* Temperature that is in equilibrium with the radiation density */ 744 static PetscScalar RDRadiationTemperature(RD rd, PetscScalar E) { 745 return PetscPowScalar(E * rd->c / (4. * rd->sigma_b), 0.25); 746 } 747 748 static PetscErrorCode RDInitialState(RD rd, Vec X) { 749 DMDALocalInfo info; 750 PetscInt i; 751 RDNode *x; 752 753 PetscFunctionBeginUser; 754 PetscCall(DMDAGetLocalInfo(rd->da, &info)); 755 PetscCall(DMDAVecGetArray(rd->da, X, &x)); 756 for (i = info.xs; i < info.xs + info.xm; i++) { 757 PetscReal coord = i * rd->L / (info.mx - 1); 758 switch (rd->initial) { 759 case 1: 760 x[i].E = 0.001; 761 x[i].T = RDRadiationTemperature(rd, x[i].E); 762 break; 763 case 2: 764 x[i].E = 0.001 + 100. * PetscExpReal(-PetscSqr(coord / 0.1)); 765 x[i].T = RDRadiationTemperature(rd, x[i].E); 766 break; 767 case 3: 768 x[i].E = 7.56e-2 * rd->unit.Joule / PetscPowScalarInt(rd->unit.meter, 3); 769 x[i].T = RDRadiationTemperature(rd, x[i].E); 770 break; 771 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No initial state %" PetscInt_FMT, rd->initial); 772 } 773 } 774 PetscCall(DMDAVecRestoreArray(rd->da, X, &x)); 775 PetscFunctionReturn(0); 776 } 777 778 static PetscErrorCode RDView(RD rd, Vec X, PetscViewer viewer) { 779 Vec Y; 780 const RDNode *x; 781 PetscScalar *y; 782 PetscInt i, m, M; 783 const PetscInt *lx; 784 DM da; 785 MPI_Comm comm; 786 787 PetscFunctionBeginUser; 788 /* 789 Create a DMDA (one dof per node, zero stencil width, same layout) to hold Trad 790 (radiation temperature). It is not necessary to create a DMDA for this, but this way 791 output and visualization will have meaningful variable names and correct scales. 792 */ 793 PetscCall(DMDAGetInfo(rd->da, 0, &M, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)); 794 PetscCall(DMDAGetOwnershipRanges(rd->da, &lx, 0, 0)); 795 PetscCall(PetscObjectGetComm((PetscObject)rd->da, &comm)); 796 PetscCall(DMDACreate1d(comm, DM_BOUNDARY_NONE, M, 1, 0, lx, &da)); 797 PetscCall(DMSetFromOptions(da)); 798 PetscCall(DMSetUp(da)); 799 PetscCall(DMDASetUniformCoordinates(da, 0., rd->L, 0., 0., 0., 0.)); 800 PetscCall(DMDASetFieldName(da, 0, "T_rad")); 801 PetscCall(DMCreateGlobalVector(da, &Y)); 802 803 /* Compute the radiation temperature from the solution at each node */ 804 PetscCall(VecGetLocalSize(Y, &m)); 805 PetscCall(VecGetArrayRead(X, (const PetscScalar **)&x)); 806 PetscCall(VecGetArray(Y, &y)); 807 for (i = 0; i < m; i++) y[i] = RDRadiationTemperature(rd, x[i].E); 808 PetscCall(VecRestoreArrayRead(X, (const PetscScalar **)&x)); 809 PetscCall(VecRestoreArray(Y, &y)); 810 811 PetscCall(VecView(Y, viewer)); 812 PetscCall(VecDestroy(&Y)); 813 PetscCall(DMDestroy(&da)); 814 PetscFunctionReturn(0); 815 } 816 817 static PetscErrorCode RDTestDifferentiation(RD rd) { 818 MPI_Comm comm; 819 RDNode n, nx; 820 PetscScalar epsilon; 821 822 PetscFunctionBeginUser; 823 PetscCall(PetscObjectGetComm((PetscObject)rd->da, &comm)); 824 epsilon = 1e-8; 825 { 826 RDNode dEm, fdEm; 827 PetscScalar T0 = 1000., T1 = T0 * (1. + epsilon), Em0, Em1; 828 n.E = 1.; 829 n.T = T0; 830 rd->MaterialEnergy(rd, &n, &Em0, &dEm); 831 n.E = 1. + epsilon; 832 n.T = T0; 833 rd->MaterialEnergy(rd, &n, &Em1, 0); 834 fdEm.E = (Em1 - Em0) / epsilon; 835 n.E = 1.; 836 n.T = T1; 837 rd->MaterialEnergy(rd, &n, &Em1, 0); 838 fdEm.T = (Em1 - Em0) / (T0 * epsilon); 839 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), 840 (double)PetscRealPart(dEm.T - fdEm.T))); 841 } 842 { 843 PetscScalar D0, D; 844 RDNode dD, dxD, fdD, fdxD; 845 n.E = 1.; 846 n.T = 1.; 847 nx.E = 1.; 848 n.T = 1.; 849 RDDiffusionCoefficient(rd, rd->bclimit, &n, &nx, &D0, &dD, &dxD); 850 n.E = 1. + epsilon; 851 n.T = 1.; 852 nx.E = 1.; 853 n.T = 1.; 854 RDDiffusionCoefficient(rd, rd->bclimit, &n, &nx, &D, 0, 0); 855 fdD.E = (D - D0) / epsilon; 856 n.E = 1; 857 n.T = 1. + epsilon; 858 nx.E = 1.; 859 n.T = 1.; 860 RDDiffusionCoefficient(rd, rd->bclimit, &n, &nx, &D, 0, 0); 861 fdD.T = (D - D0) / epsilon; 862 n.E = 1; 863 n.T = 1.; 864 nx.E = 1. + epsilon; 865 n.T = 1.; 866 RDDiffusionCoefficient(rd, rd->bclimit, &n, &nx, &D, 0, 0); 867 fdxD.E = (D - D0) / epsilon; 868 n.E = 1; 869 n.T = 1.; 870 nx.E = 1.; 871 n.T = 1. + epsilon; 872 RDDiffusionCoefficient(rd, rd->bclimit, &n, &nx, &D, 0, 0); 873 fdxD.T = (D - D0) / epsilon; 874 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), 875 (double)PetscRealPart(dD.T - fdD.T))); 876 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), 877 (double)PetscRealPart(dxD.T - fdxD.T))); 878 } 879 { 880 PetscInt i; 881 PetscReal hx = 1.; 882 PetscScalar a0; 883 RDNode n0[3], n1[3], d[3], fd[3]; 884 885 n0[0].E = 1.; 886 n0[0].T = 1.; 887 n0[1].E = 5.; 888 n0[1].T = 3.; 889 n0[2].E = 4.; 890 n0[2].T = 2.; 891 a0 = RDDiffusion(rd, hx, n0, 1, d); 892 for (i = 0; i < 3; i++) { 893 PetscCall(PetscMemcpy(n1, n0, sizeof(n0))); 894 n1[i].E += epsilon; 895 fd[i].E = (RDDiffusion(rd, hx, n1, 1, 0) - a0) / epsilon; 896 PetscCall(PetscMemcpy(n1, n0, sizeof(n0))); 897 n1[i].T += epsilon; 898 fd[i].T = (RDDiffusion(rd, hx, n1, 1, 0) - a0) / epsilon; 899 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), 900 (double)PetscRealPart(d[i].E - fd[i].E), (double)PetscRealPart(d[i].T - fd[i].T))); 901 } 902 } 903 { 904 PetscScalar rad0, rad; 905 RDNode drad, fdrad; 906 n.E = 1.; 907 n.T = 1.; 908 rad0 = RDRadiation(rd, &n, &drad); 909 n.E = 1. + epsilon; 910 n.T = 1.; 911 rad = RDRadiation(rd, &n, 0); 912 fdrad.E = (rad - rad0) / epsilon; 913 n.E = 1.; 914 n.T = 1. + epsilon; 915 rad = RDRadiation(rd, &n, 0); 916 fdrad.T = (rad - rad0) / epsilon; 917 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), 918 (double)PetscRealPart(drad.T - fdrad.T))); 919 } 920 PetscFunctionReturn(0); 921 } 922 923 static PetscErrorCode RDCreate(MPI_Comm comm, RD *inrd) { 924 RD rd; 925 PetscReal meter = 0, kilogram = 0, second = 0, Kelvin = 0, Joule = 0, Watt = 0; 926 927 PetscFunctionBeginUser; 928 *inrd = 0; 929 PetscCall(PetscNew(&rd)); 930 931 PetscOptionsBegin(comm, NULL, "Options for nonequilibrium radiation-diffusion with RD ionization", NULL); 932 { 933 rd->initial = 1; 934 PetscCall(PetscOptionsInt("-rd_initial", "Initial condition (1=Marshak, 2=Blast, 3=Marshak+)", "", rd->initial, &rd->initial, 0)); 935 switch (rd->initial) { 936 case 1: 937 case 2: 938 rd->unit.kilogram = 1.; 939 rd->unit.meter = 1.; 940 rd->unit.second = 1.; 941 rd->unit.Kelvin = 1.; 942 break; 943 case 3: 944 rd->unit.kilogram = 1.e12; 945 rd->unit.meter = 1.; 946 rd->unit.second = 1.e9; 947 rd->unit.Kelvin = 1.; 948 break; 949 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unknown initial condition %" PetscInt_FMT, rd->initial); 950 } 951 /* Fundamental units */ 952 PetscCall(PetscOptionsReal("-rd_unit_meter", "Length of 1 meter in nondimensional units", "", rd->unit.meter, &rd->unit.meter, 0)); 953 PetscCall(PetscOptionsReal("-rd_unit_kilogram", "Mass of 1 kilogram in nondimensional units", "", rd->unit.kilogram, &rd->unit.kilogram, 0)); 954 PetscCall(PetscOptionsReal("-rd_unit_second", "Time of a second in nondimensional units", "", rd->unit.second, &rd->unit.second, 0)); 955 PetscCall(PetscOptionsReal("-rd_unit_Kelvin", "Temperature of a Kelvin in nondimensional units", "", rd->unit.Kelvin, &rd->unit.Kelvin, 0)); 956 /* Derived units */ 957 rd->unit.Joule = rd->unit.kilogram * PetscSqr(rd->unit.meter / rd->unit.second); 958 rd->unit.Watt = rd->unit.Joule / rd->unit.second; 959 /* Local aliases */ 960 meter = rd->unit.meter; 961 kilogram = rd->unit.kilogram; 962 second = rd->unit.second; 963 Kelvin = rd->unit.Kelvin; 964 Joule = rd->unit.Joule; 965 Watt = rd->unit.Watt; 966 967 PetscCall(PetscOptionsBool("-rd_monitor_residual", "Display residuals every time they are evaluated", "", rd->monitor_residual, &rd->monitor_residual, NULL)); 968 PetscCall(PetscOptionsEnum("-rd_discretization", "Discretization type", "", DiscretizationTypes, (PetscEnum)rd->discretization, (PetscEnum *)&rd->discretization, NULL)); 969 if (rd->discretization == DISCRETIZATION_FE) { 970 rd->quadrature = QUADRATURE_GAUSS2; 971 PetscCall(PetscOptionsEnum("-rd_quadrature", "Finite element quadrature", "", QuadratureTypes, (PetscEnum)rd->quadrature, (PetscEnum *)&rd->quadrature, NULL)); 972 } 973 PetscCall(PetscOptionsEnum("-rd_jacobian", "Type of finite difference Jacobian", "", JacobianTypes, (PetscEnum)rd->jacobian, (PetscEnum *)&rd->jacobian, NULL)); 974 switch (rd->initial) { 975 case 1: 976 rd->leftbc = BC_ROBIN; 977 rd->Eapplied = 4 * rd->unit.Joule / PetscPowRealInt(rd->unit.meter, 3); 978 rd->L = 1. * rd->unit.meter; 979 rd->beta = 3.0; 980 rd->gamma = 3.0; 981 rd->final_time = 3 * second; 982 break; 983 case 2: 984 rd->leftbc = BC_NEUMANN; 985 rd->Eapplied = 0.; 986 rd->L = 1. * rd->unit.meter; 987 rd->beta = 3.0; 988 rd->gamma = 3.0; 989 rd->final_time = 1 * second; 990 break; 991 case 3: 992 rd->leftbc = BC_ROBIN; 993 rd->Eapplied = 7.503e6 * rd->unit.Joule / PetscPowRealInt(rd->unit.meter, 3); 994 rd->L = 5. * rd->unit.meter; 995 rd->beta = 3.5; 996 rd->gamma = 3.5; 997 rd->final_time = 20e-9 * second; 998 break; 999 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Initial %" PetscInt_FMT, rd->initial); 1000 } 1001 PetscCall(PetscOptionsEnum("-rd_leftbc", "Left boundary condition", "", BCTypes, (PetscEnum)rd->leftbc, (PetscEnum *)&rd->leftbc, NULL)); 1002 PetscCall(PetscOptionsReal("-rd_E_applied", "Radiation flux at left end of domain", "", rd->Eapplied, &rd->Eapplied, NULL)); 1003 PetscCall(PetscOptionsReal("-rd_beta", "Thermal exponent for photon absorption", "", rd->beta, &rd->beta, NULL)); 1004 PetscCall(PetscOptionsReal("-rd_gamma", "Thermal exponent for diffusion coefficient", "", rd->gamma, &rd->gamma, NULL)); 1005 PetscCall(PetscOptionsBool("-rd_view_draw", "Draw final solution", "", rd->view_draw, &rd->view_draw, NULL)); 1006 PetscCall(PetscOptionsBool("-rd_endpoint", "Discretize using endpoints (like trapezoid rule) instead of midpoint", "", rd->endpoint, &rd->endpoint, NULL)); 1007 PetscCall(PetscOptionsBool("-rd_bcmidpoint", "Impose the boundary condition at the midpoint (Theta) of the interval", "", rd->bcmidpoint, &rd->bcmidpoint, NULL)); 1008 PetscCall(PetscOptionsBool("-rd_bclimit", "Limit diffusion coefficient in definition of Robin boundary condition", "", rd->bclimit, &rd->bclimit, NULL)); 1009 PetscCall(PetscOptionsBool("-rd_test_diff", "Test differentiation in constitutive relations", "", rd->test_diff, &rd->test_diff, NULL)); 1010 PetscCall(PetscOptionsString("-rd_view_binary", "File name to hold final solution", "", rd->view_binary, rd->view_binary, sizeof(rd->view_binary), NULL)); 1011 } 1012 PetscOptionsEnd(); 1013 1014 switch (rd->initial) { 1015 case 1: 1016 case 2: 1017 rd->rho = 1.; 1018 rd->c = 1.; 1019 rd->K_R = 1.; 1020 rd->K_p = 1.; 1021 rd->sigma_b = 0.25; 1022 rd->MaterialEnergy = RDMaterialEnergy_Reduced; 1023 break; 1024 case 3: 1025 /* Table 2 */ 1026 rd->rho = 1.17e-3 * kilogram / (meter * meter * meter); /* density */ 1027 rd->K_R = 7.44e18 * PetscPowRealInt(meter, 5) * PetscPowReal(Kelvin, 3.5) * PetscPowRealInt(kilogram, -2); /* */ 1028 rd->K_p = 2.33e20 * PetscPowRealInt(meter, 5) * PetscPowReal(Kelvin, 3.5) * PetscPowRealInt(kilogram, -2); /* */ 1029 rd->I_H = 2.179e-18 * Joule; /* Hydrogen ionization potential */ 1030 rd->m_p = 1.673e-27 * kilogram; /* proton mass */ 1031 rd->m_e = 9.109e-31 * kilogram; /* electron mass */ 1032 rd->h = 6.626e-34 * Joule * second; /* Planck's constant */ 1033 rd->k = 1.381e-23 * Joule / Kelvin; /* Boltzman constant */ 1034 rd->c = 3.00e8 * meter / second; /* speed of light */ 1035 rd->sigma_b = 5.67e-8 * Watt * PetscPowRealInt(meter, -2) * PetscPowRealInt(Kelvin, -4); /* Stefan-Boltzman constant */ 1036 rd->MaterialEnergy = RDMaterialEnergy_Saha; 1037 break; 1038 } 1039 1040 PetscCall(DMDACreate1d(comm, DM_BOUNDARY_NONE, 20, sizeof(RDNode) / sizeof(PetscScalar), 1, NULL, &rd->da)); 1041 PetscCall(DMSetFromOptions(rd->da)); 1042 PetscCall(DMSetUp(rd->da)); 1043 PetscCall(DMDASetFieldName(rd->da, 0, "E")); 1044 PetscCall(DMDASetFieldName(rd->da, 1, "T")); 1045 PetscCall(DMDASetUniformCoordinates(rd->da, 0., 1., 0., 0., 0., 0.)); 1046 1047 *inrd = rd; 1048 PetscFunctionReturn(0); 1049 } 1050 1051 int main(int argc, char *argv[]) { 1052 RD rd; 1053 TS ts; 1054 SNES snes; 1055 Vec X; 1056 Mat A, B; 1057 PetscInt steps; 1058 PetscReal ftime; 1059 1060 PetscFunctionBeginUser; 1061 PetscCall(PetscInitialize(&argc, &argv, 0, NULL)); 1062 PetscCall(RDCreate(PETSC_COMM_WORLD, &rd)); 1063 PetscCall(DMCreateGlobalVector(rd->da, &X)); 1064 PetscCall(DMSetMatType(rd->da, MATAIJ)); 1065 PetscCall(DMCreateMatrix(rd->da, &B)); 1066 PetscCall(RDInitialState(rd, X)); 1067 1068 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 1069 PetscCall(TSSetProblemType(ts, TS_NONLINEAR)); 1070 PetscCall(TSSetType(ts, TSTHETA)); 1071 PetscCall(TSSetDM(ts, rd->da)); 1072 switch (rd->discretization) { 1073 case DISCRETIZATION_FD: 1074 PetscCall(TSSetIFunction(ts, NULL, RDIFunction_FD, rd)); 1075 if (rd->jacobian == JACOBIAN_ANALYTIC) PetscCall(TSSetIJacobian(ts, B, B, RDIJacobian_FD, rd)); 1076 break; 1077 case DISCRETIZATION_FE: 1078 PetscCall(TSSetIFunction(ts, NULL, RDIFunction_FE, rd)); 1079 if (rd->jacobian == JACOBIAN_ANALYTIC) PetscCall(TSSetIJacobian(ts, B, B, RDIJacobian_FE, rd)); 1080 break; 1081 } 1082 PetscCall(TSSetMaxTime(ts, rd->final_time)); 1083 PetscCall(TSSetTimeStep(ts, 1e-3)); 1084 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); 1085 PetscCall(TSSetFromOptions(ts)); 1086 1087 A = B; 1088 PetscCall(TSGetSNES(ts, &snes)); 1089 switch (rd->jacobian) { 1090 case JACOBIAN_ANALYTIC: break; 1091 case JACOBIAN_MATRIXFREE: break; 1092 case JACOBIAN_FD_COLORING: { 1093 PetscCall(SNESSetJacobian(snes, A, B, SNESComputeJacobianDefaultColor, 0)); 1094 } break; 1095 case JACOBIAN_FD_FULL: PetscCall(SNESSetJacobian(snes, A, B, SNESComputeJacobianDefault, ts)); break; 1096 } 1097 1098 if (rd->test_diff) PetscCall(RDTestDifferentiation(rd)); 1099 PetscCall(TSSolve(ts, X)); 1100 PetscCall(TSGetSolveTime(ts, &ftime)); 1101 PetscCall(TSGetStepNumber(ts, &steps)); 1102 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Steps %" PetscInt_FMT " final time %g\n", steps, (double)ftime)); 1103 if (rd->view_draw) PetscCall(RDView(rd, X, PETSC_VIEWER_DRAW_WORLD)); 1104 if (rd->view_binary[0]) { 1105 PetscViewer viewer; 1106 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, rd->view_binary, FILE_MODE_WRITE, &viewer)); 1107 PetscCall(RDView(rd, X, viewer)); 1108 PetscCall(PetscViewerDestroy(&viewer)); 1109 } 1110 PetscCall(VecDestroy(&X)); 1111 PetscCall(MatDestroy(&B)); 1112 PetscCall(RDDestroy(&rd)); 1113 PetscCall(TSDestroy(&ts)); 1114 PetscCall(PetscFinalize()); 1115 return 0; 1116 } 1117 /*TEST 1118 1119 test: 1120 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 1121 requires: !single 1122 1123 test: 1124 suffix: 2 1125 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 1126 requires: !single 1127 1128 test: 1129 suffix: 3 1130 nsize: 2 1131 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 1132 requires: !single 1133 1134 TEST*/ 1135