1 static char help[] = "Runaway electron model with Landau collision operator\n\n"; 2 3 #include <petscdmplex.h> 4 #include <petsclandau.h> 5 #include <petscts.h> 6 #include <petscds.h> 7 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; 33 34 static const PetscReal kev_joul = 6.241506479963235e+15; /* 1/1000e */ 35 static PetscBool s_quarter3DDomain_notused = PETSC_FALSE; 36 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 } 60 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 } 76 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 } 90 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 } 104 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 } 118 119 static PetscErrorCode getTe_kev(DM plex, Vec X, PetscReal *a_n, PetscReal *a_Tkev) 120 { 121 PetscErrorCode ierr; 122 LandauCtx *ctx; 123 124 PetscFunctionBeginUser; 125 ierr = DMGetApplicationContext(plex, &ctx);CHKERRQ(ierr); 126 { 127 PetscDS prob; 128 PetscReal v2, v, n; 129 PetscScalar tt[LANDAU_MAX_SPECIES],user[2] = {0.,ctx->charges[0]}, vz; 130 ierr = DMGetDS(plex, &prob);CHKERRQ(ierr); 131 ierr = PetscDSSetConstants(prob, 2, user);CHKERRQ(ierr); 132 ierr = PetscDSSetObjective(prob, 0, &f0_n);CHKERRQ(ierr); 133 ierr = DMPlexComputeIntegralFEM(plex,X,tt,NULL);CHKERRQ(ierr); 134 n = ctx->n_0*PetscRealPart(tt[0]); 135 ierr = PetscDSSetObjective(prob, 0, &f0_vz);CHKERRQ(ierr); 136 ierr = DMPlexComputeIntegralFEM(plex,X,tt,NULL);CHKERRQ(ierr); 137 vz = ctx->n_0*PetscRealPart(tt[0])/n; /* non-dimensional */ 138 /* remove drift */ 139 ierr = PetscDSSetConstants(prob, 1, &vz);CHKERRQ(ierr); 140 ierr = PetscDSSetObjective(prob, 0, &f0_ve_shift);CHKERRQ(ierr); 141 ierr = DMPlexComputeIntegralFEM(plex,X,tt,NULL);CHKERRQ(ierr); 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 PetscFunctionReturn(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; 158 159 PetscFunctionBegin; 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 PetscFunctionReturn(0); 175 } 176 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 } 183 184 /* */ 185 static PetscErrorCode testNone(TS ts, Vec X, DM plex, PetscInt stepi, PetscReal time, PetscBool islast, LandauCtx *ctx, REctx *rectx) 186 { 187 PetscFunctionBeginUser; 188 PetscFunctionReturn(0); 189 } 190 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]; 201 202 PetscFunctionBeginUser; 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 ierr = DMGetDS(plex, &prob);CHKERRQ(ierr); 207 ierr = PetscDSSetConstants(prob, 2, user);CHKERRQ(ierr); 208 ierr = PetscDSSetObjective(prob, 0, &f0_n);CHKERRQ(ierr); 209 ierr = DMPlexComputeIntegralFEM(plex,X,tt,NULL);CHKERRQ(ierr); 210 n_e = PetscRealPart(tt[0])*ctx->n_0; 211 ierr = PetscDSSetConstants(prob, ctx->num_species, constants);CHKERRQ(ierr); 212 ierr = PetscDSSetObjective(prob, 0, &f0_jz_sum);CHKERRQ(ierr); 213 ierr = DMPlexComputeIntegralFEM(plex,X,tt,NULL);CHKERRQ(ierr); 214 J = -ctx->n_0*ctx->v_0*PetscRealPart(tt[0]); 215 ierr = PetscDSSetConstants(prob, 1, &constants[0]);CHKERRQ(ierr); 216 ierr = PetscDSSetObjective(prob, 0, &f0_j_re);CHKERRQ(ierr); 217 ierr = DMPlexComputeIntegralFEM(plex,X,tt,NULL);CHKERRQ(ierr); 218 J_re = -ctx->n_0*ctx->v_0*PetscRealPart(tt[0]); 219 ierr = getTe_kev(plex, X, NULL, &Te_kev);CHKERRQ(ierr); 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 ierr = 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);CHKERRQ(ierr); 226 if (done) { 227 ierr = TSSetConvergedReason(ts,TS_CONVERGED_USER);CHKERRQ(ierr); 228 old_ratio = 0; 229 } else { 230 TSConvergedReason reason; 231 ierr = TSGetConvergedReason(ts,&reason);CHKERRQ(ierr); 232 old_ratio = ratio; 233 if (reason) done = PETSC_TRUE; 234 } 235 PetscFunctionReturn(0); 236 } 237 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 } 270 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; 280 281 PetscFunctionBeginUser; 282 ierr = VecGetDM(X, &dm);CHKERRQ(ierr); 283 ierr = DMGetDS(plex, &prob);CHKERRQ(ierr); 284 ierr = VecDuplicate(X,&X2);CHKERRQ(ierr); 285 ierr = VecCopy(X,X2);CHKERRQ(ierr); 286 if (!rectx->X_0) { 287 ierr = VecDuplicate(X,&rectx->X_0);CHKERRQ(ierr); 288 ierr = VecCopy(X,rectx->X_0);CHKERRQ(ierr); 289 } 290 ierr = VecAXPY(X,-1.0,rectx->X_0);CHKERRQ(ierr); 291 ierr = PetscDSSetConstants(prob, sizeof(LandauCtx)/sizeof(PetscScalar), (PetscScalar*)ctx);CHKERRQ(ierr); 292 rectx->idx = 0; 293 ierr = PetscDSSetObjective(prob, 0, &f0_0_diff_lp);CHKERRQ(ierr); 294 ierr = DMPlexComputeIntegralFEM(plex,X2,tt,NULL);CHKERRQ(ierr); 295 ediff = PetscPowReal(PetscRealPart(tt[0]),1./ppp); 296 ierr = PetscDSSetObjective(prob, 0, &f0_0_maxwellian_lp);CHKERRQ(ierr); 297 ierr = DMPlexComputeIntegralFEM(plex,X2,tt,NULL);CHKERRQ(ierr); 298 lpm0 = PetscPowReal(PetscRealPart(tt[0]),1./ppp); 299 if (ctx->num_species>1) { 300 rectx->idx = 1; 301 ierr = PetscDSSetObjective(prob, 0, &f0_0_diff_lp);CHKERRQ(ierr); 302 ierr = DMPlexComputeIntegralFEM(plex,X2,tt,NULL);CHKERRQ(ierr); 303 idiff = PetscPowReal(PetscRealPart(tt[0]),1./ppp); 304 ierr = PetscDSSetObjective(prob, 0, &f0_0_maxwellian_lp);CHKERRQ(ierr); 305 ierr = DMPlexComputeIntegralFEM(plex,X2,tt,NULL);CHKERRQ(ierr); 306 lpm1 = PetscPowReal(PetscRealPart(tt[0]),1./ppp); 307 } 308 ierr = 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);CHKERRQ(ierr); 309 /* view */ 310 ierr = VecViewFromOptions(X,NULL,"-vec_view_diff");CHKERRQ(ierr); 311 ierr = VecCopy(X2,X);CHKERRQ(ierr); 312 ierr = VecDestroy(&X2);CHKERRQ(ierr); 313 if (islast) { 314 ierr = VecDestroy(&rectx->X_0);CHKERRQ(ierr); 315 rectx->X_0 = NULL; 316 } 317 PetscFunctionReturn(0); 318 } 319 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; 330 331 PetscFunctionBeginUser; 332 for (ii=0;ii<ctx->num_species;ii++) constants[ii] = ctx->charges[ii]; 333 ierr = VecGetDM(X, &dm);CHKERRQ(ierr); 334 ierr = DMConvert(dm, DMPLEX, &plex);CHKERRQ(ierr); 335 ierr = DMGetDS(plex, &prob);CHKERRQ(ierr); 336 /* J */ 337 ierr = PetscDSSetConstants(prob, ctx->num_species, constants);CHKERRQ(ierr); 338 ierr = PetscDSSetObjective(prob, 0, &f0_jz_sum);CHKERRQ(ierr); 339 ierr = DMPlexComputeIntegralFEM(plex,X,tt,NULL);CHKERRQ(ierr); 340 J = -ctx->n_0*ctx->v_0*PetscRealPart(tt[0]); 341 ierr = getTe_kev(plex, X, NULL, &Te_kev);CHKERRQ(ierr); 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 ierr = 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");CHKERRQ(ierr); 353 old_ratio = ratio; 354 } else if (rectx->use_spitzer_eta) { 355 /* set E */ 356 *a_E = spit_eta*J; 357 ierr = 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);CHKERRQ(ierr); 358 } else { 359 ierr = 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);CHKERRQ(ierr); 360 } 361 /* cleanup */ 362 ierr = DMDestroy(&plex);CHKERRQ(ierr); 363 PetscFunctionReturn(0); 364 } 365 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; 375 376 PetscFunctionBeginUser; 377 for (ii=0;ii<ctx->num_species;ii++) constants[ii] = ctx->charges[ii]; 378 ierr = VecGetDM(X, &dm);CHKERRQ(ierr); 379 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 380 ierr = DMConvert(dm, DMPLEX, &plex);CHKERRQ(ierr); 381 /* get d current / dt */ 382 ierr = PetscDSSetConstants(prob, ctx->num_species, constants);CHKERRQ(ierr); 383 ierr = PetscDSSetObjective(prob, 0, &f0_jz_sum);CHKERRQ(ierr); 384 if (!X_t) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "X_t"); 385 ierr = DMPlexComputeIntegralFEM(plex,X_t,tt,NULL);CHKERRQ(ierr); 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 ierr = DMDestroy(&plex);CHKERRQ(ierr); 390 PetscFunctionReturn(0); 391 } 392 393 static PetscErrorCode EConstant(Vec X, Vec X_t, PetscInt step, PetscReal time, LandauCtx *ctx, PetscReal *a_E) 394 { 395 PetscFunctionBeginUser; 396 *a_E = ctx->Ez; 397 PetscFunctionReturn(0); 398 } 399 400 static PetscErrorCode ENone(Vec X, Vec X_t, PetscInt step, PetscReal time, LandauCtx *ctx, PetscReal *a_E) 401 { 402 PetscFunctionBeginUser; 403 *a_E = 0; 404 PetscFunctionReturn(0); 405 } 406 407 /* ------------------------------------------------------------------- */ 408 /* 409 FormSource - Evaluates source terms F(t). 410 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() 416 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; 425 PetscErrorCode ierr; 426 REctx *rectx; 427 428 PetscFunctionBeginUser; 429 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 430 ierr = DMGetApplicationContext(dm, &ctx);CHKERRQ(ierr); 431 rectx = (REctx*)ctx->data; 432 /* check for impurities */ 433 ierr = rectx->impuritySrcRate(ftime,&new_imp_rate,ctx);CHKERRQ(ierr); 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 ierr = DMConvert(dm, DMPLEX, &plex);CHKERRQ(ierr); 441 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 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 ierr = 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);CHKERRQ(ierr); 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 ierr = DMCreateGlobalVector(dm, &rectx->imp_src);CHKERRQ(ierr); 452 ierr = PetscObjectSetName((PetscObject)rectx->imp_src, "source");CHKERRQ(ierr); 453 } 454 ierr = VecZeroEntries(rectx->imp_src);CHKERRQ(ierr); 455 ierr = LandauAddMaxwellians(plex,rectx->imp_src,ftime,temps,tilda_ns,ctx);CHKERRQ(ierr); 456 /* clean up */ 457 ierr = DMDestroy(&plex);CHKERRQ(ierr); 458 if (0) { 459 KSP ksp; 460 Vec S; 461 ierr = KSPCreate(PETSC_COMM_SELF,&ksp);CHKERRQ(ierr); 462 ierr = KSPSetOptionsPrefix(ksp,"mass_");CHKERRQ(ierr); /* stokes */ 463 ierr = KSPSetOperators(ksp,ctx->M,ctx->M);CHKERRQ(ierr); 464 ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr); 465 ierr = DMCreateGlobalVector(dm, &S);CHKERRQ(ierr); 466 ierr = KSPSolve(ksp,rectx->imp_src,S);CHKERRQ(ierr); 467 ierr = VecCopy(S,rectx->imp_src);CHKERRQ(ierr); 468 ierr = VecDestroy(&S);CHKERRQ(ierr); 469 } 470 ierr = VecViewFromOptions(rectx->imp_src,NULL,"-vec_view_sources");CHKERRQ(ierr); 471 } 472 ierr = VecCopy(rectx->imp_src,F);CHKERRQ(ierr); 473 } else { 474 if (rectx->current_rate != 0 && rectx->imp_src) { 475 ierr = VecZeroEntries(rectx->imp_src);CHKERRQ(ierr); 476 } 477 rectx->current_rate = 0; 478 } 479 PetscFunctionReturn(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; 489 490 PetscFunctionBeginUser; 491 ierr = VecGetDM(X, &dm);CHKERRQ(ierr); 492 if (stepi > rectx->plotStep && rectx->plotting) { 493 rectx->plotting = PETSC_FALSE; /* was doing diagnostics, now done */ 494 rectx->plotIdx++; 495 } 496 /* view */ 497 ierr = TSGetConvergedReason(ts,&reason);CHKERRQ(ierr); 498 if (time/rectx->plotDt >= (PetscReal)rectx->plotIdx || reason) { 499 /* print norms */ 500 ierr = LandauPrintNorms(X, stepi);CHKERRQ(ierr); 501 ierr = DMConvert(dm, DMPLEX, &plex);CHKERRQ(ierr); 502 ierr = DMGetDS(plex, &prob);CHKERRQ(ierr); 503 /* diagnostics */ 504 ierr = rectx->test(ts,X,plex,stepi,time,reason ? PETSC_TRUE : PETSC_FALSE, ctx,rectx);CHKERRQ(ierr); 505 ierr = DMDestroy(&plex);CHKERRQ(ierr); 506 /* view */ 507 ierr = DMSetOutputSequenceNumber(dm, rectx->plotIdx, time*ctx->t_0);CHKERRQ(ierr); 508 ierr = VecViewFromOptions(X,NULL,"-vec_view");CHKERRQ(ierr); 509 rectx->plotStep = stepi; 510 rectx->plotting = PETSC_TRUE; 511 } 512 /* parallel check */ 513 if (reason) { 514 PetscReal val,rval; 515 PetscMPIInt rank; 516 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank);CHKERRQ(ierr); 517 ierr = TSGetSolution(ts, &X);CHKERRQ(ierr); 518 ierr = VecNorm(X,NORM_2,&val);CHKERRQ(ierr); 519 ierr = MPIU_Allreduce(&val,&rval,1,MPIU_REAL,MPIU_MAX,PETSC_COMM_WORLD);CHKERRQ(ierr); 520 if (rval != val) { 521 ierr = PetscPrintf(PETSC_COMM_SELF, " ***** [%D] ERROR max |x| = %e, my |x| = %20.13e\n",rank,rval,val);CHKERRQ(ierr); 522 } else { 523 ierr = PetscPrintf(PETSC_COMM_WORLD, "[%D] parallel consistency check OK\n",rank);CHKERRQ(ierr); 524 } 525 } 526 rectx->idx = 0; 527 PetscFunctionReturn(0); 528 } 529 530 PetscErrorCode PreStep(TS ts) 531 { 532 PetscErrorCode ierr; 533 LandauCtx *ctx; 534 REctx *rectx; 535 DM dm; 536 PetscInt stepi; 537 PetscReal time; 538 Vec X; 539 540 PetscFunctionBeginUser; 541 /* check seed RE run */ 542 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 543 ierr = TSGetTime(ts,&time);CHKERRQ(ierr); 544 ierr = TSGetSolution(ts,&X);CHKERRQ(ierr); 545 ierr = DMGetApplicationContext(dm, &ctx);CHKERRQ(ierr); 546 rectx = (REctx*)ctx->data; 547 ierr = TSGetStepNumber(ts, &stepi);CHKERRQ(ierr); 548 /* update E */ 549 ierr = rectx->E(X, NULL, stepi, time, ctx, &ctx->Ez);CHKERRQ(ierr); 550 PetscFunctionReturn(0); 551 } 552 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; 557 558 PetscFunctionBeginUser; 559 if (time >= rectx->pulse_start) *rho = rectx->pulse_rate; 560 else *rho = 0.; 561 PetscFunctionReturn(0); 562 } 563 static PetscErrorCode zeroSrc(PetscReal time, PetscReal *rho, LandauCtx *ctx) 564 { 565 PetscFunctionBeginUser; 566 *rho = 0.; 567 PetscFunctionReturn(0); 568 } 569 static PetscErrorCode pulseSrc(PetscReal time, PetscReal *rho, LandauCtx *ctx) 570 { 571 REctx *rectx = (REctx*)ctx->data; 572 573 PetscFunctionBeginUser; 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 PetscFunctionReturn(0); 587 } 588 589 #undef __FUNCT__ 590 #define __FUNCT__ "ProcessREOptions" 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; 597 598 PetscFunctionBeginUser; 599 ierr = DMCreate(PETSC_COMM_WORLD,&dummy);CHKERRQ(ierr); 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 ierr = PetscFunctionListAdd(&plist,"step",&stepSrc);CHKERRQ(ierr); 620 ierr = PetscFunctionListAdd(&plist,"none",&zeroSrc);CHKERRQ(ierr); 621 ierr = PetscFunctionListAdd(&plist,"pulse",&pulseSrc);CHKERRQ(ierr); 622 ierr = PetscStrcpy(pname,"none");CHKERRQ(ierr); 623 ierr = PetscFunctionListAdd(&testlist,"none",&testNone);CHKERRQ(ierr); 624 ierr = PetscFunctionListAdd(&testlist,"spitzer",&testSpitzer);CHKERRQ(ierr); 625 ierr = PetscFunctionListAdd(&testlist,"stable",&testStable);CHKERRQ(ierr); 626 ierr = PetscStrcpy(testname,"none");CHKERRQ(ierr); 627 ierr = PetscFunctionListAdd(&elist,"none",&ENone);CHKERRQ(ierr); 628 ierr = PetscFunctionListAdd(&elist,"spitzer",&ESpitzer);CHKERRQ(ierr); 629 ierr = PetscFunctionListAdd(&elist,"induction",&EInduction);CHKERRQ(ierr); 630 ierr = PetscFunctionListAdd(&elist,"constant",&EConstant);CHKERRQ(ierr); 631 ierr = PetscStrcpy(ename,"constant");CHKERRQ(ierr); 632 633 ierr = PetscOptionsBegin(PETSC_COMM_SELF, prefix, "Options for Runaway/seed electron model", "none");CHKERRQ(ierr); 634 ierr = PetscOptionsReal("-ex2_plot_dt", "Plotting interval", "xgc_dmplex.c", rectx->plotDt, &rectx->plotDt, NULL);CHKERRQ(ierr); 635 if (rectx->plotDt < 0) rectx->plotDt = 1e30; 636 if (rectx->plotDt == 0) rectx->plotDt = 1e-30; 637 ierr = PetscOptionsFList("-ex2_impurity_source_type","Name of impurity source to run","",plist,pname,pname,sizeof(pname),NULL);CHKERRQ(ierr); 638 ierr = PetscOptionsFList("-ex2_test_type","Name of test to run","",testlist,testname,testname,sizeof(testname),NULL);CHKERRQ(ierr); 639 ierr = PetscOptionsInt("-ex2_impurity_index", "index of sink for impurities", "none", rectx->imp_idx, &rectx->imp_idx, NULL);CHKERRQ(ierr); 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 ierr = PetscOptionsFList("-ex2_e_field_type","Electric field type","",elist,ename,ename,sizeof(ename),NULL);CHKERRQ(ierr); 642 rectx->Ne_ion = -ctx->charges[rectx->imp_idx]/ctx->charges[0]; 643 ierr = PetscOptionsReal("-ex2_t_cold","Temperature of cold electron and ions after ionization in keV","none",rectx->T_cold,&rectx->T_cold, NULL);CHKERRQ(ierr); 644 ierr = PetscOptionsReal("-ex2_pulse_start_time","Time at which pulse happens for 'pulse' source","none",rectx->pulse_start,&rectx->pulse_start, NULL);CHKERRQ(ierr); 645 ierr = PetscOptionsReal("-ex2_pulse_width_time","Width of pulse 'pulse' source","none",rectx->pulse_width,&rectx->pulse_width, NULL);CHKERRQ(ierr); 646 ierr = PetscOptionsReal("-ex2_pulse_rate","Number density of pulse for 'pulse' source","none",rectx->pulse_rate,&rectx->pulse_rate, NULL);CHKERRQ(ierr); 647 rectx->T_cold *= 1.16e7; /* convert to Kelvin */ 648 ierr = PetscOptionsReal("-ex2_ion_potential","Potential to ionize impurity (should be array) in ev","none",rectx->ion_potential,&rectx->ion_potential, NULL);CHKERRQ(ierr); 649 ierr = PetscOptionsReal("-ex2_inductance","","none",rectx->L,&rectx->L, NULL);CHKERRQ(ierr); 650 ierr = 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);CHKERRQ(ierr); 651 ierr = PetscOptionsEnd();CHKERRQ(ierr); 652 /* get impurity source rate function */ 653 ierr = PetscFunctionListFind(plist,pname,&rectx->impuritySrcRate);CHKERRQ(ierr); 654 if (!rectx->impuritySrcRate) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"No impurity source function found '%s'",pname); 655 ierr = PetscFunctionListFind(testlist,testname,&rectx->test);CHKERRQ(ierr); 656 if (!rectx->test) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"No test found '%s'",testname); 657 ierr = PetscFunctionListFind(elist,ename,&rectx->E);CHKERRQ(ierr); 658 if (!rectx->E) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"No E field function found '%s'",ename); 659 ierr = PetscFunctionListDestroy(&plist);CHKERRQ(ierr); 660 ierr = PetscFunctionListDestroy(&testlist);CHKERRQ(ierr); 661 ierr = PetscFunctionListDestroy(&elist);CHKERRQ(ierr); 662 { 663 PetscMPIInt rank; 664 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank);CHKERRQ(ierr); 665 if (rank) { /* turn off output stuff for duplicate runs */ 666 ierr = PetscOptionsClearValue(NULL,"-dm_view");CHKERRQ(ierr); 667 ierr = PetscOptionsClearValue(NULL,"-vec_view");CHKERRQ(ierr); 668 ierr = PetscOptionsClearValue(NULL,"-dm_view_diff");CHKERRQ(ierr); 669 ierr = PetscOptionsClearValue(NULL,"-vec_view_diff");CHKERRQ(ierr); 670 ierr = PetscOptionsClearValue(NULL,"-dm_view_sources");CHKERRQ(ierr); 671 ierr = PetscOptionsClearValue(NULL,"-vec_view_sources");CHKERRQ(ierr); 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 ierr = DMDestroy(&dummy);CHKERRQ(ierr); 681 PetscFunctionReturn(0); 682 } 683 684 int main(int argc, char **argv) 685 { 686 DM dm; 687 Vec X; 688 PetscErrorCode ierr; 689 PetscInt dim = 2; 690 TS ts; 691 Mat J; 692 PetscDS prob; 693 LandauCtx *ctx; 694 REctx *rectx; 695 696 ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr; 697 ierr = PetscOptionsGetInt(NULL,NULL, "-dim", &dim, NULL);CHKERRQ(ierr); 698 /* Create a mesh */ 699 ierr = LandauCreateVelocitySpace(PETSC_COMM_SELF, dim, "", &X, &J, &dm);CHKERRQ(ierr); 700 ierr = LandauCreateMassMatrix(dm, NULL);CHKERRQ(ierr); 701 ierr = DMGetApplicationContext(dm, &ctx);CHKERRQ(ierr); 702 ierr = DMSetUp(dm);CHKERRQ(ierr); 703 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 704 s_quarter3DDomain_notused = ctx->quarter3DDomain; /* ugh */ 705 /* context */ 706 ierr = PetscNew(&rectx);CHKERRQ(ierr); 707 ctx->data = rectx; 708 ierr = ProcessREOptions(rectx,ctx,dm,"");CHKERRQ(ierr); 709 ierr = DMSetOutputSequenceNumber(dm, 0, 0.0);CHKERRQ(ierr); 710 ierr = DMViewFromOptions(dm,NULL,"-dm_view");CHKERRQ(ierr); 711 ierr = DMViewFromOptions(dm,NULL,"-dm_view_sources");CHKERRQ(ierr); 712 ierr = DMViewFromOptions(dm,NULL,"-dm_view_diff");CHKERRQ(ierr); 713 /* Create timestepping solver context */ 714 ierr = TSCreate(PETSC_COMM_SELF,&ts);CHKERRQ(ierr); 715 ierr = TSSetDM(ts,dm);CHKERRQ(ierr); 716 ierr = TSSetIFunction(ts,NULL,LandauIFunction,NULL);CHKERRQ(ierr); 717 ierr = TSSetIJacobian(ts,J,J,LandauIJacobian,NULL);CHKERRQ(ierr); 718 ierr = TSSetRHSFunction(ts,NULL,FormSource,NULL);CHKERRQ(ierr); 719 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 720 ierr = TSSetSolution(ts,X);CHKERRQ(ierr); 721 ierr = TSSetApplicationContext(ts, ctx);CHKERRQ(ierr); 722 ierr = TSMonitorSet(ts,Monitor,ctx,NULL);CHKERRQ(ierr); 723 ierr = TSSetPreStep(ts,PreStep);CHKERRQ(ierr); 724 rectx->Ez_initial = ctx->Ez; /* cache for induction caclulation - applied E field */ 725 ierr = MatSetOption(J, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE);CHKERRQ(ierr); 726 { /* warm up an test just LandauIJacobian */ 727 #if defined(PETSC_USE_LOG) 728 PetscLogStage stage; 729 #endif 730 Vec vec; 731 PetscRandom rctx; 732 PetscInt ii; 733 734 ierr = PetscRandomCreate(PETSC_COMM_SELF,&rctx);CHKERRQ(ierr); 735 ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr); 736 ierr = VecDuplicate(X,&vec);CHKERRQ(ierr); 737 /* warm up */ 738 ierr = VecSetRandom(vec,rctx);CHKERRQ(ierr); 739 ierr = PetscLogStageRegister("Warmup", &stage);CHKERRQ(ierr); 740 ierr = PetscLogStagePush(stage);CHKERRQ(ierr); 741 ierr = LandauIJacobian(ts,0.0,vec,vec,1.0,J,J,ctx);CHKERRQ(ierr); 742 ierr = PetscLogStagePop();CHKERRQ(ierr); 743 /* LandauIJacobian */ 744 ierr = PetscLogStageRegister("LandauIJacobian", &stage);CHKERRQ(ierr); 745 ierr = PetscLogStagePush(stage);CHKERRQ(ierr); 746 for (ii=0;ii<10;ii++){ 747 ierr = VecSetRandom(vec,rctx);CHKERRQ(ierr); 748 ierr = LandauIJacobian(ts,0.0,vec,vec,1.0,J,J,ctx);CHKERRQ(ierr); 749 } 750 ierr = PetscLogStagePop();CHKERRQ(ierr); 751 ierr = VecDestroy(&vec);CHKERRQ(ierr); 752 ierr = PetscRandomDestroy(&rctx);CHKERRQ(ierr); 753 } 754 ierr = VecViewFromOptions(X,NULL,"-vec_view");CHKERRQ(ierr); // inital condition (monitor plots after step) 755 /* go */ 756 ierr = TSSolve(ts,X);CHKERRQ(ierr); 757 /* clean up */ 758 ierr = LandauDestroyVelocitySpace(&dm);CHKERRQ(ierr); 759 ierr = TSDestroy(&ts);CHKERRQ(ierr); 760 ierr = VecDestroy(&X);CHKERRQ(ierr); 761 if (rectx->imp_src) { 762 ierr = VecDestroy(&rectx->imp_src);CHKERRQ(ierr); 763 } 764 ierr = PetscFree(rectx);CHKERRQ(ierr); 765 ierr = PetscFinalize(); 766 return ierr; 767 } 768 769 /*TEST 770 test: 771 suffix: 0 772 requires: p4est !complex double 773 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 774 775 TEST*/ 776