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