Actual source code: ex2.c

petsc-master 2020-10-26
Report Typos and Errors
  1: static char help[] = "Runaway electron model with Landau collision operator\n\n";

  3: #include <petscdmplex.h>
  4: #include <petsclandau.h>
  5: #include <petscts.h>
  6: #include <petscds.h>

  8: /* data for runaway electron model */
  9: typedef struct REctx_struct {
 10:   PetscErrorCode (*test)(TS, Vec, DM, PetscInt, PetscReal, PetscBool,  LandauCtx *, struct REctx_struct *);
 11:   PetscErrorCode (*impuritySrcRate)(PetscReal, PetscReal *, LandauCtx*);
 12:   PetscErrorCode (*E)(Vec, Vec, PetscInt, PetscReal, LandauCtx*, PetscReal *);
 13:   PetscReal     T_cold;        /* temperature of newly ionized electrons and impurity ions */
 14:   PetscReal     ion_potential; /* ionization potential of impurity */
 15:   PetscReal     Ne_ion;        /* effective number of electrons shed in ioization of impurity */
 16:   PetscReal     Ez_initial;
 17:   PetscReal     L;             /* inductance */
 18:   Vec           X_0;
 19:   PetscInt      imp_idx;       /* index for impurity ionizing sink */
 20:   PetscReal     pulse_start;
 21:   PetscReal     pulse_width;
 22:   PetscReal     pulse_rate;
 23:   PetscReal     current_rate;
 24:   PetscInt      plotIdx;
 25:   PetscInt      plotStep;
 26:   PetscInt      idx; /* cache */
 27:   PetscReal     j; /* cache */
 28:   PetscReal     plotDt;
 29:   PetscBool     plotting;
 30:   PetscBool     use_spitzer_eta;
 31:   Vec           imp_src;
 32: } REctx;

 34: static const PetscReal kev_joul = 6.241506479963235e+15; /* 1/1000e */
 35: static PetscBool s_quarter3DDomain_notused = PETSC_FALSE;

 37: #define RE_CUT 5.
 38: /* < v, u_re * v * q > */
 39: static void f0_j_re(PetscInt dim, PetscInt Nf, PetscInt NfAux,
 40:                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
 41:                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
 42:                     PetscReal t, const PetscReal x[],  PetscInt numConstants, const PetscScalar constants[], PetscScalar *f0)
 43: {
 44:   PetscReal n_e = PetscRealPart(u[0]);
 45:   if (dim==2) {
 46:     if (x[1] > RE_CUT || x[1] < -RE_CUT) { /* simply a cutoff for REs. v_|| > 3 v(T_e) */
 47:       *f0 = n_e * 2.*PETSC_PI*x[0] * x[1] * constants[0]; /* n * r * v_|| * q */
 48:     } else {
 49:       *f0 = 0;
 50:     }
 51:   } else {
 52:     if (x[2] > RE_CUT || x[2] < -RE_CUT) { /* simply a cutoff for REs. v_|| > 3 v(T_e) */
 53:       *f0 = n_e                * x[2] * constants[0];
 54:       if (s_quarter3DDomain_notused) *f0 *= 4.0;
 55:     } else {
 56:       *f0 = 0;
 57:     }
 58:   }
 59: }

 61: /* sum < v, u*v*q > */
 62: static void f0_jz_sum(PetscInt dim, PetscInt Nf, PetscInt NfAux,
 63:                    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
 64:                    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
 65:                    PetscReal t, const PetscReal x[],  PetscInt numConstants, const PetscScalar constants[], PetscScalar *f0)
 66: {
 67:   PetscInt ii;
 68:   f0[0] = 0;
 69:   if (dim==2) {
 70:     for (ii=0;ii<numConstants;ii++) f0[0] += u[ii] * 2.*PETSC_PI*x[0] * x[1] * constants[ii]; /* n * r * v_|| * q */
 71:   } else {
 72:     for (ii=0;ii<numConstants;ii++) f0[0] += u[ii]                * x[2] * constants[ii]; /* n * v_|| * q  */
 73:     if (s_quarter3DDomain_notused) *f0 *= 4.0;
 74:   }
 75: }

 77: /* < v, n_e > */
 78: static void f0_n(PetscInt dim, PetscInt Nf, PetscInt NfAux,
 79:                   const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
 80:                   const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
 81:                   PetscReal t, const PetscReal x[],  PetscInt numConstants, const PetscScalar constants[], PetscScalar *f0)
 82: {
 83:   PetscInt ii = (PetscInt)PetscRealPart(constants[0]);
 84:   if (dim==2) f0[0] = 2.*PETSC_PI*x[0]*u[ii];
 85:   else {
 86:     f0[0] =                        u[ii];
 87:     if (s_quarter3DDomain_notused) f0[0] *= 4.0;
 88:   }
 89: }

 91: /* < v, n_e v_|| > */
 92: static void f0_vz(PetscInt dim, PetscInt Nf, PetscInt NfAux,
 93:                    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
 94:                    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
 95:                    PetscReal t, const PetscReal x[],  PetscInt numConstants, const PetscScalar constants[], PetscScalar *f0)
 96: {
 97:   PetscInt ii = (PetscInt)PetscRealPart(constants[0]);
 98:   if (dim==2) f0[0] = u[ii] * 2.*PETSC_PI*x[0] * x[1]; /* n r v_|| */
 99:   else {
100:     f0[0] =           u[ii] *                x[2]; /* n v_|| */
101:     if (s_quarter3DDomain_notused) f0[0] *= 4.0;
102:   }
103: }

105: /* < v, n_e (v-shift) > */
106: static void f0_ve_shift(PetscInt dim, PetscInt Nf, PetscInt NfAux,
107:                          const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
108:                          const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
109:                          PetscReal t, const PetscReal x[],  PetscInt numConstants, const PetscScalar constants[], PetscScalar *f0)
110: {
111:   PetscReal vz = numConstants==1 ? PetscRealPart(constants[0]) : 0;
112:   if (dim==2) *f0 = u[0] * 2.*PETSC_PI*x[0] * PetscSqrtReal(x[0]*x[0] + (x[1]-vz)*(x[1]-vz));             /* n r v */
113:   else {
114:     *f0 =           u[0] *                PetscSqrtReal(x[0]*x[0] + x[1]*x[1] + (x[2]-vz)*(x[2]-vz)); /* n v */
115:     if (s_quarter3DDomain_notused) *f0 *= 4.0;
116:   }
117: }

119: static PetscErrorCode getTe_kev(DM plex, Vec X, PetscReal *a_n, PetscReal *a_Tkev)
120: {
122:   LandauCtx      *ctx;

125:   DMGetApplicationContext(plex, &ctx);
126:   {
127:     PetscDS        prob;
128:     PetscReal      v2, v, n;
129:     PetscScalar    tt[LANDAU_MAX_SPECIES],user[2] = {0.,ctx->charges[0]}, vz;
130:     DMGetDS(plex, &prob);
131:     PetscDSSetConstants(prob, 2, user);
132:     PetscDSSetObjective(prob, 0, &f0_n);
133:     DMPlexComputeIntegralFEM(plex,X,tt,NULL);
134:     n = ctx->n_0*PetscRealPart(tt[0]);
135:     PetscDSSetObjective(prob, 0, &f0_vz);
136:     DMPlexComputeIntegralFEM(plex,X,tt,NULL);
137:     vz = ctx->n_0*PetscRealPart(tt[0])/n; /* non-dimensional */
138:     /* remove drift */
139:     PetscDSSetConstants(prob, 1, &vz);
140:     PetscDSSetObjective(prob, 0, &f0_ve_shift);
141:     DMPlexComputeIntegralFEM(plex,X,tt,NULL);
142:     v = ctx->n_0*ctx->v_0*PetscRealPart(tt[0])/n;         /* remove number density to get velocity */
143:     v2 = PetscSqr(v);                      /* use real space: m^2 / s^2 */
144:     if (a_Tkev) *a_Tkev = (v2*ctx->masses[0]*PETSC_PI/8)*kev_joul; /* temperature in kev */
145:     if (a_n) *a_n = n;
146:   }
147:   return(0);
148: }
149:  /* CalculateE - Calculate the electric field  */
150:  /*  T        -- Electron temperature  */
151:  /*  n        -- Electron density  */
152:  /*  lnLambda --   */
153:  /*  eps0     --  */
154:  /*  E        -- output E, input \hat E */
155: static PetscReal CalculateE(PetscReal Tev, PetscReal n, PetscReal lnLambda, PetscReal eps0, PetscReal *E)
156: {
157:   PetscReal c,e,m;

160:   c = 299792458;
161:   e = 1.602176e-19;
162:   m = 9.10938e-31;
163:   if (1) {
164:     double Ec, Ehat = *E, betath = PetscSqrtReal(2*Tev*e/(m*c*c)), j0 = Ehat * 7/(PetscSqrtReal(2)*2) * PetscPowReal(betath,3) * n * e * c;
165:     Ec = n*lnLambda*PetscPowReal(e,3) / (4*PETSC_PI*PetscPowReal(eps0,2)*m*c*c);
166:     *E = Ec;
167:     PetscPrintf(PETSC_COMM_WORLD, "CalculateE j0=%g Ec = %g\n",j0,Ec);
168:   } else {
169:     PetscReal Ed, vth;
170:     vth = PetscSqrtReal(8*Tev*e/(m*PETSC_PI));
171:     Ed =  n*lnLambda*PetscPowReal(e,3) / (4*PETSC_PI*PetscPowReal(eps0,2)*m*vth*vth);
172:     *E = Ed;
173:   }
174:   return(0);
175: }

177: static PetscReal Spitzer(PetscReal m_e, PetscReal e, PetscReal Z, PetscReal epsilon0,  PetscReal lnLam, PetscReal kTe_joules)
178: {
179:   PetscReal Fz = (1+1.198*Z+0.222*Z*Z)/(1+2.966*Z+0.753*Z*Z), eta;
180:   eta = Fz*4./3.*PetscSqrtReal(2.*PETSC_PI)*Z*PetscSqrtReal(m_e)*PetscSqr(e)*lnLam*PetscPowReal(4*PETSC_PI*epsilon0,-2.)*PetscPowReal(kTe_joules,-1.5);
181:   return eta;
182: }

184: /*  */
185: static PetscErrorCode testNone(TS ts, Vec X, DM plex, PetscInt stepi, PetscReal time, PetscBool islast, LandauCtx *ctx, REctx *rectx)
186: {
188:   return(0);
189: }

191: /*  */
192: static PetscErrorCode testSpitzer(TS ts, Vec X, DM plex, PetscInt stepi, PetscReal time, PetscBool islast, LandauCtx *ctx, REctx *rectx)
193: {
194:   PetscErrorCode    ierr;
195:   PetscInt          ii;
196:   PetscDS           prob;
197:   static PetscReal  old_ratio = 0;
198:   PetscBool         done=PETSC_FALSE;
199:   PetscReal         J,J_re,spit_eta,Te_kev=0,E,ratio,Z,n_e;
200:   PetscScalar       user[2] = {0.,ctx->charges[0]}, constants[LANDAU_MAX_SPECIES],tt[LANDAU_MAX_SPECIES];

203:   if (ctx->num_species<2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "ctx->num_species %D < 2",ctx->num_species);
204:   for (ii=0;ii<ctx->num_species;ii++) constants[ii] = ctx->charges[ii];
205:   Z = -ctx->charges[1]/ctx->charges[0];
206:   DMGetDS(plex, &prob);
207:   PetscDSSetConstants(prob, 2, user);
208:   PetscDSSetObjective(prob, 0, &f0_n);
209:   DMPlexComputeIntegralFEM(plex,X,tt,NULL);
210:   n_e = PetscRealPart(tt[0])*ctx->n_0;
211:   PetscDSSetConstants(prob, ctx->num_species, constants);
212:   PetscDSSetObjective(prob, 0, &f0_jz_sum);
213:   DMPlexComputeIntegralFEM(plex,X,tt,NULL);
214:   J = -ctx->n_0*ctx->v_0*PetscRealPart(tt[0]);
215:   PetscDSSetConstants(prob, 1, &constants[0]);
216:   PetscDSSetObjective(prob, 0, &f0_j_re);
217:   DMPlexComputeIntegralFEM(plex,X,tt,NULL);
218:   J_re = -ctx->n_0*ctx->v_0*PetscRealPart(tt[0]);
219:   getTe_kev(plex, X, NULL, &Te_kev);
220:   spit_eta = Spitzer(ctx->masses[0],-ctx->charges[0],Z,ctx->epsilon0,ctx->lnLam,Te_kev/kev_joul); /* kev --> J (kT) */
221:   E = ctx->Ez; /* keep real E */
222:   ratio = E/J/spit_eta;
223:   done = PETSC_FALSE;
224:   PetscPrintf(PETSC_COMM_WORLD, "%s %4D) time=%10.3e n_e= %10.3e E= %10.3e J= %10.3e J_re= %10.3e %.3g %% Te_kev= %10.3e E/J to eta ratio=%g (diff=%g)\n",
225:                      done ? "DONE" : "testSpitzer",stepi,time,n_e/ctx->n_0,E,J,J_re,100*J_re/J,Te_kev,ratio,old_ratio-ratio);
226:   if (done) {
227:     TSSetConvergedReason(ts,TS_CONVERGED_USER);
228:     old_ratio = 0;
229:   } else {
230:     TSConvergedReason reason;
231:     TSGetConvergedReason(ts,&reason);
232:     old_ratio = ratio;
233:     if (reason) done = PETSC_TRUE;
234:   }
235:   return(0);
236: }

238: static const double ppp = 2;
239: static void f0_0_diff_lp(PetscInt dim, PetscInt Nf, PetscInt NfAux,
240:                           const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
241:                           const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
242:                           PetscReal t, const PetscReal x[],  PetscInt numConstants, const PetscScalar constants[], PetscScalar *f0)
243: {
244:   LandauCtx       *ctx = (LandauCtx *)constants;
245:   REctx           *rectx = (REctx*)ctx->data;
246:   PetscInt        ii = rectx->idx, i;
247:   const PetscReal kT_m = ctx->k*ctx->thermal_temps[ii]/ctx->masses[ii]; /* kT/m */
248:   const PetscReal n = ctx->n[ii];
249:   PetscReal       diff, f_maxwell, v2 = 0, theta = 2*kT_m/(ctx->v_0*ctx->v_0); /* theta = 2kT/mc^2 */
250:   for (i = 0; i < dim; ++i) v2 += x[i]*x[i];
251:   f_maxwell = n*PetscPowReal(PETSC_PI*theta,-1.5)*(PetscExpReal(-v2/theta));
252:   diff = 2.*PETSC_PI*x[0]*(PetscRealPart(u[ii]) - f_maxwell);
253:   f0[0] = PetscPowReal(diff,ppp);
254: }
255: static void f0_0_maxwellian_lp(PetscInt dim, PetscInt Nf, PetscInt NfAux,
256:                           const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
257:                           const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
258:                           PetscReal t, const PetscReal x[],  PetscInt numConstants, const PetscScalar constants[], PetscScalar *f0)
259: {
260:   LandauCtx       *ctx = (LandauCtx *)constants;
261:   REctx           *rectx = (REctx*)ctx->data;
262:   PetscInt        ii = rectx->idx, i;
263:   const PetscReal kT_m = ctx->k*ctx->thermal_temps[ii]/ctx->masses[ii]; /* kT/m */
264:   const PetscReal n = ctx->n[ii];
265:   PetscReal       f_maxwell, v2 = 0, theta = 2*kT_m/(ctx->v_0*ctx->v_0); /* theta = 2kT/mc^2 */
266:   for (i = 0; i < dim; ++i) v2 += x[i]*x[i];
267:   f_maxwell = 2.*PETSC_PI*x[0] * n*PetscPowReal(PETSC_PI*theta,-1.5)*(PetscExpReal(-v2/theta));
268:   f0[0] = PetscPowReal(f_maxwell,ppp);
269: }

271: /*  */
272: static PetscErrorCode testStable(TS ts, Vec X, DM plex, PetscInt stepi, PetscReal time, PetscBool islast, LandauCtx *ctx, REctx *rectx)
273: {
274:   PetscErrorCode    ierr;
275:   PetscDS           prob;
276:   Vec               X2;
277:   PetscReal         ediff,idiff=0,lpm0,lpm1=1;
278:   PetscScalar       tt[LANDAU_MAX_SPECIES];
279:   DM                dm;

282:   VecGetDM(X, &dm);
283:   DMGetDS(plex, &prob);
284:   VecDuplicate(X,&X2);
285:   VecCopy(X,X2);
286:   if (!rectx->X_0) {
287:     VecDuplicate(X,&rectx->X_0);
288:     VecCopy(X,rectx->X_0);
289:   }
290:   VecAXPY(X,-1.0,rectx->X_0);
291:   PetscDSSetConstants(prob, sizeof(LandauCtx)/sizeof(PetscScalar), (PetscScalar*)ctx);
292:   rectx->idx = 0;
293:   PetscDSSetObjective(prob, 0, &f0_0_diff_lp);
294:   DMPlexComputeIntegralFEM(plex,X2,tt,NULL);
295:   ediff = PetscPowReal(PetscRealPart(tt[0]),1./ppp);
296:   PetscDSSetObjective(prob, 0, &f0_0_maxwellian_lp);
297:   DMPlexComputeIntegralFEM(plex,X2,tt,NULL);
298:   lpm0 = PetscPowReal(PetscRealPart(tt[0]),1./ppp);
299:   if (ctx->num_species>1) {
300:     rectx->idx = 1;
301:     PetscDSSetObjective(prob, 0, &f0_0_diff_lp);
302:     DMPlexComputeIntegralFEM(plex,X2,tt,NULL);
303:     idiff = PetscPowReal(PetscRealPart(tt[0]),1./ppp);
304:     PetscDSSetObjective(prob, 0, &f0_0_maxwellian_lp);
305:     DMPlexComputeIntegralFEM(plex,X2,tt,NULL);
306:     lpm1 = PetscPowReal(PetscRealPart(tt[0]),1./ppp);
307:   }
308:   PetscPrintf(PETSC_COMM_WORLD, "%s %D) time=%10.3e n-%d norm electrons/max=%20.13e ions/max=%20.13e\n", "----",stepi,time,(int)ppp,ediff/lpm0,idiff/lpm1);
309:   /* view */
310:   VecViewFromOptions(X,NULL,"-vec_view_diff");
311:   VecCopy(X2,X);
312:   VecDestroy(&X2);
313:   if (islast) {
314:     VecDestroy(&rectx->X_0);
315:     rectx->X_0 = NULL;
316:   }
317:   return(0);
318: }

320: /* E = eta_spitzer(Te-vz)*J */
321: static PetscErrorCode ESpitzer(Vec X,  Vec X_t,  PetscInt stepi, PetscReal time, LandauCtx *ctx, PetscReal *a_E)
322: {
323:   PetscErrorCode    ierr;
324:   PetscInt          ii;
325:   PetscReal         spit_eta,Te_kev,J,ratio;
326:   PetscScalar       tt[LANDAU_MAX_SPECIES], constants[LANDAU_MAX_SPECIES];
327:   PetscDS           prob;
328:   DM                dm,plex;
329:   REctx             *rectx = (REctx*)ctx->data;

332:   for (ii=0;ii<ctx->num_species;ii++) constants[ii] = ctx->charges[ii];
333:   VecGetDM(X, &dm);
334:   DMConvert(dm, DMPLEX, &plex);
335:   DMGetDS(plex, &prob);
336:   /* J */
337:   PetscDSSetConstants(prob, ctx->num_species, constants);
338:   PetscDSSetObjective(prob, 0, &f0_jz_sum);
339:   DMPlexComputeIntegralFEM(plex,X,tt,NULL);
340:   J = -ctx->n_0*ctx->v_0*PetscRealPart(tt[0]);
341:   getTe_kev(plex, X, NULL, &Te_kev);
342:   spit_eta = Spitzer(ctx->masses[0],-ctx->charges[0],-ctx->charges[1]/ctx->charges[0],ctx->epsilon0,ctx->lnLam,Te_kev/kev_joul); /* kev --> J (kT) */
343:   *a_E = ctx->Ez; /* no change */
344:   if (!rectx->use_spitzer_eta && time > 10) {
345:     static PetscReal  old_ratio = 1e10;
346:     ratio = *a_E/J/spit_eta;
347:     if ((ratio < 1.01 && ratio > 0.99) || (old_ratio <= ratio && ratio < 1.03 && ratio > 0.97)) {
348:       rectx->use_spitzer_eta = PETSC_TRUE; /* use it next time */
349:       rectx->j = J;
350:       rectx->pulse_start = time + 1.; /* start quench now */
351:     }
352:     PetscPrintf(PETSC_COMM_WORLD,"\t\t%D) t=%10.3e ESpitzer E/J vs spitzer ratio=%20.13e J=%10.3e E=%10.3e spit_eta=%10.3e Te_kev=%10.3e %s\n",stepi,time,ratio, J, *a_E, spit_eta, Te_kev, rectx->use_spitzer_eta ? " switch to Spitzer E" : " keep testing");
353:     old_ratio = ratio;
354:   } else if (rectx->use_spitzer_eta) {
355:     /* set E */
356:     *a_E = spit_eta*J;
357:     PetscPrintf(PETSC_COMM_WORLD,"\t%D) use ESpitzer E=%10.3e J=%10.3e Te_kev=%10.3e spit_eta=%10.3e t=%g\n",stepi,*a_E,J,Te_kev,spit_eta,time);
358:   } else {
359:     PetscPrintf(PETSC_COMM_WORLD,"\t\t\t%D) ESpitzer delay E=%10.3e J=%10.3e Te_kev=%10.3e spit_eta=%10.3e t=%g ratio=%20.13e\n",stepi,*a_E,J,Te_kev,spit_eta,time,ctx->Ez/J/spit_eta);
360:   }
361:   /* cleanup */
362:   DMDestroy(&plex);
363:   return(0);
364: }

366: static PetscErrorCode EInduction(Vec X, Vec X_t, PetscInt step, PetscReal time, LandauCtx *ctx, PetscReal *a_E)
367: {
368:   REctx             *rectx = (REctx*)ctx->data;
369:   PetscErrorCode    ierr;
370:   PetscInt          ii;
371:   DM                dm,plex;
372:   PetscScalar       tt[LANDAU_MAX_SPECIES], constants[LANDAU_MAX_SPECIES];
373:   PetscReal         dJ_dt;
374:   PetscDS           prob;

377:   for (ii=0;ii<ctx->num_species;ii++) constants[ii] = ctx->charges[ii];
378:   VecGetDM(X, &dm);
379:   DMGetDS(dm, &prob);
380:   DMConvert(dm, DMPLEX, &plex);
381:   /* get d current / dt */
382:   PetscDSSetConstants(prob, ctx->num_species, constants);
383:   PetscDSSetObjective(prob, 0, &f0_jz_sum);
384:   if (!X_t) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "X_t");
385:   DMPlexComputeIntegralFEM(plex,X_t,tt,NULL);
386:   dJ_dt = -ctx->n_0*ctx->v_0*PetscRealPart(tt[0])/ctx->t_0;
387:   /* E induction */
388:   *a_E = -rectx->L*dJ_dt + rectx->Ez_initial;
389:   DMDestroy(&plex);
390:   return(0);
391: }

393: static PetscErrorCode EConstant(Vec X,  Vec X_t, PetscInt step, PetscReal time, LandauCtx *ctx, PetscReal *a_E)
394: {
396:   *a_E = ctx->Ez;
397:   return(0);
398: }

400: static PetscErrorCode ENone(Vec X,  Vec X_t, PetscInt step, PetscReal time, LandauCtx *ctx, PetscReal *a_E)
401: {
403:   *a_E = 0;
404:   return(0);
405: }

407: /* ------------------------------------------------------------------- */
408: /*
409:    FormSource - Evaluates source terms F(t).

411:    Input Parameters:
412: .  ts - the TS context
413: .  time -
414: .  X_dummmy - input vector
415: .  dummy - optional user-defined context, as set by SNESSetFunction()

417:    Output Parameter:
418: .  F - function vector
419:  */
420: static PetscErrorCode FormSource(TS ts,PetscReal ftime,Vec X_dummmy, Vec F,void *dummy)
421: {
422:   PetscReal      new_imp_rate;
423:   LandauCtx      *ctx;
424:   DM             dm,plex;
426:   REctx          *rectx;

429:   TSGetDM(ts,&dm);
430:   DMGetApplicationContext(dm, &ctx);
431:   rectx = (REctx*)ctx->data;
432:   /* check for impurities */
433:   rectx->impuritySrcRate(ftime,&new_imp_rate,ctx);
434:   if (new_imp_rate != 0) {
435:     if (new_imp_rate != rectx->current_rate) {
436:       PetscInt       ii;
437:       PetscReal      dne_dt,dni_dt,tilda_ns[LANDAU_MAX_SPECIES],temps[LANDAU_MAX_SPECIES];
438:       PetscDS        prob; /* diagnostics only */
439:       rectx->current_rate = new_imp_rate;
440:       DMConvert(dm, DMPLEX, &plex);
441:       DMGetDS(dm, &prob);
442:       dni_dt = new_imp_rate              /* *ctx->t_0 */; /* fully ionized immediately, no normalize, stay in non-dim */
443:       dne_dt = new_imp_rate*rectx->Ne_ion/* *ctx->t_0 */;
444:       PetscPrintf(PETSC_COMM_SELF, "\tFormSource: have new_imp_rate= %10.3e time= %10.3e de/dt= %10.3e di/dt= %10.3e ***\n",new_imp_rate,ftime,dne_dt,dni_dt);
445:       for (ii=1;ii<LANDAU_MAX_SPECIES;ii++) tilda_ns[ii] = 0;
446:       for (ii=1;ii<LANDAU_MAX_SPECIES;ii++)    temps[ii] = 1;
447:       tilda_ns[0] = dne_dt;        tilda_ns[rectx->imp_idx] = dni_dt;
448:       temps[0]    = rectx->T_cold;    temps[rectx->imp_idx] = rectx->T_cold;
449:       /* add it */
450:       if (!rectx->imp_src) {
451:         DMCreateGlobalVector(dm, &rectx->imp_src);
452:         PetscObjectSetName((PetscObject)rectx->imp_src, "source");
453:       }
454:       VecZeroEntries(rectx->imp_src);
455:       LandauAddMaxwellians(plex,rectx->imp_src,ftime,temps,tilda_ns,ctx);
456:       /* clean up */
457:       DMDestroy(&plex);
458:       if (0) {
459:         KSP ksp;
460:         Vec S;
461:         KSPCreate(PETSC_COMM_SELF,&ksp);
462:         KSPSetOptionsPrefix(ksp,"mass_"); /* stokes */
463:         KSPSetOperators(ksp,ctx->M,ctx->M);
464:         KSPSetFromOptions(ksp);
465:         DMCreateGlobalVector(dm, &S);
466:         KSPSolve(ksp,rectx->imp_src,S);
467:         VecCopy(S,rectx->imp_src);
468:         VecDestroy(&S);
469:       }
470:       VecViewFromOptions(rectx->imp_src,NULL,"-vec_view_sources");
471:     }
472:     VecCopy(rectx->imp_src,F);
473:   } else {
474:     if (rectx->current_rate != 0 && rectx->imp_src) {
475:       VecZeroEntries(rectx->imp_src);
476:     }
477:     rectx->current_rate = 0;
478:   }
479:   return(0);
480: }
481: PetscErrorCode Monitor(TS ts, PetscInt stepi, PetscReal time, Vec X, void *actx)
482: {
483:   LandauCtx         *ctx = (LandauCtx*) actx;   /* user-defined application context */
484:   REctx             *rectx = (REctx*)ctx->data;
485:   DM                dm,plex;
486:   PetscDS           prob;
487:   TSConvergedReason reason;
488:   PetscErrorCode    ierr;

491:   VecGetDM(X, &dm);
492:   if (stepi > rectx->plotStep && rectx->plotting) {
493:     rectx->plotting = PETSC_FALSE; /* was doing diagnostics, now done */
494:     rectx->plotIdx++;
495:   }
496:   /* view */
497:   TSGetConvergedReason(ts,&reason);
498:   if (time/rectx->plotDt >= (PetscReal)rectx->plotIdx || reason) {
499:     /* print norms */
500:     LandauPrintNorms(X, stepi);
501:     DMConvert(dm, DMPLEX, &plex);
502:     DMGetDS(plex, &prob);
503:     /* diagnostics */
504:     rectx->test(ts,X,plex,stepi,time,reason ? PETSC_TRUE : PETSC_FALSE, ctx,rectx);
505:     DMDestroy(&plex);
506:     /* view */
507:     DMSetOutputSequenceNumber(dm, rectx->plotIdx, time*ctx->t_0);
508:     VecViewFromOptions(X,NULL,"-vec_view");
509:     rectx->plotStep = stepi;
510:     rectx->plotting = PETSC_TRUE;
511:   }
512:   /* parallel check */
513:   if (reason) {
514:     PetscReal    val,rval;
515:     PetscMPIInt    rank;
516:     MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
517:     TSGetSolution(ts, &X);
518:     VecNorm(X,NORM_2,&val);
519:     MPIU_Allreduce(&val,&rval,1,MPIU_REAL,MPIU_MAX,PETSC_COMM_WORLD);
520:     if (rval != val) {
521:       PetscPrintf(PETSC_COMM_SELF, " ***** [%D] ERROR max |x| = %e, my |x| = %20.13e\n",rank,rval,val);
522:     } else {
523:       PetscPrintf(PETSC_COMM_WORLD, "[%D] parallel consistency check OK\n",rank);
524:     }
525:   }
526:   rectx->idx = 0;
527:   return(0);
528: }

530: PetscErrorCode PreStep(TS ts)
531: {
533:   LandauCtx      *ctx;
534:   REctx          *rectx;
535:   DM             dm;
536:   PetscInt       stepi;
537:   PetscReal      time;
538:   Vec            X;

541:   /* check seed RE run */
542:   TSGetDM(ts,&dm);
543:   TSGetTime(ts,&time);
544:   TSGetSolution(ts,&X);
545:   DMGetApplicationContext(dm, &ctx);
546:   rectx = (REctx*)ctx->data;
547:   TSGetStepNumber(ts, &stepi);
548:   /* update E */
549:   rectx->E(X, NULL, stepi, time, ctx, &ctx->Ez);
550:   return(0);
551: }

553: /* model for source of non-ionized impurities, profile provided by model, in du/dt form in normalized units (tricky because n_0 is normalized with electrons) */
554: static PetscErrorCode stepSrc(PetscReal time, PetscReal *rho, LandauCtx *ctx)
555: {
556:   REctx         *rectx = (REctx*)ctx->data;

559:   if (time >= rectx->pulse_start) *rho = rectx->pulse_rate;
560:   else *rho = 0.;
561:   return(0);
562: }
563: static PetscErrorCode zeroSrc(PetscReal time, PetscReal *rho, LandauCtx *ctx)
564: {
566:   *rho = 0.;
567:   return(0);
568: }
569: static PetscErrorCode pulseSrc(PetscReal time, PetscReal *rho, LandauCtx *ctx)
570: {
571:   REctx *rectx = (REctx*)ctx->data;

574:   if (time < rectx->pulse_start || time > rectx->pulse_start + 3*rectx->pulse_width) *rho = 0;
575:   else if (0) {
576:     double t = time - rectx->pulse_start, start = rectx->pulse_width, stop = 2*rectx->pulse_width, cycle = 3*rectx->pulse_width, steep = 5, xi = 0.75 - (stop - start)/(2* cycle);
577:     *rho = rectx->pulse_rate * (cycle / (stop - start)) / (1 + PetscExpReal(steep*(PetscSinReal(2*PETSC_PI*((t - start)/cycle + xi)) - PetscSinReal(2*PETSC_PI*xi))));
578:   } else if (0) {
579:     double x = 2*(time - rectx->pulse_start)/(3*rectx->pulse_width) - 1;
580:     if (x==1 || x==-1) *rho = 0;
581:     else *rho = rectx->pulse_rate * PetscExpReal(-1/(1-x*x));
582:   } else {
583:     double x = PetscSinReal((time-rectx->pulse_start)/(3*rectx->pulse_width)*2*PETSC_PI - PETSC_PI/2) + 1; /* 0:2, integrates to 1.0 */
584:     *rho = rectx->pulse_rate * x / (3*rectx->pulse_width);
585:   }
586:   return(0);
587: }

591: static PetscErrorCode ProcessREOptions(REctx *rectx, const LandauCtx *ctx, DM dm, const char prefix[])
592: {
593:   PetscErrorCode    ierr;
594:   PetscFunctionList plist = NULL, testlist = NULL, elist = NULL;
595:   char              pname[256],testname[256],ename[256];
596:   DM                dummy;

599:   DMCreate(PETSC_COMM_WORLD,&dummy);
600:   rectx->Ne_ion = 1;                 /* number of electrons given up by impurity ion */
601:   rectx->T_cold = .005;              /* kev */
602:   rectx->ion_potential = 15;         /* ev */
603:   rectx->L = 2;
604:   rectx->X_0 = NULL;
605:   rectx->imp_idx = ctx->num_species - 1; /* default ionized impurity as last one */
606:   rectx->pulse_start = 1;
607:   rectx->pulse_width = 1;
608:   rectx->plotStep = PETSC_MAX_INT;
609:   rectx->pulse_rate = 1.e-1;
610:   rectx->current_rate = 0;
611:   rectx->plotIdx = 0;
612:   rectx->imp_src = 0;
613:   rectx->j = 0;
614:   rectx->plotDt = 1.0;
615:   rectx->plotting = PETSC_TRUE;
616:   rectx->use_spitzer_eta = PETSC_FALSE;
617:   rectx->idx = 0;
618:   /* Register the available impurity sources */
619:   PetscFunctionListAdd(&plist,"step",&stepSrc);
620:   PetscFunctionListAdd(&plist,"none",&zeroSrc);
621:   PetscFunctionListAdd(&plist,"pulse",&pulseSrc);
622:   PetscStrcpy(pname,"none");
623:   PetscFunctionListAdd(&testlist,"none",&testNone);
624:   PetscFunctionListAdd(&testlist,"spitzer",&testSpitzer);
625:   PetscFunctionListAdd(&testlist,"stable",&testStable);
626:   PetscStrcpy(testname,"none");
627:   PetscFunctionListAdd(&elist,"none",&ENone);
628:   PetscFunctionListAdd(&elist,"spitzer",&ESpitzer);
629:   PetscFunctionListAdd(&elist,"induction",&EInduction);
630:   PetscFunctionListAdd(&elist,"constant",&EConstant);
631:   PetscStrcpy(ename,"constant");

633:   PetscOptionsBegin(PETSC_COMM_SELF, prefix, "Options for Runaway/seed electron model", "none");
634:   PetscOptionsReal("-ex2_plot_dt", "Plotting interval", "xgc_dmplex.c", rectx->plotDt, &rectx->plotDt, NULL);
635:   if (rectx->plotDt < 0) rectx->plotDt = 1e30;
636:   if (rectx->plotDt == 0) rectx->plotDt = 1e-30;
637:   PetscOptionsFList("-ex2_impurity_source_type","Name of impurity source to run","",plist,pname,pname,sizeof(pname),NULL);
638:   PetscOptionsFList("-ex2_test_type","Name of test to run","",testlist,testname,testname,sizeof(testname),NULL);
639:   PetscOptionsInt("-ex2_impurity_index", "index of sink for impurities", "none", rectx->imp_idx, &rectx->imp_idx, NULL);
640:   if ((rectx->imp_idx >= ctx->num_species || rectx->imp_idx < 1) && ctx->num_species > 1) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"index of sink for impurities ions is out of range (%D), must be > 0 && < NS",rectx->imp_idx);
641:   PetscOptionsFList("-ex2_e_field_type","Electric field type","",elist,ename,ename,sizeof(ename),NULL);
642:   rectx->Ne_ion = -ctx->charges[rectx->imp_idx]/ctx->charges[0];
643:   PetscOptionsReal("-ex2_t_cold","Temperature of cold electron and ions after ionization in keV","none",rectx->T_cold,&rectx->T_cold, NULL);
644:   PetscOptionsReal("-ex2_pulse_start_time","Time at which pulse happens for 'pulse' source","none",rectx->pulse_start,&rectx->pulse_start, NULL);
645:   PetscOptionsReal("-ex2_pulse_width_time","Width of pulse 'pulse' source","none",rectx->pulse_width,&rectx->pulse_width, NULL);
646:   PetscOptionsReal("-ex2_pulse_rate","Number density of pulse for 'pulse' source","none",rectx->pulse_rate,&rectx->pulse_rate, NULL);
647:   rectx->T_cold *= 1.16e7; /* convert to Kelvin */
648:   PetscOptionsReal("-ex2_ion_potential","Potential to ionize impurity (should be array) in ev","none",rectx->ion_potential,&rectx->ion_potential, NULL);
649:   PetscOptionsReal("-ex2_inductance","","none",rectx->L,&rectx->L, NULL);
650:   PetscInfo5(dummy, "Num electrons from ions=%g, T_cold=%10.3e, ion potential=%10.3e, E_z=%10.3e v_0=%10.3e\n",rectx->Ne_ion,rectx->T_cold,rectx->ion_potential,ctx->Ez,ctx->v_0);
651:   PetscOptionsEnd();
652:   /* get impurity source rate function */
653:   PetscFunctionListFind(plist,pname,&rectx->impuritySrcRate);
654:   if (!rectx->impuritySrcRate) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"No impurity source function found '%s'",pname);
655:   PetscFunctionListFind(testlist,testname,&rectx->test);
656:   if (!rectx->test) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"No test found '%s'",testname);
657:   PetscFunctionListFind(elist,ename,&rectx->E);
658:   if (!rectx->E) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"No E field function found '%s'",ename);
659:   PetscFunctionListDestroy(&plist);
660:   PetscFunctionListDestroy(&testlist);
661:   PetscFunctionListDestroy(&elist);
662:   {
663:     PetscMPIInt    rank;
664:     MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
665:     if (rank) { /* turn off output stuff for duplicate runs */
666:       PetscOptionsClearValue(NULL,"-dm_view");
667:       PetscOptionsClearValue(NULL,"-vec_view");
668:       PetscOptionsClearValue(NULL,"-dm_view_diff");
669:       PetscOptionsClearValue(NULL,"-vec_view_diff");
670:       PetscOptionsClearValue(NULL,"-dm_view_sources");
671:       PetscOptionsClearValue(NULL,"-vec_view_sources");
672:     }
673:   }
674:   /* convert E from Conner-Hastie E_c units to real */
675:   {
676:     PetscReal E = ctx->Ez, Tev = ctx->thermal_temps[0]*8.621738e-5, n = ctx->n_0*ctx->n[0];
677:     CalculateE(Tev, n, ctx->lnLam, ctx->epsilon0, &E);
678:     ((LandauCtx*)ctx)->Ez *= E;
679:   }
680:   DMDestroy(&dummy);
681:   return(0);
682: }

684: int main(int argc, char **argv)
685: {
686:   DM             dm;
687:   Vec            X;
689:   PetscInt       dim = 2;
690:   TS             ts;
691:   Mat            J;
692:   PetscDS        prob;
693:   LandauCtx      *ctx;
694:   REctx          *rectx;

696:   PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr;
697:   PetscOptionsGetInt(NULL,NULL, "-dim", &dim, NULL);
698:   /* Create a mesh */
699:   LandauCreateVelocitySpace(PETSC_COMM_SELF, dim, "", &X, &J, &dm);
700:   LandauCreateMassMatrix(dm, NULL);
701:   DMGetApplicationContext(dm, &ctx);
702:   DMSetUp(dm);
703:   DMGetDS(dm, &prob);
704:   s_quarter3DDomain_notused = ctx->quarter3DDomain; /* ugh */
705:   /* context */
706:   PetscNew(&rectx);
707:   ctx->data = rectx;
708:   ProcessREOptions(rectx,ctx,dm,"");
709:   DMSetOutputSequenceNumber(dm, 0, 0.0);
710:   DMViewFromOptions(dm,NULL,"-dm_view");
711:   DMViewFromOptions(dm,NULL,"-dm_view_sources");
712:   DMViewFromOptions(dm,NULL,"-dm_view_diff");
713:   /* Create timestepping solver context */
714:   TSCreate(PETSC_COMM_SELF,&ts);
715:   TSSetDM(ts,dm);
716:   TSSetIFunction(ts,NULL,LandauIFunction,NULL);
717:   TSSetIJacobian(ts,J,J,LandauIJacobian,NULL);
718:   TSSetRHSFunction(ts,NULL,FormSource,NULL);
719:   TSSetFromOptions(ts);
720:   TSSetSolution(ts,X);
721:   TSSetApplicationContext(ts, ctx);
722:   TSMonitorSet(ts,Monitor,ctx,NULL);
723:   TSSetPreStep(ts,PreStep);
724:   rectx->Ez_initial = ctx->Ez;       /* cache for induction caclulation - applied E field */
725:   if (1) { /* warm up an test just LandauIJacobian */
726:     #if defined PETSC_USE_LOG
727:     PetscLogStage stage;
728:     #endif
729:     Vec           vec;
730:     VecDuplicate(X,&vec);
731:     VecSet(vec,1.);
732:     /* warm up */
733:     PetscLogStageRegister("Warmup", &stage);
734:     PetscLogStagePush(stage);
735:     LandauIJacobian(ts, 0.0, vec, vec, 0.0, J, J, ctx);
736:     PetscObjectSetName((PetscObject)J, "Jacobian");
737:     MatViewFromOptions(J,NULL,"-initial_mat_view");
738:     PetscObjectSetName((PetscObject)J, "Mass + Jacobian");
739:     //VecScale(vec,10.*PETSC_MACHINE_EPSILON);
740:     LandauIJacobian(ts, 0.0, vec, vec, 1.0, J, J, ctx);
741:     MatViewFromOptions(J,NULL,"-initial_mat_view");
742:     PetscObjectSetName((PetscObject)J, "Jacobian");
743:     PetscLogStagePop();
744:     VecDestroy(&vec);
745:   }
746:   VecViewFromOptions(X,NULL,"-vec_view"); // inital condition (monitor plots after step)
747:   /* go */
748:   TSSolve(ts,X);
749:   /* clean up */
750:   LandauDestroyVelocitySpace(&dm);
751:   TSDestroy(&ts);
752:   VecDestroy(&X);
753:   if (rectx->imp_src) {
754:     VecDestroy(&rectx->imp_src);
755:   }
756:   PetscFree(rectx);
757:   PetscFinalize();
758:   return ierr;
759: }

761: /*TEST
762:   test:
763:     suffix: 0
764:     requires: p4est !complex double
765:     args: -dm_landau_Ez 0 -petscspace_degree 2 -petscspace_poly_tensor 1 -dm_landau_type p4est -info :dm,tsadapt -dm_landau_ion_masses 2 -dm_landau_ion_charges 1 -dm_landau_thermal_temps 5,5 -dm_landau_n 2,2 -dm_landau_n_0 5e19 -ts_monitor -snes_rtol 1.e-10 -snes_stol 1.e-14 -snes_monitor -snes_converged_reason -snes_max_it 10 -ts_type arkimex -ts_arkimex_type 1bee -ts_max_snes_failures -1 -ts_rtol 1e-3 -ts_dt 1.e-1 -ts_max_time 1 -ts_adapt_clip .5,1.25 -ts_max_steps 2 -ts_adapt_scale_solve_failed 0.75 -ts_adapt_time_step_increase_delay 5 -pc_type lu -ksp_type preonly -dm_landau_amr_levels_max 9 -dm_landau_domain_radius -.75 -ex2_impurity_source_type pulse -ex2_pulse_start_time 1e-1 -ex2_pulse_width_time 10 -ex2_pulse_rate 1e-2 -ex2_t_cold .05 -ex2_plot_dt 1e-1 -dm_refine 1

767: TEST*/