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