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 #include <petscdmcomposite.h> 8 9 /* data for runaway electron model */ 10 typedef struct REctx_struct { 11 PetscErrorCode (*test)(TS, Vec, PetscInt, PetscReal, PetscBool, LandauCtx *, struct REctx_struct *); 12 PetscErrorCode (*impuritySrcRate)(PetscReal, PetscReal *, LandauCtx*); 13 PetscErrorCode (*E)(Vec, Vec, PetscInt, PetscReal, LandauCtx*, PetscReal *); 14 PetscReal T_cold; /* temperature of newly ionized electrons and impurity ions */ 15 PetscReal ion_potential; /* ionization potential of impurity */ 16 PetscReal Ne_ion; /* effective number of electrons shed in ioization of impurity */ 17 PetscReal Ez_initial; 18 PetscReal L; /* inductance */ 19 Vec X_0; 20 PetscInt imp_idx; /* index for impurity ionizing sink */ 21 PetscReal pulse_start; 22 PetscReal pulse_width; 23 PetscReal pulse_rate; 24 PetscReal current_rate; 25 PetscInt plotIdx; 26 PetscInt plotStep; 27 PetscInt idx; /* cache */ 28 PetscReal j; /* cache */ 29 PetscReal plotDt; 30 PetscBool plotting; 31 PetscBool use_spitzer_eta; 32 PetscInt print_period; 33 PetscInt grid_view_idx; 34 } REctx; 35 36 static const PetscReal kev_joul = 6.241506479963235e+15; /* 1/1000e */ 37 38 #define RE_CUT 3. 39 /* < v, u_re * v * q > */ 40 static void f0_j_re(PetscInt dim, PetscInt Nf, PetscInt NfAux, 41 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 42 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 43 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar *f0) 44 { 45 PetscReal n_e = PetscRealPart(u[0]); 46 if (dim==2) { 47 if (x[1] > RE_CUT || x[1] < -RE_CUT) { /* simply a cutoff for REs. v_|| > 3 v(T_e) */ 48 *f0 = n_e * 2.*PETSC_PI*x[0] * x[1] * constants[0]; /* n * r * v_|| * q */ 49 } else { 50 *f0 = 0; 51 } 52 } else { 53 if (x[2] > RE_CUT || x[2] < -RE_CUT) { /* simply a cutoff for REs. v_|| > 3 v(T_e) */ 54 *f0 = n_e * x[2] * constants[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 q[], PetscScalar *f0) 66 { 67 PetscInt ii; 68 f0[0] = 0; 69 if (dim==2) { 70 for (ii=0;ii<Nf;ii++) f0[0] += u[ii] * 2.*PETSC_PI*x[0] * x[1] * q[ii]; /* n * r * v_|| * q * v_0 */ 71 } else { 72 for (ii=0;ii<Nf;ii++) f0[0] += u[ii] * x[2] * q[ii]; /* n * v_|| * q * v_0 */ 73 } 74 } 75 76 /* < v, n_e > */ 77 static void f0_n(PetscInt dim, PetscInt Nf, PetscInt NfAux, 78 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 79 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 80 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar *f0) 81 { 82 PetscInt ii = (PetscInt)PetscRealPart(constants[0]); 83 if (dim==2) f0[0] = 2.*PETSC_PI*x[0]*u[ii]; 84 else { 85 f0[0] = u[ii]; 86 } 87 } 88 89 /* < v, n_e v_|| > */ 90 static void f0_vz(PetscInt dim, PetscInt Nf, PetscInt NfAux, 91 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 92 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 93 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar *f0) 94 { 95 PetscInt ii = (PetscInt)PetscRealPart(constants[0]); 96 if (dim==2) f0[0] = u[ii] * 2.*PETSC_PI*x[0] * x[1]; /* n r v_|| */ 97 else { 98 f0[0] = u[ii] * x[2]; /* n v_|| */ 99 } 100 } 101 102 /* < v, n_e (v-shift) > */ 103 static void f0_ve_shift(PetscInt dim, PetscInt Nf, PetscInt NfAux, 104 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 105 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 106 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar *f0) 107 { 108 PetscReal vz = numConstants>0 ? PetscRealPart(constants[0]) : 0; 109 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 */ 110 else { 111 *f0 = u[0] * PetscSqrtReal(x[0]*x[0] + x[1]*x[1] + (x[2]-vz)*(x[2]-vz)); /* n v */ 112 } 113 } 114 115 /* CalculateE - Calculate the electric field */ 116 /* T -- Electron temperature */ 117 /* n -- Electron density */ 118 /* lnLambda -- */ 119 /* eps0 -- */ 120 /* E -- output E, input \hat E */ 121 static PetscReal CalculateE(PetscReal Tev, PetscReal n, PetscReal lnLambda, PetscReal eps0, PetscReal *E) 122 { 123 PetscReal c,e,m; 124 125 PetscFunctionBegin; 126 c = 299792458.0; 127 e = 1.602176e-19; 128 m = 9.10938e-31; 129 if (1) { 130 double Ec, Ehat = *E, betath = PetscSqrtReal(2*Tev*e/(m*c*c)), j0 = Ehat * 7/(PetscSqrtReal(2)*2) * PetscPowReal(betath,3) * n * e * c; 131 Ec = n*lnLambda*PetscPowReal(e,3) / (4*PETSC_PI*PetscPowReal(eps0,2)*m*c*c); 132 *E = Ec; 133 PetscPrintf(PETSC_COMM_WORLD, "CalculateE j0=%g Ec = %g\n",j0,Ec); 134 } else { 135 PetscReal Ed, vth; 136 vth = PetscSqrtReal(8*Tev*e/(m*PETSC_PI)); 137 Ed = n*lnLambda*PetscPowReal(e,3) / (4*PETSC_PI*PetscPowReal(eps0,2)*m*vth*vth); 138 *E = Ed; 139 } 140 PetscFunctionReturn(0); 141 } 142 143 static PetscReal Spitzer(PetscReal m_e, PetscReal e, PetscReal Z, PetscReal epsilon0, PetscReal lnLam, PetscReal kTe_joules) 144 { 145 PetscReal Fz = (1+1.198*Z+0.222*Z*Z)/(1+2.966*Z+0.753*Z*Z), eta; 146 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); 147 return eta; 148 } 149 150 /* */ 151 static PetscErrorCode testNone(TS ts, Vec X, PetscInt stepi, PetscReal time, PetscBool islast, LandauCtx *ctx, REctx *rectx) 152 { 153 PetscFunctionBeginUser; 154 PetscFunctionReturn(0); 155 } 156 157 /* */ 158 static PetscErrorCode testSpitzer(TS ts, Vec X, PetscInt stepi, PetscReal time, PetscBool islast, LandauCtx *ctx, REctx *rectx) 159 { 160 PetscErrorCode ierr; 161 PetscInt ii,nDMs; 162 PetscDS prob; 163 static PetscReal old_ratio = 1e10; 164 TSConvergedReason reason; 165 PetscReal J,J_re,spit_eta,Te_kev=0,E,ratio,Z,n_e,v,v2; 166 PetscScalar user[2] = {0.,ctx->charges[0]}, q[LANDAU_MAX_SPECIES],tt[LANDAU_MAX_SPECIES],vz; 167 PetscReal dt; 168 DM pack, plexe = ctx->plex[0], plexi = (ctx->num_grids==1) ? NULL : ctx->plex[1]; 169 Vec *XsubArray; 170 171 PetscFunctionBeginUser; 172 if (ctx->num_species!=2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "ctx->num_species %D != 2",ctx->num_species); 173 ierr = VecGetDM(X, &pack);CHKERRQ(ierr); 174 if (!pack) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "no DM"); 175 ierr = DMCompositeGetNumberDM(pack,&nDMs);CHKERRQ(ierr); 176 if (nDMs != ctx->num_grids*ctx->batch_sz) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "nDMs != ctx->num_grids*ctx->batch_sz %D != %D",nDMs,ctx->num_grids*ctx->batch_sz); 177 ierr = PetscMalloc(sizeof(*XsubArray)*nDMs, &XsubArray);CHKERRQ(ierr); 178 ierr = DMCompositeGetAccessArray(pack, X, nDMs, NULL, XsubArray);CHKERRQ(ierr); // read only 179 ierr = TSGetTimeStep(ts,&dt);CHKERRQ(ierr); 180 /* get current for each grid */ 181 for (ii=0;ii<ctx->num_species;ii++) q[ii] = ctx->charges[ii]; 182 ierr = DMGetDS(plexe, &prob);CHKERRQ(ierr); 183 ierr = PetscDSSetConstants(prob, 2, &q[0]);CHKERRQ(ierr); 184 ierr = PetscDSSetObjective(prob, 0, &f0_jz_sum);CHKERRQ(ierr); 185 ierr = DMPlexComputeIntegralFEM(plexe,XsubArray[ LAND_PACK_IDX(ctx->batch_view_idx,0) ],tt,NULL);CHKERRQ(ierr); 186 J = -ctx->n_0*ctx->v_0*PetscRealPart(tt[0]); 187 if (plexi) { // add first (only) ion 188 ierr = DMGetDS(plexi, &prob);CHKERRQ(ierr); 189 ierr = PetscDSSetConstants(prob, 1, &q[1]);CHKERRQ(ierr); 190 ierr = PetscDSSetObjective(prob, 0, &f0_jz_sum);CHKERRQ(ierr); 191 ierr = DMPlexComputeIntegralFEM(plexi,XsubArray[LAND_PACK_IDX(ctx->batch_view_idx,1)],tt,NULL);CHKERRQ(ierr); 192 J += -ctx->n_0*ctx->v_0*PetscRealPart(tt[0]); 193 } 194 /* get N_e */ 195 ierr = DMGetDS(plexe, &prob);CHKERRQ(ierr); 196 ierr = PetscDSSetConstants(prob, 1, user);CHKERRQ(ierr); 197 ierr = PetscDSSetObjective(prob, 0, &f0_n);CHKERRQ(ierr); 198 ierr = DMPlexComputeIntegralFEM(plexe,XsubArray[LAND_PACK_IDX(ctx->batch_view_idx,0)],tt,NULL);CHKERRQ(ierr); 199 n_e = PetscRealPart(tt[0])*ctx->n_0; 200 /* Z */ 201 Z = -ctx->charges[1]/ctx->charges[0]; 202 /* remove drift */ 203 if (0) { 204 user[0] = 0; // electrons 205 ierr = DMGetDS(plexe, &prob);CHKERRQ(ierr); 206 ierr = PetscDSSetConstants(prob, 1, user);CHKERRQ(ierr); 207 ierr = PetscDSSetObjective(prob, 0, &f0_vz);CHKERRQ(ierr); 208 ierr = DMPlexComputeIntegralFEM(plexe,XsubArray[LAND_PACK_IDX(ctx->batch_view_idx,0)],tt,NULL);CHKERRQ(ierr); 209 vz = ctx->n_0*PetscRealPart(tt[0])/n_e; /* non-dimensional */ 210 } else vz = 0; 211 /* thermal velocity */ 212 ierr = DMGetDS(plexe, &prob);CHKERRQ(ierr); 213 ierr = PetscDSSetConstants(prob, 1, &vz);CHKERRQ(ierr); 214 ierr = PetscDSSetObjective(prob, 0, &f0_ve_shift);CHKERRQ(ierr); 215 ierr = DMPlexComputeIntegralFEM(plexe,XsubArray[LAND_PACK_IDX(ctx->batch_view_idx,0)],tt,NULL);CHKERRQ(ierr); 216 v = ctx->n_0*ctx->v_0*PetscRealPart(tt[0])/n_e; /* remove number density to get velocity */ 217 v2 = PetscSqr(v); /* use real space: m^2 / s^2 */ 218 Te_kev = (v2*ctx->masses[0]*PETSC_PI/8)*kev_joul; /* temperature in kev */ 219 //Te_kev = ctx->thermal_temps[0]/1.1604525e7; 220 spit_eta = Spitzer(ctx->masses[0],-ctx->charges[0],Z,ctx->epsilon0,ctx->lnLam,Te_kev/kev_joul); /* kev --> J (kT) */ 221 if (0) { 222 ierr = DMGetDS(plexe, &prob);CHKERRQ(ierr); 223 ierr = PetscDSSetConstants(prob, 1, q);CHKERRQ(ierr); 224 ierr = PetscDSSetObjective(prob, 0, &f0_j_re);CHKERRQ(ierr); 225 ierr = DMPlexComputeIntegralFEM(plexe,XsubArray[LAND_PACK_IDX(ctx->batch_view_idx,0)],tt,NULL);CHKERRQ(ierr); 226 } else tt[0] = 0; 227 J_re = -ctx->n_0*ctx->v_0*PetscRealPart(tt[0]); 228 ierr = DMCompositeRestoreAccessArray(pack, X, nDMs, NULL, XsubArray);CHKERRQ(ierr); // read only 229 ierr = PetscFree(XsubArray);CHKERRQ(ierr); 230 231 if (rectx->use_spitzer_eta) { 232 E = ctx->Ez = spit_eta*(rectx->j-J_re); 233 } else { 234 E = ctx->Ez; /* keep real E */ 235 rectx->j = J; /* cache */ 236 } 237 238 ratio = E/J/spit_eta; 239 if (stepi>10 && !rectx->use_spitzer_eta && ( 240 //(old_ratio-ratio < 1.e-3 && ratio > 0.99 && ratio < 1.01) || 241 //(old_ratio-ratio < 1.e-4 && ratio > 0.98 && ratio < 1.02) || 242 (old_ratio-ratio < 1.e-6))) { 243 rectx->pulse_start = time + 0.98*dt; 244 rectx->use_spitzer_eta = PETSC_TRUE; 245 } 246 ierr = TSGetConvergedReason(ts,&reason);CHKERRQ(ierr); 247 ierr = TSGetConvergedReason(ts,&reason);CHKERRQ(ierr); 248 if ((rectx->plotting) || stepi == 0 || reason || rectx->pulse_start == time + 0.98*dt) { 249 ierr = 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 spit_eta=%g\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",spit_eta);CHKERRQ(ierr); 250 if (rectx->pulse_start == time + 0.98*dt) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Spitzer complete ratio=%g",ratio); 251 } 252 old_ratio = ratio; 253 PetscFunctionReturn(0); 254 } 255 256 static const double ppp = 2; 257 static void f0_0_diff_lp(PetscInt dim, PetscInt Nf, PetscInt NfAux, 258 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 259 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 260 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar *f0) 261 { 262 LandauCtx *ctx = (LandauCtx *)constants; 263 REctx *rectx = (REctx*)ctx->data; 264 PetscInt ii = rectx->idx, i; 265 const PetscReal kT_m = ctx->k*ctx->thermal_temps[ii]/ctx->masses[ii]; /* kT/m */ 266 const PetscReal n = ctx->n[ii]; 267 PetscReal diff, f_maxwell, v2 = 0, theta = 2*kT_m/(ctx->v_0*ctx->v_0); /* theta = 2kT/mc^2 */ 268 for (i = 0; i < dim; ++i) v2 += x[i]*x[i]; 269 f_maxwell = n*PetscPowReal(PETSC_PI*theta,-1.5)*(PetscExpReal(-v2/theta)); 270 diff = 2.*PETSC_PI*x[0]*(PetscRealPart(u[ii]) - f_maxwell); 271 f0[0] = PetscPowReal(diff,ppp); 272 } 273 static void f0_0_maxwellian_lp(PetscInt dim, PetscInt Nf, PetscInt NfAux, 274 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 275 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 276 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar *f0) 277 { 278 LandauCtx *ctx = (LandauCtx *)constants; 279 REctx *rectx = (REctx*)ctx->data; 280 PetscInt ii = rectx->idx, i; 281 const PetscReal kT_m = ctx->k*ctx->thermal_temps[ii]/ctx->masses[ii]; /* kT/m */ 282 const PetscReal n = ctx->n[ii]; 283 PetscReal f_maxwell, v2 = 0, theta = 2*kT_m/(ctx->v_0*ctx->v_0); /* theta = 2kT/mc^2 */ 284 for (i = 0; i < dim; ++i) v2 += x[i]*x[i]; 285 f_maxwell = 2.*PETSC_PI*x[0] * n*PetscPowReal(PETSC_PI*theta,-1.5)*(PetscExpReal(-v2/theta)); 286 f0[0] = PetscPowReal(f_maxwell,ppp); 287 } 288 289 /* */ 290 static PetscErrorCode testStable(TS ts, Vec X, PetscInt stepi, PetscReal time, PetscBool islast, LandauCtx *ctx, REctx *rectx) 291 { 292 PetscErrorCode ierr; 293 PetscDS prob; 294 Vec X2; 295 PetscReal ediff,idiff=0,lpm0,lpm1=1; 296 PetscScalar tt[LANDAU_MAX_SPECIES]; 297 DM dm, plex = ctx->plex[0]; 298 299 PetscFunctionBeginUser; 300 ierr = VecGetDM(X, &dm);CHKERRQ(ierr); 301 ierr = DMGetDS(plex, &prob);CHKERRQ(ierr); 302 ierr = VecDuplicate(X,&X2);CHKERRQ(ierr); 303 ierr = VecCopy(X,X2);CHKERRQ(ierr); 304 if (!rectx->X_0) { 305 ierr = VecDuplicate(X,&rectx->X_0);CHKERRQ(ierr); 306 ierr = VecCopy(X,rectx->X_0);CHKERRQ(ierr); 307 } 308 ierr = VecAXPY(X,-1.0,rectx->X_0);CHKERRQ(ierr); 309 ierr = PetscDSSetConstants(prob, sizeof(LandauCtx)/sizeof(PetscScalar), (PetscScalar*)ctx);CHKERRQ(ierr); 310 rectx->idx = 0; 311 ierr = PetscDSSetObjective(prob, 0, &f0_0_diff_lp);CHKERRQ(ierr); 312 ierr = DMPlexComputeIntegralFEM(plex,X2,tt,NULL);CHKERRQ(ierr); 313 ediff = PetscPowReal(PetscRealPart(tt[0]),1./ppp); 314 ierr = PetscDSSetObjective(prob, 0, &f0_0_maxwellian_lp);CHKERRQ(ierr); 315 ierr = DMPlexComputeIntegralFEM(plex,X2,tt,NULL);CHKERRQ(ierr); 316 lpm0 = PetscPowReal(PetscRealPart(tt[0]),1./ppp); 317 if (ctx->num_species>1) { 318 rectx->idx = 1; 319 ierr = PetscDSSetObjective(prob, 0, &f0_0_diff_lp);CHKERRQ(ierr); 320 ierr = DMPlexComputeIntegralFEM(plex,X2,tt,NULL);CHKERRQ(ierr); 321 idiff = PetscPowReal(PetscRealPart(tt[0]),1./ppp); 322 ierr = PetscDSSetObjective(prob, 0, &f0_0_maxwellian_lp);CHKERRQ(ierr); 323 ierr = DMPlexComputeIntegralFEM(plex,X2,tt,NULL);CHKERRQ(ierr); 324 lpm1 = PetscPowReal(PetscRealPart(tt[0]),1./ppp); 325 } 326 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); 327 /* view */ 328 ierr = VecCopy(X2,X);CHKERRQ(ierr); 329 ierr = VecDestroy(&X2);CHKERRQ(ierr); 330 if (islast) { 331 ierr = VecDestroy(&rectx->X_0);CHKERRQ(ierr); 332 rectx->X_0 = NULL; 333 } 334 PetscFunctionReturn(0); 335 } 336 337 static PetscErrorCode EInduction(Vec X, Vec X_t, PetscInt step, PetscReal time, LandauCtx *ctx, PetscReal *a_E) 338 { 339 REctx *rectx = (REctx*)ctx->data; 340 PetscErrorCode ierr; 341 PetscInt ii; 342 DM dm,plex; 343 PetscScalar tt[LANDAU_MAX_SPECIES], qv0[LANDAU_MAX_SPECIES]; 344 PetscReal dJ_dt; 345 PetscDS prob; 346 347 PetscFunctionBeginUser; 348 for (ii=0;ii<ctx->num_species;ii++) qv0[ii] = ctx->charges[ii]*ctx->v_0; 349 ierr = VecGetDM(X, &dm);CHKERRQ(ierr); 350 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 351 ierr = DMConvert(dm, DMPLEX, &plex);CHKERRQ(ierr); 352 /* get d current / dt */ 353 ierr = PetscDSSetConstants(prob, ctx->num_species, qv0);CHKERRQ(ierr); 354 ierr = PetscDSSetObjective(prob, 0, &f0_jz_sum);CHKERRQ(ierr); 355 if (!X_t) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "X_t"); 356 ierr = DMPlexComputeIntegralFEM(plex,X_t,tt,NULL);CHKERRQ(ierr); 357 dJ_dt = -ctx->n_0*PetscRealPart(tt[0])/ctx->t_0; 358 /* E induction */ 359 *a_E = -rectx->L*dJ_dt + rectx->Ez_initial; 360 ierr = DMDestroy(&plex);CHKERRQ(ierr); 361 PetscFunctionReturn(0); 362 } 363 364 static PetscErrorCode EConstant(Vec X, Vec X_t, PetscInt step, PetscReal time, LandauCtx *ctx, PetscReal *a_E) 365 { 366 PetscFunctionBeginUser; 367 *a_E = ctx->Ez; 368 PetscFunctionReturn(0); 369 } 370 371 static PetscErrorCode ENone(Vec X, Vec X_t, PetscInt step, PetscReal time, LandauCtx *ctx, PetscReal *a_E) 372 { 373 PetscFunctionBeginUser; 374 *a_E = 0; 375 PetscFunctionReturn(0); 376 } 377 378 /* ------------------------------------------------------------------- */ 379 /* 380 FormSource - Evaluates source terms F(t). 381 382 Input Parameters: 383 . ts - the TS context 384 . time - 385 . X_dummmy - input vector 386 . dummy - optional user-defined context, as set by SNESSetFunction() 387 388 Output Parameter: 389 . F - function vector 390 */ 391 static PetscErrorCode FormSource(TS ts, PetscReal ftime, Vec X_dummmy, Vec F, void *dummy) 392 { 393 PetscReal new_imp_rate; 394 LandauCtx *ctx; 395 DM pack; 396 PetscErrorCode ierr; 397 REctx *rectx; 398 399 PetscFunctionBeginUser; 400 ierr = TSGetDM(ts,&pack);CHKERRQ(ierr); 401 ierr = DMGetApplicationContext(pack, &ctx);CHKERRQ(ierr); 402 rectx = (REctx*)ctx->data; 403 /* check for impurities */ 404 ierr = rectx->impuritySrcRate(ftime,&new_imp_rate,ctx);CHKERRQ(ierr); 405 if (new_imp_rate != 0) { 406 if (new_imp_rate != rectx->current_rate) { 407 PetscInt ii; 408 PetscReal dne_dt,dni_dt,tilda_ns[LANDAU_MAX_SPECIES],temps[LANDAU_MAX_SPECIES]; 409 Vec globFarray[LANDAU_MAX_GRIDS*LANDAU_MAX_BATCH_SZ]; 410 rectx->current_rate = new_imp_rate; 411 for (ii=1;ii<LANDAU_MAX_SPECIES;ii++) tilda_ns[ii] = 0; 412 for (ii=1;ii<LANDAU_MAX_SPECIES;ii++) temps[ii] = 1; 413 dni_dt = new_imp_rate /* *ctx->t_0 */; /* fully ionized immediately, no normalize, stay in non-dim */ 414 dne_dt = new_imp_rate*rectx->Ne_ion /* *ctx->t_0 */; 415 tilda_ns[0] = dne_dt; tilda_ns[rectx->imp_idx] = dni_dt; 416 temps[0] = rectx->T_cold; temps[rectx->imp_idx] = rectx->T_cold; 417 ierr = PetscInfo4(ctx->plex[0], "\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);CHKERRQ(ierr); 418 ierr = DMCompositeGetAccessArray(pack, F, ctx->num_grids*ctx->batch_sz, NULL, globFarray);CHKERRQ(ierr); 419 for (PetscInt grid=0 ; grid<ctx->num_grids ; grid++) { 420 /* add it */ 421 ierr = LandauAddMaxwellians(ctx->plex[grid],globFarray[ LAND_PACK_IDX(0,grid) ],ftime,temps,tilda_ns,grid,0,ctx);CHKERRQ(ierr); 422 ierr = VecViewFromOptions(globFarray[ LAND_PACK_IDX(0,grid) ],NULL,"-vec_view_sources");CHKERRQ(ierr); 423 } 424 // Does DMCompositeRestoreAccessArray copy the data back? (no) 425 ierr = DMCompositeRestoreAccessArray(pack, F, ctx->num_grids*ctx->batch_sz, NULL, globFarray);CHKERRQ(ierr); 426 } 427 } else { 428 ierr = VecZeroEntries(F);CHKERRQ(ierr); 429 rectx->current_rate = 0; 430 } 431 PetscFunctionReturn(0); 432 } 433 PetscErrorCode Monitor(TS ts, PetscInt stepi, PetscReal time, Vec X, void *actx) 434 { 435 LandauCtx *ctx = (LandauCtx*) actx; /* user-defined application context */ 436 REctx *rectx = (REctx*)ctx->data; 437 DM pack; 438 Vec globXArray[LANDAU_MAX_GRIDS*LANDAU_MAX_BATCH_SZ]; 439 TSConvergedReason reason; 440 PetscErrorCode ierr; 441 PetscFunctionBeginUser; 442 ierr = VecGetDM(X, &pack);CHKERRQ(ierr); 443 ierr = DMCompositeGetAccessArray(pack, X, ctx->num_grids*ctx->batch_sz, NULL, globXArray);CHKERRQ(ierr); 444 if (stepi > rectx->plotStep && rectx->plotting) { 445 rectx->plotting = PETSC_FALSE; /* was doing diagnostics, now done */ 446 rectx->plotIdx++; 447 } 448 /* view */ 449 ierr = TSGetConvergedReason(ts,&reason);CHKERRQ(ierr); 450 if (time/rectx->plotDt >= (PetscReal)rectx->plotIdx || reason) { 451 if ((reason || stepi==0 || rectx->plotIdx%rectx->print_period==0) && ctx->verbose > 0) { 452 /* print norms */ 453 ierr = LandauPrintNorms(X, stepi);CHKERRQ(ierr); 454 } 455 if (!rectx->plotting) { /* first step of possible backtracks */ 456 rectx->plotting = PETSC_TRUE; 457 /* diagnostics + change E field with Sptizer (not just a monitor) */ 458 ierr = rectx->test(ts,X,stepi,time,reason ? PETSC_TRUE : PETSC_FALSE, ctx, rectx);CHKERRQ(ierr); 459 } else { 460 PetscPrintf(PETSC_COMM_WORLD, "\t\t ERROR SKIP test spit ------\n"); 461 rectx->plotting = PETSC_TRUE; 462 } 463 ierr = PetscObjectSetName((PetscObject) globXArray[ LAND_PACK_IDX(ctx->batch_view_idx,rectx->grid_view_idx) ], rectx->grid_view_idx==0 ? "ue" : "ui");CHKERRQ(ierr); 464 /* view, overwrite step when back tracked */ 465 ierr = DMSetOutputSequenceNumber(pack, rectx->plotIdx, time*ctx->t_0);CHKERRQ(ierr); 466 ierr = VecViewFromOptions(globXArray[ LAND_PACK_IDX(ctx->batch_view_idx, rectx->grid_view_idx) ],NULL,"-vec_view");CHKERRQ(ierr); 467 468 rectx->plotStep = stepi; 469 } else { 470 if (rectx->plotting) PetscPrintf(PETSC_COMM_WORLD," ERROR rectx->plotting=%D step %D\n",rectx->plotting,stepi); 471 /* diagnostics + change E field with Sptizer (not just a monitor) - can we lag this? */ 472 ierr = rectx->test(ts,X,stepi,time,reason ? PETSC_TRUE : PETSC_FALSE, ctx, rectx);CHKERRQ(ierr); 473 } 474 /* parallel check that only works of all batches are identical */ 475 if (reason && ctx->verbose > 3) { 476 PetscReal val,rval; 477 PetscMPIInt rank; 478 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank);CHKERRMPI(ierr); 479 for (PetscInt grid=0;grid<ctx->num_grids;grid++) { 480 PetscInt nerrors=0; 481 for (PetscInt i=0; i<ctx->batch_sz;i++) { 482 ierr = VecNorm(globXArray[ LAND_PACK_IDX(i,grid) ],NORM_2,&val);CHKERRQ(ierr); 483 if (i==0) rval = val; 484 else if ((val=PetscAbs(val-rval)/rval) > 1000*PETSC_MACHINE_EPSILON) { 485 PetscPrintf(PETSC_COMM_SELF, " [%D] Warning %D.%D) diff = %2.15e\n",rank,grid,i,val);CHKERRQ(ierr); 486 nerrors++; 487 } 488 } 489 if (nerrors) { 490 ierr = PetscPrintf(PETSC_COMM_SELF, " ***** [%D] ERROR max %D errors\n",rank,nerrors);CHKERRQ(ierr); 491 } else { 492 ierr = PetscPrintf(PETSC_COMM_WORLD, "[%D] %D) batch consistency check OK\n",rank,grid);CHKERRQ(ierr); 493 } 494 } 495 } 496 rectx->idx = 0; 497 ierr = DMCompositeRestoreAccessArray(pack, X, ctx->num_grids*ctx->batch_sz, NULL, globXArray);CHKERRQ(ierr); 498 PetscFunctionReturn(0); 499 } 500 501 PetscErrorCode PreStep(TS ts) 502 { 503 PetscErrorCode ierr; 504 LandauCtx *ctx; 505 REctx *rectx; 506 DM dm; 507 PetscInt stepi; 508 PetscReal time; 509 Vec X; 510 511 PetscFunctionBeginUser; 512 /* not used */ 513 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 514 ierr = TSGetTime(ts,&time);CHKERRQ(ierr); 515 ierr = TSGetSolution(ts,&X);CHKERRQ(ierr); 516 ierr = DMGetApplicationContext(dm, &ctx);CHKERRQ(ierr); 517 rectx = (REctx*)ctx->data; 518 ierr = TSGetStepNumber(ts, &stepi);CHKERRQ(ierr); 519 /* update E */ 520 ierr = rectx->E(X, NULL, stepi, time, ctx, &ctx->Ez);CHKERRQ(ierr); 521 PetscFunctionReturn(0); 522 } 523 524 /* 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) */ 525 static PetscErrorCode stepSrc(PetscReal time, PetscReal *rho, LandauCtx *ctx) 526 { 527 REctx *rectx = (REctx*)ctx->data; 528 529 PetscFunctionBeginUser; 530 if (time >= rectx->pulse_start) *rho = rectx->pulse_rate; 531 else *rho = 0.; 532 PetscFunctionReturn(0); 533 } 534 static PetscErrorCode zeroSrc(PetscReal time, PetscReal *rho, LandauCtx *ctx) 535 { 536 PetscFunctionBeginUser; 537 *rho = 0.; 538 PetscFunctionReturn(0); 539 } 540 static PetscErrorCode pulseSrc(PetscReal time, PetscReal *rho, LandauCtx *ctx) 541 { 542 REctx *rectx = (REctx*)ctx->data; 543 544 PetscFunctionBeginUser; 545 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'"); 546 if (time < rectx->pulse_start || time > rectx->pulse_start + 3*rectx->pulse_width) *rho = 0; 547 /* else if (0) { */ 548 /* 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); */ 549 /* *rho = rectx->pulse_rate * (cycle / (stop - start)) / (1 + PetscExpReal(steep*(PetscSinReal(2*PETSC_PI*((t - start)/cycle + xi)) - PetscSinReal(2*PETSC_PI*xi)))); */ 550 /* } else if (0) { */ 551 /* double x = 2*(time - rectx->pulse_start)/(3*rectx->pulse_width) - 1; */ 552 /* if (x==1 || x==-1) *rho = 0; */ 553 /* else *rho = rectx->pulse_rate * PetscExpReal(-1/(1-x*x)); */ 554 /* } */ 555 else { 556 double x = PetscSinReal((time-rectx->pulse_start)/(3*rectx->pulse_width)*2*PETSC_PI - PETSC_PI/2) + 1; /* 0:2, integrates to 1.0 */ 557 *rho = rectx->pulse_rate * x / (3*rectx->pulse_width); 558 if (!rectx->use_spitzer_eta) rectx->use_spitzer_eta = PETSC_TRUE; /* use it next time */ 559 } 560 PetscFunctionReturn(0); 561 } 562 563 #undef __FUNCT__ 564 #define __FUNCT__ "ProcessREOptions" 565 static PetscErrorCode ProcessREOptions(REctx *rectx, const LandauCtx *ctx, DM dm, const char prefix[]) 566 { 567 PetscErrorCode ierr; 568 PetscFunctionList plist = NULL, testlist = NULL, elist = NULL; 569 char pname[256],testname[256],ename[256]; 570 DM dm_dummy; 571 PetscBool Connor_E = PETSC_FALSE; 572 573 PetscFunctionBeginUser; 574 ierr = DMCreate(PETSC_COMM_WORLD,&dm_dummy);CHKERRQ(ierr); 575 rectx->Ne_ion = 1; /* number of electrons given up by impurity ion */ 576 rectx->T_cold = .005; /* kev */ 577 rectx->ion_potential = 15; /* ev */ 578 rectx->L = 2; 579 rectx->X_0 = NULL; 580 rectx->imp_idx = ctx->num_species - 1; /* default ionized impurity as last one */ 581 rectx->pulse_start = PETSC_MAX_REAL; 582 rectx->pulse_width = 1; 583 rectx->plotStep = PETSC_MAX_INT; 584 rectx->pulse_rate = 1.e-1; 585 rectx->current_rate = 0; 586 rectx->plotIdx = 0; 587 rectx->j = 0; 588 rectx->plotDt = 1.0; 589 rectx->plotting = PETSC_FALSE; 590 rectx->use_spitzer_eta = PETSC_FALSE; 591 rectx->idx = 0; 592 rectx->print_period = 10; 593 rectx->grid_view_idx = 0; 594 /* Register the available impurity sources */ 595 ierr = PetscFunctionListAdd(&plist,"step",&stepSrc);CHKERRQ(ierr); 596 ierr = PetscFunctionListAdd(&plist,"none",&zeroSrc);CHKERRQ(ierr); 597 ierr = PetscFunctionListAdd(&plist,"pulse",&pulseSrc);CHKERRQ(ierr); 598 ierr = PetscStrcpy(pname,"none");CHKERRQ(ierr); 599 ierr = PetscFunctionListAdd(&testlist,"none",&testNone);CHKERRQ(ierr); 600 ierr = PetscFunctionListAdd(&testlist,"spitzer",&testSpitzer);CHKERRQ(ierr); 601 ierr = PetscFunctionListAdd(&testlist,"stable",&testStable);CHKERRQ(ierr); 602 ierr = PetscStrcpy(testname,"none");CHKERRQ(ierr); 603 ierr = PetscFunctionListAdd(&elist,"none",&ENone);CHKERRQ(ierr); 604 ierr = PetscFunctionListAdd(&elist,"induction",&EInduction);CHKERRQ(ierr); 605 ierr = PetscFunctionListAdd(&elist,"constant",&EConstant);CHKERRQ(ierr); 606 ierr = PetscStrcpy(ename,"constant");CHKERRQ(ierr); 607 608 ierr = PetscOptionsBegin(PETSC_COMM_SELF, prefix, "Options for Runaway/seed electron model", "none");CHKERRQ(ierr); 609 ierr = PetscOptionsReal("-ex2_plot_dt", "Plotting interval", "ex2.c", rectx->plotDt, &rectx->plotDt, NULL);CHKERRQ(ierr); 610 if (rectx->plotDt < 0) rectx->plotDt = 1e30; 611 if (rectx->plotDt == 0) rectx->plotDt = 1e-30; 612 ierr = PetscOptionsInt("-ex2_print_period", "Plotting interval", "ex2.c", rectx->print_period, &rectx->print_period, NULL);CHKERRQ(ierr); 613 ierr = PetscOptionsInt("-ex2_grid_view_idx", "grid_view_idx", "ex2.c", rectx->grid_view_idx, &rectx->grid_view_idx, NULL);CHKERRQ(ierr); 614 if (rectx->grid_view_idx >= ctx->num_grids) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"rectx->grid_view_idx (%D) >= ctx->num_grids (%D)",rectx->imp_idx,ctx->num_grids); 615 ierr = PetscOptionsFList("-ex2_impurity_source_type","Name of impurity source to run","",plist,pname,pname,sizeof(pname),NULL);CHKERRQ(ierr); 616 ierr = PetscOptionsFList("-ex2_test_type","Name of test to run","",testlist,testname,testname,sizeof(testname),NULL);CHKERRQ(ierr); 617 ierr = PetscOptionsInt("-ex2_impurity_index", "index of sink for impurities", "none", rectx->imp_idx, &rectx->imp_idx, NULL);CHKERRQ(ierr); 618 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); 619 ierr = PetscOptionsFList("-ex2_e_field_type","Electric field type","",elist,ename,ename,sizeof(ename),NULL);CHKERRQ(ierr); 620 rectx->Ne_ion = -ctx->charges[rectx->imp_idx]/ctx->charges[0]; 621 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); 622 ierr = PetscOptionsReal("-ex2_pulse_start_time","Time at which pulse happens for 'pulse' source","none",rectx->pulse_start,&rectx->pulse_start, NULL);CHKERRQ(ierr); 623 ierr = PetscOptionsReal("-ex2_pulse_width_time","Width of pulse 'pulse' source","none",rectx->pulse_width,&rectx->pulse_width, NULL);CHKERRQ(ierr); 624 ierr = PetscOptionsReal("-ex2_pulse_rate","Number density of pulse for 'pulse' source","none",rectx->pulse_rate,&rectx->pulse_rate, NULL);CHKERRQ(ierr); 625 rectx->T_cold *= 1.16e7; /* convert to Kelvin */ 626 ierr = PetscOptionsReal("-ex2_ion_potential","Potential to ionize impurity (should be array) in ev","none",rectx->ion_potential,&rectx->ion_potential, NULL);CHKERRQ(ierr); 627 ierr = PetscOptionsReal("-ex2_inductance","Inductance E feild","none",rectx->L,&rectx->L, NULL);CHKERRQ(ierr); 628 ierr = PetscOptionsBool("-ex2_connor_e_field_units","Scale Ex but Connor-Hastie E_c","none",Connor_E,&Connor_E, NULL);CHKERRQ(ierr); 629 ierr = 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);CHKERRQ(ierr); 630 ierr = PetscOptionsEnd();CHKERRQ(ierr); 631 /* get impurity source rate function */ 632 ierr = PetscFunctionListFind(plist,pname,&rectx->impuritySrcRate);CHKERRQ(ierr); 633 if (!rectx->impuritySrcRate) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"No impurity source function found '%s'",pname); 634 ierr = PetscFunctionListFind(testlist,testname,&rectx->test);CHKERRQ(ierr); 635 if (!rectx->test) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"No test found '%s'",testname); 636 ierr = PetscFunctionListFind(elist,ename,&rectx->E);CHKERRQ(ierr); 637 if (!rectx->E) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"No E field function found '%s'",ename); 638 ierr = PetscFunctionListDestroy(&plist);CHKERRQ(ierr); 639 ierr = PetscFunctionListDestroy(&testlist);CHKERRQ(ierr); 640 ierr = PetscFunctionListDestroy(&elist);CHKERRQ(ierr); 641 642 /* convert E from Connor-Hastie E_c units to real if doing Spitzer E */ 643 if (Connor_E) { 644 PetscReal E = ctx->Ez, Tev = ctx->thermal_temps[0]*8.621738e-5, n = ctx->n_0*ctx->n[0]; 645 CalculateE(Tev, n, ctx->lnLam, ctx->epsilon0, &E); 646 ((LandauCtx*)ctx)->Ez *= E; 647 } 648 ierr = DMDestroy(&dm_dummy);CHKERRQ(ierr); 649 PetscFunctionReturn(0); 650 } 651 652 static PetscErrorCode MatrixNfDestroy(void *ptr) 653 { 654 PetscInt *nf = (PetscInt *)ptr; 655 PetscErrorCode ierr; 656 PetscFunctionBegin; 657 ierr = PetscFree(nf);CHKERRQ(ierr); 658 PetscFunctionReturn(0); 659 } 660 661 int main(int argc, char **argv) 662 { 663 DM pack; 664 Vec X,*XsubArray; 665 PetscErrorCode ierr; 666 PetscInt dim = 2, nDMs; 667 TS ts; 668 Mat J; 669 PetscDS prob; 670 LandauCtx *ctx; 671 REctx *rectx; 672 #if defined PETSC_USE_LOG 673 PetscLogStage stage; 674 #endif 675 PetscMPIInt rank; 676 #if defined(PETSC_HAVE_THREADSAFETY) 677 double starttime, endtime; 678 #endif 679 PetscContainer container; 680 ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr; 681 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank);CHKERRMPI(ierr); 682 if (rank) { /* turn off output stuff for duplicate runs */ 683 ierr = PetscOptionsClearValue(NULL,"-dm_view");CHKERRQ(ierr); 684 ierr = PetscOptionsClearValue(NULL,"-vec_view");CHKERRQ(ierr); 685 ierr = PetscOptionsClearValue(NULL,"-dm_view_diff");CHKERRQ(ierr); 686 ierr = PetscOptionsClearValue(NULL,"-vec_view_diff");CHKERRQ(ierr); 687 ierr = PetscOptionsClearValue(NULL,"-dm_view_sources");CHKERRQ(ierr); 688 ierr = PetscOptionsClearValue(NULL,"-vec_view_0");CHKERRQ(ierr); 689 ierr = PetscOptionsClearValue(NULL,"-dm_view_0");CHKERRQ(ierr); 690 ierr = PetscOptionsClearValue(NULL,"-vec_view_sources");CHKERRQ(ierr); 691 ierr = PetscOptionsClearValue(NULL,"-info");CHKERRQ(ierr); /* this does not work */ 692 } 693 ierr = PetscOptionsGetInt(NULL,NULL, "-dim", &dim, NULL);CHKERRQ(ierr); 694 /* Create a mesh */ 695 ierr = LandauCreateVelocitySpace(PETSC_COMM_WORLD, dim, "", &X, &J, &pack);CHKERRQ(ierr); 696 ierr = DMCompositeGetNumberDM(pack,&nDMs);CHKERRQ(ierr); 697 ierr = PetscMalloc(sizeof(*XsubArray)*nDMs, &XsubArray);CHKERRQ(ierr); 698 ierr = PetscObjectSetName((PetscObject)J, "Jacobian");CHKERRQ(ierr); 699 ierr = PetscObjectSetName((PetscObject)X, "f");CHKERRQ(ierr); 700 ierr = LandauCreateMassMatrix(pack, NULL);CHKERRQ(ierr); 701 ierr = DMGetApplicationContext(pack, &ctx);CHKERRQ(ierr); 702 ierr = DMSetUp(pack);CHKERRQ(ierr); 703 /* context */ 704 ierr = PetscNew(&rectx);CHKERRQ(ierr); 705 ctx->data = rectx; 706 ierr = ProcessREOptions(rectx,ctx,pack,"");CHKERRQ(ierr); 707 ierr = DMGetDS(pack, &prob);CHKERRQ(ierr); 708 ierr = DMCompositeGetAccessArray(pack, X, nDMs, NULL, XsubArray);CHKERRQ(ierr); // read only 709 ierr = PetscObjectSetName((PetscObject) XsubArray[ LAND_PACK_IDX(ctx->batch_view_idx, rectx->grid_view_idx) ], rectx->grid_view_idx==0 ? "ue" : "ui");CHKERRQ(ierr); 710 ierr = DMViewFromOptions(ctx->plex[rectx->grid_view_idx],NULL,"-dm_view");CHKERRQ(ierr); 711 ierr = DMViewFromOptions(ctx->plex[rectx->grid_view_idx], NULL,"-dm_view_0");CHKERRQ(ierr); 712 ierr = VecViewFromOptions(XsubArray[ LAND_PACK_IDX(ctx->batch_view_idx,rectx->grid_view_idx) ], NULL,"-vec_view_0");CHKERRQ(ierr); // initial condition (monitor plots after step) 713 ierr = DMCompositeRestoreAccessArray(pack, X, nDMs, NULL, XsubArray);CHKERRQ(ierr); // read only 714 ierr = PetscFree(XsubArray);CHKERRQ(ierr); 715 ierr = VecViewFromOptions(X, NULL,"-vec_view_global");CHKERRQ(ierr); // initial condition (monitor plots after step) 716 ierr = DMSetOutputSequenceNumber(pack, 0, 0.0);CHKERRQ(ierr); 717 /* Create timestepping solver context */ 718 ierr = TSCreate(PETSC_COMM_SELF,&ts);CHKERRQ(ierr); 719 ierr = TSSetDM(ts,pack);CHKERRQ(ierr); 720 ierr = TSSetIFunction(ts,NULL,LandauIFunction,NULL);CHKERRQ(ierr); 721 ierr = TSSetIJacobian(ts,J,J,LandauIJacobian,NULL);CHKERRQ(ierr); 722 ierr = TSSetRHSFunction(ts,NULL,FormSource,NULL);CHKERRQ(ierr); 723 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 724 ierr = TSSetSolution(ts,X);CHKERRQ(ierr); 725 ierr = TSSetApplicationContext(ts, ctx);CHKERRQ(ierr); 726 ierr = TSMonitorSet(ts,Monitor,ctx,NULL);CHKERRQ(ierr); 727 ierr = TSSetPreStep(ts,PreStep);CHKERRQ(ierr); 728 { 729 PetscInt *pNf; 730 ierr = PetscContainerCreate(PETSC_COMM_SELF, &container);CHKERRQ(ierr); 731 ierr = PetscMalloc(sizeof(*pNf), &pNf);CHKERRQ(ierr); 732 *pNf = ctx->batch_sz + 100000*ctx->numConcurrency; 733 ierr = PetscContainerSetPointer(container, (void *)pNf);CHKERRQ(ierr); 734 ierr = PetscContainerSetUserDestroy(container, MatrixNfDestroy);CHKERRQ(ierr); 735 ierr = PetscObjectCompose((PetscObject)ctx->J, "batch size", (PetscObject) container);CHKERRQ(ierr); 736 ierr = PetscContainerDestroy(&container);CHKERRQ(ierr); 737 } 738 rectx->Ez_initial = ctx->Ez; /* cache for induction caclulation - applied E field */ 739 if (1) { /* warm up an test just LandauIJacobian */ 740 Vec vec; 741 PetscInt nsteps; 742 PetscReal dt; 743 ierr = PetscLogStageRegister("Warmup", &stage);CHKERRQ(ierr); 744 ierr = PetscLogStagePush(stage);CHKERRQ(ierr); 745 ierr = VecDuplicate(X,&vec);CHKERRQ(ierr); 746 ierr = VecCopy(X,vec);CHKERRQ(ierr); 747 ierr = TSGetMaxSteps(ts,&nsteps);CHKERRQ(ierr); 748 ierr = TSGetTimeStep(ts,&dt);CHKERRQ(ierr); 749 ierr = TSSetMaxSteps(ts,1);CHKERRQ(ierr); 750 ierr = TSSolve(ts,X);CHKERRQ(ierr); 751 ierr = TSSetMaxSteps(ts,nsteps);CHKERRQ(ierr); 752 ierr = TSSetStepNumber(ts,0);CHKERRQ(ierr); 753 ierr = TSSetTime(ts,0);CHKERRQ(ierr); 754 ierr = TSSetTimeStep(ts,dt);CHKERRQ(ierr); 755 rectx->plotIdx = 0; 756 rectx->plotting = PETSC_FALSE; 757 ierr = PetscLogStagePop();CHKERRQ(ierr); 758 ierr = VecCopy(vec,X);CHKERRQ(ierr); 759 ierr = VecDestroy(&vec);CHKERRQ(ierr); 760 ctx->aux_bool = PETSC_FALSE; // flag for not a clean Jacobian 761 } 762 /* go */ 763 ierr = PetscLogStageRegister("Solve", &stage);CHKERRQ(ierr); 764 ctx->stage = 0; // lets not use this stage 765 #if defined(PETSC_HAVE_THREADSAFETY) 766 ctx->stage = 1; // not set with thread safty 767 #endif 768 ierr = TSSetSolution(ts,X);CHKERRQ(ierr); 769 ierr = PetscLogStagePush(stage);CHKERRQ(ierr); 770 #if defined(PETSC_HAVE_THREADSAFETY) 771 starttime = MPI_Wtime(); 772 #endif 773 ierr = TSSolve(ts,X);CHKERRQ(ierr); 774 ierr = PetscLogStagePop();CHKERRQ(ierr); 775 #if defined(PETSC_HAVE_THREADSAFETY) 776 endtime = MPI_Wtime(); 777 ctx->times[LANDAU_EX2_TSSOLVE] += (endtime - starttime); 778 #endif 779 ierr = VecViewFromOptions(X, NULL,"-vec_view_global");CHKERRQ(ierr); 780 /* clean up */ 781 ierr = LandauDestroyVelocitySpace(&pack);CHKERRQ(ierr); 782 ierr = TSDestroy(&ts);CHKERRQ(ierr); 783 ierr = VecDestroy(&X);CHKERRQ(ierr); 784 ierr = PetscFree(rectx);CHKERRQ(ierr); 785 ierr = PetscContainerDestroy(&container);CHKERRQ(ierr); 786 ierr = PetscFinalize(); 787 return ierr; 788 } 789 790 /*TEST 791 792 testset: 793 requires: p4est !complex double 794 output_file: output/ex2_0.out 795 args: -dm_landau_num_species_grid 1,1 -dm_landau_Ez 0 -petscspace_degree 3 -petscspace_poly_tensor 1 -dm_landau_type p4est -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 2,2 -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_landau_gpu_assembly true -dm_landau_batch_size 2 796 test: 797 suffix: cpu 798 args: -dm_landau_device_type cpu 799 test: 800 suffix: kokkos 801 requires: kokkos_kernels 802 args: -dm_landau_device_type kokkos -dm_mat_type aijkokkos -dm_vec_type kokkos 803 test: 804 suffix: cuda 805 requires: cuda 806 args: -dm_landau_device_type cuda -dm_mat_type aijcusparse -dm_vec_type cuda -mat_cusparse_use_cpu_solve 807 808 TEST*/ 809