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