Actual source code: ex2.c

petsc-main 2021-04-20
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 */

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

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

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

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

100: /* < v, n_e (v-shift) > */
101: static void f0_ve_shift(PetscInt dim, PetscInt Nf, PetscInt NfAux,
102:                          const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
103:                          const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
104:                          PetscReal t, const PetscReal x[],  PetscInt numConstants, const PetscScalar constants[], PetscScalar *f0)
105: {
106:   PetscReal vz = numConstants==1 ? PetscRealPart(constants[0]) : 0;
107:   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 */
108:   else {
109:     *f0 =           u[0] *                PetscSqrtReal(x[0]*x[0] + x[1]*x[1] + (x[2]-vz)*(x[2]-vz)); /* n v */
110:   }
111: }

113:  /* CalculateE - Calculate the electric field  */
114:  /*  T        -- Electron temperature  */
115:  /*  n        -- Electron density  */
116:  /*  lnLambda --   */
117:  /*  eps0     --  */
118:  /*  E        -- output E, input \hat E */
119: static PetscReal CalculateE(PetscReal Tev, PetscReal n, PetscReal lnLambda, PetscReal eps0, PetscReal *E)
120: {
121:   PetscReal c,e,m;

124:   c = 299792458;
125:   e = 1.602176e-19;
126:   m = 9.10938e-31;
127:   if (1) {
128:     double Ec, Ehat = *E, betath = PetscSqrtReal(2*Tev*e/(m*c*c)), j0 = Ehat * 7/(PetscSqrtReal(2)*2) * PetscPowReal(betath,3) * n * e * c;
129:     Ec = n*lnLambda*PetscPowReal(e,3) / (4*PETSC_PI*PetscPowReal(eps0,2)*m*c*c);
130:     *E = Ec;
131:     PetscPrintf(PETSC_COMM_WORLD, "CalculateE j0=%g Ec = %g\n",j0,Ec);
132:   } else {
133:     PetscReal Ed, vth;
134:     vth = PetscSqrtReal(8*Tev*e/(m*PETSC_PI));
135:     Ed =  n*lnLambda*PetscPowReal(e,3) / (4*PETSC_PI*PetscPowReal(eps0,2)*m*vth*vth);
136:     *E = Ed;
137:   }
138:   return(0);
139: }

141: static PetscReal Spitzer(PetscReal m_e, PetscReal e, PetscReal Z, PetscReal epsilon0,  PetscReal lnLam, PetscReal kTe_joules)
142: {
143:   PetscReal Fz = (1+1.198*Z+0.222*Z*Z)/(1+2.966*Z+0.753*Z*Z), eta;
144:   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);
145:   return eta;
146: }

148: /*  */
149: static PetscErrorCode testNone(TS ts, Vec X, DM plex, PetscInt stepi, PetscReal time, PetscBool islast, LandauCtx *ctx, REctx *rectx)
150: {
152:   return(0);
153: }

155: /*  */
156: static PetscErrorCode testSpitzer(TS ts, Vec X, DM plex, PetscInt stepi, PetscReal time, PetscBool islast, LandauCtx *ctx, REctx *rectx)
157: {
158:   PetscErrorCode    ierr;
159:   PetscInt          ii;
160:   PetscDS           prob;
161:   static PetscReal  old_ratio = 1e10;
162:   TSConvergedReason reason;
163:   PetscReal         J,J_re,spit_eta,Te_kev=0,E,ratio,Z,n_e,v,v2;
164:   PetscScalar       user[2] = {0.,ctx->charges[0]}, charges[LANDAU_MAX_SPECIES],tt[LANDAU_MAX_SPECIES],vz;
165:   PetscReal         dt;

168:   if (ctx->num_species<2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "ctx->num_species %D < 2",ctx->num_species);
169:   DMGetDS(plex, &prob);
170:   TSGetTimeStep(ts,&dt);
171:   /* get current */
172:   for (ii=0;ii<ctx->num_species;ii++) charges[ii] = ctx->charges[ii];
173:   PetscDSSetConstants(prob, ctx->num_species, charges);
174:   PetscDSSetObjective(prob, 0, &f0_jz_sum);
175:   DMPlexComputeIntegralFEM(plex,X,tt,NULL);
176:   J = -ctx->n_0*ctx->v_0*PetscRealPart(tt[0]);
177:   /* get N_e */
178:   PetscDSSetConstants(prob, 2, user);
179:   PetscDSSetObjective(prob, 0, &f0_n);
180:   DMPlexComputeIntegralFEM(plex,X,tt,NULL);
181:   n_e = PetscRealPart(tt[0])*ctx->n_0;
182:   /* Z */
183:   Z = -ctx->charges[1]/ctx->charges[0];
184:   if (ctx->charges[rectx->imp_idx] != ctx->charges[1]) {
185:     PetscReal   Znew, n_i1,n_ix;
186:     user[0] = 1.0;
187:     PetscDSSetConstants(prob, 2, user);
188:     DMPlexComputeIntegralFEM(plex,X,tt,NULL);
189:     n_i1 = PetscRealPart(tt[0])*ctx->n_0;
190:     user[0] = (PetscScalar)rectx->imp_idx;
191:     PetscDSSetConstants(prob, 2, user);
192:     DMPlexComputeIntegralFEM(plex,X,tt,NULL);
193:     n_ix = PetscRealPart(tt[0])*ctx->n_0;
194:     Znew = -(ctx->charges[1]*n_i1 + ctx->charges[rectx->imp_idx]*n_ix)/(ctx->charges[0]*(n_i1 + n_ix));
195:     Z = Znew;
196:   }
197:   /* remove drift */
198:   if (1) {
199:     user[0] = .0;
200:     PetscDSSetConstants(prob, 2, user);
201:     PetscDSSetObjective(prob, 0, &f0_vz);
202:     DMPlexComputeIntegralFEM(plex,X,tt,NULL);
203:     vz = ctx->n_0*PetscRealPart(tt[0])/n_e; /* non-dimensional */
204:   } else vz = 0;
205:   /* thermal velocity */
206:   PetscDSSetConstants(prob, 1, &vz);
207:   PetscDSSetObjective(prob, 0, &f0_ve_shift);
208:   DMPlexComputeIntegralFEM(plex,X,tt,NULL);
209:   v = ctx->n_0*ctx->v_0*PetscRealPart(tt[0])/n_e;   /* remove number density to get velocity */
210:   v2 = PetscSqr(v);                                 /* use real space: m^2 / s^2 */
211:   Te_kev = (v2*ctx->masses[0]*PETSC_PI/8)*kev_joul; /* temperature in kev */
212:   spit_eta = Spitzer(ctx->masses[0],-ctx->charges[0],Z,ctx->epsilon0,ctx->lnLam,Te_kev/kev_joul); /* kev --> J (kT) */
213:   if (1) {
214:     PetscDSSetConstants(prob, 1, charges);
215:     PetscDSSetObjective(prob, 0, &f0_j_re);
216:     DMPlexComputeIntegralFEM(plex,X,tt,NULL);
217:   } else tt[0] = 0;
218:   J_re = -ctx->n_0*ctx->v_0*PetscRealPart(tt[0]);

220:   if (rectx->use_spitzer_eta) {
221:     E = ctx->Ez = spit_eta*(rectx->j-J_re);
222:   } else {
223:     E = ctx->Ez; /* keep real E */
224:     rectx->j = J; /* cache */
225:   }

227:   ratio = E/J/spit_eta;
228:   if (stepi>10 && !rectx->use_spitzer_eta && ((old_ratio-ratio < 1.e-3 && ratio > 0.99 && ratio < 1.01) || (old_ratio-ratio < 1.e-4 && ratio > 0.98 && ratio < 1.02))) {
229:     rectx->pulse_start = time + 0.98*dt;
230:     rectx->use_spitzer_eta = PETSC_TRUE;
231:   }
232:   TSGetConvergedReason(ts,&reason);
233:   TSGetConvergedReason(ts,&reason);
234:   if ((rectx->plotting) || stepi == 0 || reason || rectx->pulse_start == time + 0.98*dt) {
235:     PetscPrintf(ctx->comm, "testSpitzer: %4D) time=%11.4e n_e= %10.3e E= %10.3e J= %10.3e J_re= %10.3e %.3g %% Te_kev= %10.3e Z_eff=%g E/J to eta ratio=%g (diff=%g) %s %s\n",stepi,time,n_e/ctx->n_0,ctx->Ez,J,J_re,100*J_re/J, Te_kev,Z,ratio,old_ratio-ratio, rectx->use_spitzer_eta ? "using Spitzer eta*J E" : "constant E",rectx->pulse_start != time + 0.98*dt ? "normal" : "transition");
236:   }
237:   old_ratio = ratio;
238:   return(0);
239: }

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

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

285:   VecGetDM(X, &dm);
286:   DMGetDS(plex, &prob);
287:   VecDuplicate(X,&X2);
288:   VecCopy(X,X2);
289:   if (!rectx->X_0) {
290:     VecDuplicate(X,&rectx->X_0);
291:     VecCopy(X,rectx->X_0);
292:   }
293:   VecAXPY(X,-1.0,rectx->X_0);
294:   PetscDSSetConstants(prob, sizeof(LandauCtx)/sizeof(PetscScalar), (PetscScalar*)ctx);
295:   rectx->idx = 0;
296:   PetscDSSetObjective(prob, 0, &f0_0_diff_lp);
297:   DMPlexComputeIntegralFEM(plex,X2,tt,NULL);
298:   ediff = PetscPowReal(PetscRealPart(tt[0]),1./ppp);
299:   PetscDSSetObjective(prob, 0, &f0_0_maxwellian_lp);
300:   DMPlexComputeIntegralFEM(plex,X2,tt,NULL);
301:   lpm0 = PetscPowReal(PetscRealPart(tt[0]),1./ppp);
302:   if (ctx->num_species>1) {
303:     rectx->idx = 1;
304:     PetscDSSetObjective(prob, 0, &f0_0_diff_lp);
305:     DMPlexComputeIntegralFEM(plex,X2,tt,NULL);
306:     idiff = PetscPowReal(PetscRealPart(tt[0]),1./ppp);
307:     PetscDSSetObjective(prob, 0, &f0_0_maxwellian_lp);
308:     DMPlexComputeIntegralFEM(plex,X2,tt,NULL);
309:     lpm1 = PetscPowReal(PetscRealPart(tt[0]),1./ppp);
310:   }
311:   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);
312:   /* view */
313:   VecViewFromOptions(X,NULL,"-vec_view_diff");
314:   VecCopy(X2,X);
315:   VecDestroy(&X2);
316:   if (islast) {
317:     VecDestroy(&rectx->X_0);
318:     rectx->X_0 = NULL;
319:   }
320:   return(0);
321: }

323: static PetscErrorCode EInduction(Vec X, Vec X_t, PetscInt step, PetscReal time, LandauCtx *ctx, PetscReal *a_E)
324: {
325:   REctx             *rectx = (REctx*)ctx->data;
326:   PetscErrorCode    ierr;
327:   PetscInt          ii;
328:   DM                dm,plex;
329:   PetscScalar       tt[LANDAU_MAX_SPECIES], constants[LANDAU_MAX_SPECIES];
330:   PetscReal         dJ_dt;
331:   PetscDS           prob;

334:   for (ii=0;ii<ctx->num_species;ii++) constants[ii] = ctx->charges[ii];
335:   VecGetDM(X, &dm);
336:   DMGetDS(dm, &prob);
337:   DMConvert(dm, DMPLEX, &plex);
338:   /* get d current / dt */
339:   PetscDSSetConstants(prob, ctx->num_species, constants);
340:   PetscDSSetObjective(prob, 0, &f0_jz_sum);
341:   if (!X_t) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "X_t");
342:   DMPlexComputeIntegralFEM(plex,X_t,tt,NULL);
343:   dJ_dt = -ctx->n_0*ctx->v_0*PetscRealPart(tt[0])/ctx->t_0;
344:   /* E induction */
345:   *a_E = -rectx->L*dJ_dt + rectx->Ez_initial;
346:   DMDestroy(&plex);
347:   return(0);
348: }

350: static PetscErrorCode EConstant(Vec X,  Vec X_t, PetscInt step, PetscReal time, LandauCtx *ctx, PetscReal *a_E)
351: {
353:   *a_E = ctx->Ez;
354:   return(0);
355: }

357: static PetscErrorCode ENone(Vec X,  Vec X_t, PetscInt step, PetscReal time, LandauCtx *ctx, PetscReal *a_E)
358: {
360:   *a_E = 0;
361:   return(0);
362: }

364: /* ------------------------------------------------------------------- */
365: /*
366:    FormSource - Evaluates source terms F(t).

368:    Input Parameters:
369: .  ts - the TS context
370: .  time -
371: .  X_dummmy - input vector
372: .  dummy - optional user-defined context, as set by SNESSetFunction()

374:    Output Parameter:
375: .  F - function vector
376:  */
377: static PetscErrorCode FormSource(TS ts,PetscReal ftime,Vec X_dummmy, Vec F,void *dummy)
378: {
379:   PetscReal      new_imp_rate;
380:   LandauCtx      *ctx;
381:   DM             dm,plex;
383:   REctx          *rectx;

386:   TSGetDM(ts,&dm);
387:   DMGetApplicationContext(dm, &ctx);
388:   rectx = (REctx*)ctx->data;
389:   /* check for impurities */
390:   rectx->impuritySrcRate(ftime,&new_imp_rate,ctx);
391:   if (new_imp_rate != 0) {
392:     if (new_imp_rate != rectx->current_rate) {
393:       PetscInt       ii;
394:       PetscReal      dne_dt,dni_dt,tilda_ns[LANDAU_MAX_SPECIES],temps[LANDAU_MAX_SPECIES];
395:       PetscDS        prob; /* diagnostics only */
396:       rectx->current_rate = new_imp_rate;
397:       DMConvert(dm, DMPLEX, &plex);
398:       DMGetDS(dm, &prob);
399:       dni_dt = new_imp_rate              /* *ctx->t_0 */; /* fully ionized immediately, no normalize, stay in non-dim */
400:       dne_dt = new_imp_rate*rectx->Ne_ion/* *ctx->t_0 */;
401:       PetscInfo4(dm, "\tHave 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);
402:       for (ii=1;ii<LANDAU_MAX_SPECIES;ii++) tilda_ns[ii] = 0;
403:       for (ii=1;ii<LANDAU_MAX_SPECIES;ii++)    temps[ii] = 1;
404:       tilda_ns[0] = dne_dt;        tilda_ns[rectx->imp_idx] = dni_dt;
405:       temps[0]    = rectx->T_cold;    temps[rectx->imp_idx] = rectx->T_cold;
406:       /* add it */
407:       if (!rectx->imp_src) {
408:         DMCreateGlobalVector(dm, &rectx->imp_src);
409:         PetscObjectSetName((PetscObject)rectx->imp_src, "source");
410:       }
411:       VecZeroEntries(rectx->imp_src);
412:       LandauAddMaxwellians(plex,rectx->imp_src,ftime,temps,tilda_ns,ctx);
413:       /* clean up */
414:       DMDestroy(&plex);
415:       VecViewFromOptions(rectx->imp_src,NULL,"-vec_view_sources");
416:     }
417:     VecCopy(rectx->imp_src,F);
418:   } else {
419:     if (rectx->current_rate != 0 && rectx->imp_src) {
420:       VecZeroEntries(rectx->imp_src);
421:     }
422:     VecZeroEntries(F);
423:     rectx->current_rate = 0;
424:   }
425:   return(0);
426: }
427: PetscErrorCode Monitor(TS ts, PetscInt stepi, PetscReal time, Vec X, void *actx)
428: {
429:   LandauCtx         *ctx = (LandauCtx*) actx;   /* user-defined application context */
430:   REctx             *rectx = (REctx*)ctx->data;
431:   DM                dm,plex;
432:   TSConvergedReason reason;
433:   PetscErrorCode    ierr;
435:   VecGetDM(X, &dm);
436:   if (stepi > rectx->plotStep && rectx->plotting) {
437:     rectx->plotting = PETSC_FALSE; /* was doing diagnostics, now done */
438:     rectx->plotIdx++;
439:   }
440:   /* view */
441:   TSGetConvergedReason(ts,&reason);
442:   if (time/rectx->plotDt >= (PetscReal)rectx->plotIdx || reason) {
443:     if ((reason || stepi==0 || rectx->plotIdx%10==0) && ctx->verbose > 0){
444:       /* print norms */
445:       LandauPrintNorms(X, stepi);
446:     }
447:     if (!rectx->plotting) { /* first step of possible backtracks */
448:       rectx->plotting = PETSC_TRUE;
449:       DMConvert(dm, DMPLEX, &plex);
450:       /* diagnostics + change E field with Sptizer (not just a monitor) */
451:       rectx->test(ts,X,plex,stepi,time,reason ? PETSC_TRUE : PETSC_FALSE, ctx, rectx);
452:       DMDestroy(&plex);
453:     } else {
454:       PetscPrintf(PETSC_COMM_WORLD, "\t\t ERROR SKIP test spit ------\n");
455:       rectx->plotting = PETSC_TRUE;
456:     }
457:     /* view, overwrite step when back tracked */
458:     DMSetOutputSequenceNumber(dm, rectx->plotIdx, time*ctx->t_0);
459:     VecViewFromOptions(X,NULL,"-vec_view");
460:     rectx->plotStep = stepi;
461:   } else {
462:     if (rectx->plotting) PetscPrintf(PETSC_COMM_WORLD," ERROR rectx->plotting=%D step %D\n",rectx->plotting,stepi);
463:     DMConvert(dm, DMPLEX, &plex);
464:     /* diagnostics + change E field with Sptizer (not just a monitor) - can we lag this? */
465:     rectx->test(ts,X,plex,stepi,time,reason ? PETSC_TRUE : PETSC_FALSE, ctx, rectx);
466:     DMDestroy(&plex);
467:   }
468:   /* parallel check */
469:   if (reason && ctx->verbose > 0) {
470:     PetscReal    val,rval;
471:     PetscMPIInt  rank;
472:     MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
473:     TSGetSolution(ts, &X);
474:     VecNorm(X,NORM_2,&val);
475:     MPIU_Allreduce(&val,&rval,1,MPIU_REAL,MPIU_MAX,PETSC_COMM_WORLD);
476:     if (rval != val) {
477:       PetscPrintf(PETSC_COMM_SELF, " ***** [%D] ERROR max |x| = %22.15e, my |x| = %22.15e diff=%e\n",rank,rval,val,rval-val);
478:     } else {
479:       PetscPrintf(PETSC_COMM_WORLD, "[%D] parallel consistency check OK\n",rank);
480:     }
481:   }
482:   rectx->idx = 0;
483:   return(0);
484: }

486: PetscErrorCode PreStep(TS ts)
487: {
489:   LandauCtx      *ctx;
490:   REctx          *rectx;
491:   DM             dm;
492:   PetscInt       stepi;
493:   PetscReal      time;
494:   Vec            X;

497:   /* not used */
498:   TSGetDM(ts,&dm);
499:   TSGetTime(ts,&time);
500:   TSGetSolution(ts,&X);
501:   DMGetApplicationContext(dm, &ctx);
502:   rectx = (REctx*)ctx->data;
503:   TSGetStepNumber(ts, &stepi);
504:   /* update E */
505:   rectx->E(X, NULL, stepi, time, ctx, &ctx->Ez);
506:   return(0);
507: }

509: /* 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) */
510: static PetscErrorCode stepSrc(PetscReal time, PetscReal *rho, LandauCtx *ctx)
511: {
512:   REctx         *rectx = (REctx*)ctx->data;

515:   if (time >= rectx->pulse_start) *rho = rectx->pulse_rate;
516:   else *rho = 0.;
517:   return(0);
518: }
519: static PetscErrorCode zeroSrc(PetscReal time, PetscReal *rho, LandauCtx *ctx)
520: {
522:   *rho = 0.;
523:   return(0);
524: }
525: static PetscErrorCode pulseSrc(PetscReal time, PetscReal *rho, LandauCtx *ctx)
526: {
527:   REctx *rectx = (REctx*)ctx->data;

530:   if (rectx->pulse_start == PETSC_MAX_REAL) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"'-ex2_pulse_start_time X' must be used with '-ex2_impurity_source_type pulse'");
531:   if (time < rectx->pulse_start || time > rectx->pulse_start + 3*rectx->pulse_width) *rho = 0;
532:   /* else if (0) { */
533:   /*   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); */
534:   /*   *rho = rectx->pulse_rate * (cycle / (stop - start)) / (1 + PetscExpReal(steep*(PetscSinReal(2*PETSC_PI*((t - start)/cycle + xi)) - PetscSinReal(2*PETSC_PI*xi)))); */
535:   /* } else if (0) { */
536:   /*   double x = 2*(time - rectx->pulse_start)/(3*rectx->pulse_width) - 1; */
537:   /*   if (x==1 || x==-1) *rho = 0; */
538:   /*   else *rho = rectx->pulse_rate * PetscExpReal(-1/(1-x*x)); */
539:   /* } */
540:   else {
541:     double x = PetscSinReal((time-rectx->pulse_start)/(3*rectx->pulse_width)*2*PETSC_PI - PETSC_PI/2) + 1; /* 0:2, integrates to 1.0 */
542:     *rho = rectx->pulse_rate * x / (3*rectx->pulse_width);
543:     if (!rectx->use_spitzer_eta) rectx->use_spitzer_eta = PETSC_TRUE; /* use it next time */
544:   }
545:   return(0);
546: }

550: static PetscErrorCode ProcessREOptions(REctx *rectx, const LandauCtx *ctx, DM dm, const char prefix[])
551: {
552:   PetscErrorCode    ierr;
553:   PetscFunctionList plist = NULL, testlist = NULL, elist = NULL;
554:   char              pname[256],testname[256],ename[256];
555:   DM                dm_dummy;
556:   PetscBool         Connor_E = PETSC_FALSE;

559:   DMCreate(PETSC_COMM_WORLD,&dm_dummy);
560:   rectx->Ne_ion = 1;                 /* number of electrons given up by impurity ion */
561:   rectx->T_cold = .005;              /* kev */
562:   rectx->ion_potential = 15;         /* ev */
563:   rectx->L = 2;
564:   rectx->X_0 = NULL;
565:   rectx->imp_idx = ctx->num_species - 1; /* default ionized impurity as last one */
566:   rectx->pulse_start = PETSC_MAX_REAL;
567:   rectx->pulse_width = 1;
568:   rectx->plotStep = PETSC_MAX_INT;
569:   rectx->pulse_rate = 1.e-1;
570:   rectx->current_rate = 0;
571:   rectx->plotIdx = 0;
572:   rectx->imp_src = 0;
573:   rectx->j = 0;
574:   rectx->plotDt = 1.0;
575:   rectx->plotting = PETSC_FALSE;
576:   rectx->use_spitzer_eta = PETSC_FALSE;
577:   rectx->idx = 0;
578:   /* Register the available impurity sources */
579:   PetscFunctionListAdd(&plist,"step",&stepSrc);
580:   PetscFunctionListAdd(&plist,"none",&zeroSrc);
581:   PetscFunctionListAdd(&plist,"pulse",&pulseSrc);
582:   PetscStrcpy(pname,"none");
583:   PetscFunctionListAdd(&testlist,"none",&testNone);
584:   PetscFunctionListAdd(&testlist,"spitzer",&testSpitzer);
585:   PetscFunctionListAdd(&testlist,"stable",&testStable);
586:   PetscStrcpy(testname,"none");
587:   PetscFunctionListAdd(&elist,"none",&ENone);
588:   PetscFunctionListAdd(&elist,"induction",&EInduction);
589:   PetscFunctionListAdd(&elist,"constant",&EConstant);
590:   PetscStrcpy(ename,"constant");

592:   PetscOptionsBegin(PETSC_COMM_SELF, prefix, "Options for Runaway/seed electron model", "none");
593:   PetscOptionsReal("-ex2_plot_dt", "Plotting interval", "xgc_dmplex.c", rectx->plotDt, &rectx->plotDt, NULL);
594:   if (rectx->plotDt < 0) rectx->plotDt = 1e30;
595:   if (rectx->plotDt == 0) rectx->plotDt = 1e-30;
596:   PetscOptionsFList("-ex2_impurity_source_type","Name of impurity source to run","",plist,pname,pname,sizeof(pname),NULL);
597:   PetscOptionsFList("-ex2_test_type","Name of test to run","",testlist,testname,testname,sizeof(testname),NULL);
598:   PetscOptionsInt("-ex2_impurity_index", "index of sink for impurities", "none", rectx->imp_idx, &rectx->imp_idx, NULL);
599:   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);
600:   PetscOptionsFList("-ex2_e_field_type","Electric field type","",elist,ename,ename,sizeof(ename),NULL);
601:   rectx->Ne_ion = -ctx->charges[rectx->imp_idx]/ctx->charges[0];
602:   PetscOptionsReal("-ex2_t_cold","Temperature of cold electron and ions after ionization in keV","none",rectx->T_cold,&rectx->T_cold, NULL);
603:   PetscOptionsReal("-ex2_pulse_start_time","Time at which pulse happens for 'pulse' source","none",rectx->pulse_start,&rectx->pulse_start, NULL);
604:   PetscOptionsReal("-ex2_pulse_width_time","Width of pulse 'pulse' source","none",rectx->pulse_width,&rectx->pulse_width, NULL);
605:   PetscOptionsReal("-ex2_pulse_rate","Number density of pulse for 'pulse' source","none",rectx->pulse_rate,&rectx->pulse_rate, NULL);
606:   rectx->T_cold *= 1.16e7; /* convert to Kelvin */
607:   PetscOptionsReal("-ex2_ion_potential","Potential to ionize impurity (should be array) in ev","none",rectx->ion_potential,&rectx->ion_potential, NULL);
608:   PetscOptionsReal("-ex2_inductance","Inductance E feild","none",rectx->L,&rectx->L, NULL);
609:   PetscOptionsBool("-ex2_connor_e_field_units","Scale Ex but Connor-Hastie E_c","none",Connor_E,&Connor_E, NULL);
610:   PetscInfo5(dm_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);
611:   PetscOptionsEnd();
612:   /* get impurity source rate function */
613:   PetscFunctionListFind(plist,pname,&rectx->impuritySrcRate);
614:   if (!rectx->impuritySrcRate) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"No impurity source function found '%s'",pname);
615:   PetscFunctionListFind(testlist,testname,&rectx->test);
616:   if (!rectx->test) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"No test found '%s'",testname);
617:   PetscFunctionListFind(elist,ename,&rectx->E);
618:   if (!rectx->E) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"No E field function found '%s'",ename);
619:   PetscFunctionListDestroy(&plist);
620:   PetscFunctionListDestroy(&testlist);
621:   PetscFunctionListDestroy(&elist);

623:   /* convert E from Connor-Hastie E_c units to real if doing Spitzer E */
624:   if (Connor_E) {
625:     PetscReal E = ctx->Ez, Tev = ctx->thermal_temps[0]*8.621738e-5, n = ctx->n_0*ctx->n[0];
626:     CalculateE(Tev, n, ctx->lnLam, ctx->epsilon0, &E);
627:     ((LandauCtx*)ctx)->Ez *= E;
628:   }
629:   DMDestroy(&dm_dummy);
630:   return(0);
631: }

633: int main(int argc, char **argv)
634: {
635:   DM             dm;
636:   Vec            X;
638:   PetscInt       dim = 2;
639:   TS             ts;
640:   Mat            J;
641:   PetscDS        prob;
642:   LandauCtx      *ctx;
643:   REctx          *rectx;
644: #if defined PETSC_USE_LOG
645:   PetscLogStage stage;
646: #endif
647:   PetscMPIInt    rank;
648:   PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr;
649:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
650:   if (rank) { /* turn off output stuff for duplicate runs */
651:     PetscOptionsClearValue(NULL,"-dm_view");
652:     PetscOptionsClearValue(NULL,"-vec_view");
653:     PetscOptionsClearValue(NULL,"-dm_view_diff");
654:     PetscOptionsClearValue(NULL,"-vec_view_diff");
655:     PetscOptionsClearValue(NULL,"-dm_view_sources");
656:     PetscOptionsClearValue(NULL,"-vec_view_sources");
657:     PetscOptionsClearValue(NULL,"-info"); /* this does not work */
658:   }
659:   PetscOptionsGetInt(NULL,NULL, "-dim", &dim, NULL);
660:   /* Create a mesh */
661:   LandauCreateVelocitySpace(PETSC_COMM_WORLD, dim, "", &X, &J, &dm);
662:   PetscObjectSetName((PetscObject)J, "Jacobian");
663:   LandauCreateMassMatrix(dm, NULL);
664:   DMGetApplicationContext(dm, &ctx);
665:   DMSetUp(dm);
666:   DMGetDS(dm, &prob);
667:   /* context */
668:   PetscNew(&rectx);
669:   ctx->data = rectx;
670:   ProcessREOptions(rectx,ctx,dm,"");
671:   DMSetOutputSequenceNumber(dm, 0, 0.0);
672:   DMViewFromOptions(dm,NULL,"-dm_view");
673:   DMViewFromOptions(dm,NULL,"-dm_view_sources");
674:   DMViewFromOptions(dm,NULL,"-dm_view_diff");
675:   /* Create timestepping solver context */
676:   TSCreate(PETSC_COMM_SELF,&ts);
677:   TSSetDM(ts,dm);
678:   TSSetIFunction(ts,NULL,LandauIFunction,NULL);
679:   TSSetIJacobian(ts,J,J,LandauIJacobian,NULL);
680:   TSSetRHSFunction(ts,NULL,FormSource,NULL);
681:   TSSetFromOptions(ts);
682:   TSSetSolution(ts,X);
683:   TSSetApplicationContext(ts, ctx);
684:   TSMonitorSet(ts,Monitor,ctx,NULL);
685:   TSSetPreStep(ts,PreStep);
686:   rectx->Ez_initial = ctx->Ez;       /* cache for induction caclulation - applied E field */
687:   if (1) { /* warm up an test just LandauIJacobian */
688:     Vec           vec;
689:     PetscInt      nsteps;
690:     PetscReal     dt;
691:     PetscLogStageRegister("Warmup", &stage);
692:     PetscLogStagePush(stage);
693:     VecDuplicate(X,&vec);
694:     VecCopy(X,vec);
695:     TSGetMaxSteps(ts,&nsteps);
696:     TSGetTimeStep(ts,&dt);
697:     TSSetMaxSteps(ts,1);
698:     TSSetSolution(ts,vec);
699:     TSSolve(ts,vec);
700:     TSSetMaxSteps(ts,nsteps);
701:     TSSetStepNumber(ts,0);
702:     TSSetTime(ts,0);
703:     TSSetTimeStep(ts,dt);
704:     TSSetSolution(ts,X);
705:     rectx->plotIdx = 0;
706:     rectx->plotting = PETSC_FALSE;
707:     PetscLogStagePop();
708:     VecDestroy(&vec);
709:     ctx->aux_bool = PETSC_FALSE; // flag for not a clean Jacobian
710:   }
711:   VecViewFromOptions(X,NULL,"-vec_view"); // inital condition (monitor plots after step)
712:   /* go */
713:   PetscLogStageRegister("Solve", &stage);
714:   MPI_Barrier(MPI_COMM_WORLD);
715:   PetscLogStagePush(stage);
716:   TSSolve(ts,X);
717:   PetscLogStagePop();
718:   /* clean up */
719:   LandauDestroyVelocitySpace(&dm);
720:   TSDestroy(&ts);
721:   VecDestroy(&X);
722:   if (rectx->imp_src) {
723:     VecDestroy(&rectx->imp_src);
724:   }
725:   PetscFree(rectx);
726:   PetscFinalize();
727:   return ierr;
728: }

730: /*TEST
731:   test:
732:     suffix: 0
733:     requires: p4est !complex double
734:     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 -dm_preallocate_only -dm_landau_gpu_assembly true -dm_landau_device_type cpu

736:   test:
737:     suffix: kokkos
738:     requires: p4est !complex double kokkos_kernels !define(PETSC_USE_CTABLE)
739:     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 -dm_preallocate_only -dm_landau_device_type kokkos -dm_landau_gpu_assembly true -dm_mat_type aijkokkos -dm_vec_type kokkos

741:   test:
742:     suffix: cuda
743:     requires: p4est !complex double cuda !define(PETSC_USE_CTABLE)
744:     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 -dm_preallocate_only -dm_landau_device_type cuda -dm_landau_gpu_assembly true -dm_mat_type aijcusparse -dm_vec_type cuda

746: TEST*/