Actual source code: plexland.c
petsc-3.14.6 2021-03-30
1: #include <petsc/private/dmpleximpl.h>
2: #include <petsclandau.h>
3: #include <petscts.h>
4: #include <petscdmforest.h>
6: /* Landau collision operator */
7: #define PETSC_THREAD_SYNC
8: #define PETSC_DEVICE_FUNC_DECL static
9: #include "land_kernel.h"
11: #define LANDAU_VL 1
12: static PetscErrorCode LandauPointDataCreate(PetscReal **IPData, PetscInt dim, PetscInt nip, PetscInt Ns)
13: {
14: PetscErrorCode ierr, d, s, jj, nip_pad = LANDAU_VL*(nip/LANDAU_VL + !!(nip%LANDAU_VL)), pnt_sz = (dim + Ns*(1+dim));
15: PetscReal *pdata;
17: PetscMalloc(nip_pad*pnt_sz*sizeof(PetscReal),IPData);
18: /* pad with zeros in case we vectorize into this */
19: for (jj=nip, pdata = *IPData + nip*pnt_sz; jj < nip_pad; jj++, pdata += pnt_sz){
20: LandauPointData *fplpt = (LandauPointData*)pdata; /* [dim + NS*(1+dim)] */
21: for (d=0;d<dim;d++) fplpt->crd[d] = -1;
22: for (s=0;s<Ns;s++) {
23: fplpt->fdf[s].f = 0;
24: for (d=0;d<dim;d++) fplpt->fdf[s].df[d] = 0;
25: }
26: }
27: return(0);
28: }
30: static PetscErrorCode LandauPointDataDestroy(PetscReal *IPData)
31: {
32: PetscErrorCode ierr;
34: PetscFree(IPData);
35: return(0);
36: }
37: /* ------------------------------------------------------------------- */
38: /*
39: LandauFormJacobian_Internal - Evaluates Jacobian matrix.
41: Input Parameters:
42: . globX - input vector
43: . actx - optional user-defined context
44: . dim - dimension
46: Output Parameters:
47: . J0acP - Jacobian matrix filled, not created
48: */
49: PetscErrorCode LandauFormJacobian_Internal(Vec a_X, Mat JacP, const PetscInt dim, void *a_ctx)
50: {
51: LandauCtx *ctx = (LandauCtx*)a_ctx;
52: PetscErrorCode ierr;
53: PetscInt cStart, cEnd, elemMatSize;
54: DM plex = NULL;
55: PetscDS prob;
56: PetscSection section,globsection;
57: PetscScalar *elemMat;
58: PetscInt numCells,totDim,ej,Nq,*Nbf,*Ncf,Nb,Ncx,Nf,d,f,fieldA,Nip,nip_pad,ipdata_sz;
59: PetscQuadrature quad;
60: PetscTabulation *Tf;
61: PetscReal *wiGlob, nu_alpha[LANDAU_MAX_SPECIES], nu_beta[LANDAU_MAX_SPECIES];
62: const PetscReal *quadWeights;
63: PetscReal *IPData,*invJ,*invJ_a;
64: PetscReal invMass[LANDAU_MAX_SPECIES],Eq_m[LANDAU_MAX_SPECIES],m_0=ctx->m_0; /* normalize mass -- not needed! */
65: PetscLogDouble flops;
66: Vec locX;
73: PetscLogEventBegin(ctx->events[1],0,0,0,0);
74: DMConvert(ctx->dmv, DMPLEX, &plex);
75: DMCreateLocalVector(plex, &locX);
76: VecZeroEntries(locX); /* zero BCs so don't set */
77: DMGlobalToLocalBegin(plex, a_X, INSERT_VALUES, locX);
78: DMGlobalToLocalEnd (plex, a_X, INSERT_VALUES, locX);
79: DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd);
80: DMGetLocalSection(plex, §ion);
81: DMGetGlobalSection(plex, &globsection);
82: DMGetDS(plex, &prob);
83: PetscDSGetTabulation(prob, &Tf); // Bf, &Df
84: PetscDSGetDimensions(prob, &Nbf); Nb = Nbf[0]; /* number of vertices*S */
85: PetscSectionGetNumFields(section, &Nf); if (Nf!=ctx->num_species) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Nf %D != S",Nf);
86: PetscDSGetComponents(prob, &Ncf); Ncx = Ncf[0]; if (Ncx!=1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Nc %D != 1",Ncx);
87: for (fieldA=0;fieldA<Nf;fieldA++) {
88: invMass[fieldA] = m_0/ctx->masses[fieldA];
89: Eq_m[fieldA] = -ctx->Ez * ctx->t_0 * ctx->charges[fieldA] / (ctx->v_0 * ctx->masses[fieldA]); /* normalize dimensionless */
90: if (dim==2) Eq_m[fieldA] *= 2 * PETSC_PI; /* add the 2pi term that is not in Landau */
91: nu_alpha[fieldA] = PetscSqr(ctx->charges[fieldA]/m_0)*m_0/ctx->masses[fieldA];
92: nu_beta[fieldA] = PetscSqr(ctx->charges[fieldA]/ctx->epsilon0)*ctx->lnLam / (8*PETSC_PI) * ctx->t_0*ctx->n_0/PetscPowReal(ctx->v_0,3);
93: }
94: PetscDSGetTotalDimension(prob, &totDim);
95: numCells = cEnd - cStart;
96: PetscFEGetQuadrature(ctx->fe[0], &quad);
97: PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, &quadWeights);
98: if (Nb!=Nq) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Nb!=Nq %D %D over integration or simplices? Tf[0]->Nb=%D dim=%D",Nb,Nq,Tf[0]->Nb,dim);
99: if (Nq >LANDAU_MAX_NQ) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Order too high. Nq = %D > LANDAU_MAX_NQ (%D)",Nq,LANDAU_MAX_NQ);
100: Nip = numCells*Nq;
101: nip_pad = LANDAU_VL*(Nip/LANDAU_VL + !!(Nip%LANDAU_VL));
102: flops = (PetscLogDouble)numCells*(PetscLogDouble)Nq*(PetscLogDouble)(5*dim*dim*Nf*Nf + 165);
103: MatZeroEntries(JacP);
104: elemMatSize = totDim*totDim;
105: {
106: static int cc = 0;
107: PetscScalar uu[LANDAU_MAX_SPECIES],u_x[LANDAU_MAX_SPECIES*LANDAU_DIM];
108: /* collect f data */
109: if (ctx->verbose > 2 || (ctx->verbose > 0 && cc++ == 0)) {
110: PetscInt N,Nloc;
111: MatGetSize(JacP,&N,NULL);
112: VecGetSize(locX,&Nloc);
113: PetscPrintf(PETSC_COMM_WORLD,"[%D]%s: %D IPs, %D cells, %s elements, totDim=%D, Nb=%D, Nq=%D, elemMatSize=%D, dim=%D, Tab: Nb=%D Nf=%D Np=%D cdim=%D N=%D N_loc=%D\n",
114: 0,"FormLandau",Nq*numCells,numCells,ctx->simplex ? "SIMPLEX" : "TENSOR", totDim, Nb, Nq, elemMatSize, dim, Tf[0]->Nb, Nf, Tf[0]->Np, Tf[0]->cdim, N, Nloc);
115: }
116: LandauPointDataCreate(&IPData, dim, Nq*numCells, Nf);
117: ipdata_sz = (dim + Nf*(1+dim));
118: PetscMalloc3(elemMatSize,&elemMat,nip_pad,&wiGlob,nip_pad*dim*dim,&invJ_a);
119: /* cache geometry and x, f and df/dx at IPs */
120: for (ej = 0, invJ = invJ_a ; ej < numCells; ++ej, invJ += Nq*dim*dim) {
121: PetscReal vj[LANDAU_MAX_NQ*LANDAU_DIM],detJj[LANDAU_MAX_NQ], Jdummy[LANDAU_MAX_NQ*LANDAU_DIM*LANDAU_DIM];
122: PetscInt qj,f;
123: PetscScalar *coef = NULL;
124: DMPlexComputeCellGeometryFEM(plex, cStart+ej, quad, vj, Jdummy, invJ, detJj);
125: DMPlexVecGetClosure(plex, section, locX, cStart+ej, NULL, &coef);
126: /* create point data for cell i for Landau tensor: x, f(x), grad f(x) */
127: for (qj = 0; qj < Nq; ++qj) {
128: PetscInt gidx = (ej*Nq + qj);
129: LandauPointData *pnt_data = (LandauPointData*)(IPData + gidx*ipdata_sz);
130: PetscScalar refSpaceDer[LANDAU_DIM];
131: PetscInt dOffset = 0, fOffset = 0;
132: for (d = 0; d < dim; ++d) pnt_data->crd[d] = vj[qj * dim + d]; /* coordinate */
133: wiGlob[gidx] = detJj[qj] * quadWeights[qj];
134: if (dim==2) wiGlob[gidx] *= pnt_data->crd[0]; /* cylindrical coordinate, w/o 2pi */
135: /* get u & du (EvaluateFieldJets) */
136: for (f = 0; f < Nf; ++f) {
137: const PetscReal *Bq = &Tf[f]->T[0][qj*Nb];
138: const PetscReal *Dq = &Tf[f]->T[1][qj*Nb*dim];
139: PetscInt b, e;
140: uu[fOffset] = 0.0;
141: for (d = 0; d < LANDAU_DIM; ++d) refSpaceDer[d] = 0.0;
142: for (b = 0; b < Nb; ++b) {
143: const PetscInt cidx = b;
144: uu[fOffset] += Bq[cidx]*coef[dOffset+cidx];
145: for (d = 0; d < dim; ++d) refSpaceDer[d] += Dq[cidx*dim+d]*coef[dOffset+cidx];
146: }
147: for (d = 0; d < dim; ++d) {
148: for (e = 0, u_x[fOffset*dim+d] = 0.0; e < dim; ++e) { // should add directly to point data here!!!
149: u_x[fOffset*dim+d] += invJ[qj * dim * dim + e*dim+d]*refSpaceDer[e];
150: }
151: }
152: fOffset += 1;
153: dOffset += Nb;
154: }
155: /* copy to IPDataLocal */
156: for (f=0;f<Nf;f++) {
157: pnt_data->fdf[f].f = PetscRealPart(uu[f]);
158: for (d = 0; d < dim; ++d) pnt_data->fdf[f].df[d] = PetscRealPart(u_x[f*dim+d]);
159: }
160: } /* q */
161: DMPlexVecRestoreClosure(plex, section, locX, cStart+ej, NULL, &coef);
162: } /* e */
163: }
164: DMRestoreLocalVector(plex, &locX);
165: PetscLogEventEnd(ctx->events[1],0,0,0,0);
167: /* outer element loop j is like a regular assembly loop */
168: if (ctx->deviceType == LANDAU_CUDA) {
169: #if defined(PETSC_HAVE_CUDA)
170: LandauCUDAJacobian(plex,Nq,nu_alpha,nu_beta,invMass,Eq_m,IPData,wiGlob,invJ_a,ctx->subThreadBlockSize,ctx->events,ctx->quarter3DDomain,JacP);
171: #else
172: SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"-landau_device_type %s not built","cuda");
173: #endif
174: } else if (ctx->deviceType == LANDAU_KOKKOS) {
175: #if defined(PETSC_HAVE_KOKKOS)
176: LandauKokkosJacobian(plex,Nq,nu_alpha,nu_beta,invMass,Eq_m,IPData,wiGlob,invJ_a,ctx->subThreadBlockSize,ctx->events,ctx->quarter3DDomain,JacP);
177: #else
178: SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"-landau_device_type %s not built","kokkos");
179: #endif
180: } else { /* CPU version */
181: for (ej = cStart, invJ = invJ_a; ej < cEnd; ++ej, invJ += Nq*dim*dim) {
182: PetscInt qj;
183: PetscLogEventBegin(ctx->events[3],0,0,0,0);
184: PetscMemzero(elemMat, totDim *totDim * sizeof(PetscScalar));
185: PetscLogEventEnd(ctx->events[3],0,0,0,0);
186: for (qj = 0; qj < Nq; ++qj) {
187: PetscReal g2[1][LANDAU_MAX_SUB_THREAD_BLOCKS][LANDAU_MAX_SPECIES][LANDAU_DIM], g3[1][LANDAU_MAX_SUB_THREAD_BLOCKS][LANDAU_MAX_SPECIES][LANDAU_DIM][LANDAU_DIM];
188: const PetscInt nip = numCells*Nq, jpidx = Nq*(ej-cStart) + qj, one = 1, zero = 0; /* length of inner global interation, outer integration point */
190: PetscLogEventBegin(ctx->events[4],0,0,0,0);
191: PetscLogFlops(flops);
192: landau_inner_integral(zero, one, zero, one, zero, nip, 1, jpidx, Nf, dim, IPData, wiGlob, &invJ[qj*dim*dim], nu_alpha, nu_beta, invMass, Eq_m, ctx->quarter3DDomain, Nq, Nb, qj, qj+1, Tf[0]->T[0], Tf[0]->T[1], elemMat, g2, g3, ej);
193: PetscLogEventEnd(ctx->events[4],0,0,0,0);
194: } /* qj loop */
195: /* assemble matrix */
196: PetscLogEventBegin(ctx->events[6],0,0,0,0);
197: DMPlexMatSetClosure(plex, section, globsection, JacP, ej, elemMat, ADD_VALUES);
198: PetscLogEventEnd(ctx->events[6],0,0,0,0);
200: if (ej==-1) {
201: PetscPrintf(PETSC_COMM_SELF, "CPU Element matrix\n");
202: for (d = 0; d < totDim; ++d){
203: for (f = 0; f < totDim; ++f) PetscPrintf(PETSC_COMM_SELF," %17.9e", PetscRealPart(elemMat[d*totDim + f]));
204: PetscPrintf(PETSC_COMM_SELF,"\n");
205: }
206: }
207: } /* ej cells loop, not cuda */
208: }
209: // PetscSleep(2); exit(13);
211: /* assemble matrix or vector */
212: PetscLogEventBegin(ctx->events[7],0,0,0,0);
213: MatAssemblyBegin(JacP, MAT_FINAL_ASSEMBLY);
214: MatAssemblyEnd(JacP, MAT_FINAL_ASSEMBLY);
215: MatScale(JacP, -1.0); /* The code reflect the papers: du/dt = C, whereas PETSc use the form G(u) = du/dt - C(u) = 0 */
216: PetscLogEventEnd(ctx->events[7],0,0,0,0);
218: /* clean up */
219: PetscFree3(elemMat,wiGlob,invJ_a);
220: DMDestroy(&plex);
221: LandauPointDataDestroy(IPData);
222: return(0);
223: }
225: #if defined(LANDAU_ADD_BCS)
226: static void zero_bc(PetscInt dim, PetscInt Nf, PetscInt NfAux,
227: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
228: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
229: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar uexact[])
230: {
231: uexact[0] = 0;
232: }
233: #endif
235: #define MATVEC2(__a,__x,__p) {int i,j; for (i=0.; i<2; i++) {__p[i] = 0; for (j=0.; j<2; j++) __p[i] += __a[i][j]*__x[j]; }}
236: static void CircleInflate(PetscReal r1, PetscReal r2, PetscReal r0, PetscInt num_sections, PetscReal x, PetscReal y,
237: PetscReal *outX, PetscReal *outY)
238: {
239: PetscReal rr = PetscSqrtReal(x*x + y*y), outfact, efact;
240: if (rr < r1 + PETSC_SQRT_MACHINE_EPSILON) {
241: *outX = x; *outY = y;
242: } else {
243: const PetscReal xy[2] = {x,y}, sinphi=y/rr, cosphi=x/rr;
244: PetscReal cth,sth,xyprime[2],Rth[2][2],rotcos,newrr;
245: if (num_sections==2) {
246: rotcos = 0.70710678118654;
247: outfact = 1.5; efact = 2.5;
248: /* rotate normalized vector into [-pi/4,pi/4) */
249: if (sinphi >= 0.) { /* top cell, -pi/2 */
250: cth = 0.707106781186548; sth = -0.707106781186548;
251: } else { /* bottom cell -pi/8 */
252: cth = 0.707106781186548; sth = .707106781186548;
253: }
254: } else if (num_sections==3) {
255: rotcos = 0.86602540378443;
256: outfact = 1.5; efact = 2.5;
257: /* rotate normalized vector into [-pi/6,pi/6) */
258: if (sinphi >= 0.5) { /* top cell, -pi/3 */
259: cth = 0.5; sth = -0.866025403784439;
260: } else if (sinphi >= -.5) { /* mid cell 0 */
261: cth = 1.; sth = .0;
262: } else { /* bottom cell +pi/3 */
263: cth = 0.5; sth = 0.866025403784439;
264: }
265: } else if (num_sections==4) {
266: rotcos = 0.9238795325112;
267: outfact = 1.5; efact = 3;
268: /* rotate normalized vector into [-pi/8,pi/8) */
269: if (sinphi >= 0.707106781186548) { /* top cell, -3pi/8 */
270: cth = 0.38268343236509; sth = -0.923879532511287;
271: } else if (sinphi >= 0.) { /* mid top cell -pi/8 */
272: cth = 0.923879532511287; sth = -.38268343236509;
273: } else if (sinphi >= -0.707106781186548) { /* mid bottom cell + pi/8 */
274: cth = 0.923879532511287; sth = 0.38268343236509;
275: } else { /* bottom cell + 3pi/8 */
276: cth = 0.38268343236509; sth = .923879532511287;
277: }
278: } else {
279: cth = 0.; sth = 0.; rotcos = 0; efact = 0;
280: }
281: Rth[0][0] = cth; Rth[0][1] =-sth;
282: Rth[1][0] = sth; Rth[1][1] = cth;
283: MATVEC2(Rth,xy,xyprime);
284: if (num_sections==2) {
285: newrr = xyprime[0]/rotcos;
286: } else {
287: PetscReal newcosphi=xyprime[0]/rr, rin = r1, rout = rr - rin;
288: PetscReal routmax = r0*rotcos/newcosphi - rin, nroutmax = r0 - rin, routfrac = rout/routmax;
289: newrr = rin + routfrac*nroutmax;
290: }
291: *outX = cosphi*newrr; *outY = sinphi*newrr;
292: /* grade */
293: PetscReal fact,tt,rs,re, rr = PetscSqrtReal(PetscSqr(*outX) + PetscSqr(*outY));
294: if (rr > r2) { rs = r2; re = r0; fact = outfact;} /* outer zone */
295: else { rs = r1; re = r2; fact = efact;} /* electron zone */
296: tt = (rs + PetscPowReal((rr - rs)/(re - rs),fact) * (re-rs)) / rr;
297: *outX *= tt;
298: *outY *= tt;
299: }
300: }
302: static PetscErrorCode GeometryDMLandau(DM base, PetscInt point, PetscInt dim, const PetscReal abc[], PetscReal xyz[], void *a_ctx)
303: {
304: LandauCtx *ctx = (LandauCtx*)a_ctx;
305: PetscReal r = abc[0], z = abc[1];
306: if (ctx->inflate) {
307: PetscReal absR, absZ;
308: absR = PetscAbsReal(r);
309: absZ = PetscAbsReal(z);
310: CircleInflate(ctx->i_radius,ctx->e_radius,ctx->radius,ctx->num_sections,absR,absZ,&absR,&absZ);
311: r = (r > 0) ? absR : -absR;
312: z = (z > 0) ? absZ : -absZ;
313: }
314: xyz[0] = r;
315: xyz[1] = z;
316: if (dim==3) xyz[2] = abc[2];
318: return(0);
319: }
321: static PetscErrorCode ErrorIndicator_Simple(PetscInt dim, PetscReal volume, PetscReal x[], PetscInt Nc, const PetscInt Nf[], const PetscScalar u[], const PetscScalar u_x[], PetscReal *error, void *actx)
322: {
323: PetscReal err = 0.0;
324: PetscInt f = *(PetscInt*)actx, j;
326: for (j = 0; j < dim; ++j) {
327: err += PetscSqr(PetscRealPart(u_x[f*dim+j]));
328: }
329: err = PetscRealPart(u[f]); /* just use rho */
330: *error = volume * err; /* * (ctx->axisymmetric ? 2.*PETSC_PI * r : 1); */
331: return(0);
332: }
334: static PetscErrorCode LandauDMCreateVMesh(MPI_Comm comm, const PetscInt dim, const char prefix[], LandauCtx *ctx, DM *dm)
335: {
337: PetscReal radius = ctx->radius;
338: size_t len;
339: char fname[128] = ""; /* we can add a file if we want */
342: /* create DM */
343: PetscStrlen(fname, &len);
344: if (len) {
345: PetscInt dim2;
346: DMPlexCreateFromFile(comm, fname, ctx->interpolate, dm);
347: DMGetDimension(*dm, &dim2);
348: } else { /* p4est, quads */
349: /* Create plex mesh of Landau domain */
350: if (!ctx->sphere) {
351: PetscInt cells[] = {2,2,2};
352: PetscReal lo[] = {-radius,-radius,-radius}, hi[] = {radius,radius,radius};
353: DMBoundaryType periodicity[3] = {DM_BOUNDARY_NONE, dim==2 ? DM_BOUNDARY_NONE : DM_BOUNDARY_NONE, DM_BOUNDARY_NONE};
354: if (dim==2) { lo[0] = 0; cells[0] = 1; }
355: else if (ctx->quarter3DDomain) { lo[0] = lo[1] = 0; cells[0] = cells[1] = 2; }
356: DMPlexCreateBoxMesh(comm, dim, PETSC_FALSE, cells, lo, hi, periodicity, PETSC_TRUE, dm);
357: DMLocalizeCoordinates(*dm); /* needed for periodic */
358: if (dim==3) {PetscObjectSetName((PetscObject) *dm, "cube");}
359: else {PetscObjectSetName((PetscObject) *dm, "half-plane");}
360: } else if (dim==2) {
361: PetscInt numCells,cells[16][4],i,j;
362: PetscInt numVerts;
363: PetscReal inner_radius1 = ctx->i_radius, inner_radius2 = ctx->e_radius;
364: PetscReal *flatCoords = NULL;
365: PetscInt *flatCells = NULL, *pcell;
366: if (ctx->num_sections==2) {
367: #if 1
368: numCells = 5;
369: numVerts = 10;
370: int cells2[][4] = { {0,1,4,3},
371: {1,2,5,4},
372: {3,4,7,6},
373: {4,5,8,7},
374: {6,7,8,9} };
375: for (i = 0; i < numCells; i++) for (j = 0; j < 4; j++) cells[i][j] = cells2[i][j];
376: PetscMalloc2(numVerts * 2, &flatCoords, numCells * 4, &flatCells);
377: {
378: PetscReal (*coords)[2] = (PetscReal (*) [2]) flatCoords;
379: for (j = 0; j < numVerts-1; j++) {
380: PetscReal z, r, theta = -PETSC_PI/2 + (j%3) * PETSC_PI/2;
381: PetscReal rad = (j >= 6) ? inner_radius1 : (j >= 3) ? inner_radius2 : ctx->radius;
382: z = rad * PetscSinReal(theta);
383: coords[j][1] = z;
384: r = rad * PetscCosReal(theta);
385: coords[j][0] = r;
386: }
387: coords[numVerts-1][0] = coords[numVerts-1][1] = 0;
388: }
389: #else
390: numCells = 4;
391: numVerts = 8;
392: static int cells2[][4] = {{0,1,2,3},
393: {4,5,1,0},
394: {5,6,2,1},
395: {6,7,3,2}};
396: for (i = 0; i < numCells; i++) for (j = 0; j < 4; j++) cells[i][j] = cells2[i][j];
397: PetscMalloc2(numVerts * 2, &flatCoords, numCells * 4, &flatCells);
398: {
399: PetscReal (*coords)[2] = (PetscReal (*) [2]) flatCoords;
400: PetscInt j;
401: for (j = 0; j < 8; j++) {
402: PetscReal z, r;
403: PetscReal theta = -PETSC_PI/2 + (j%4) * PETSC_PI/3.;
404: PetscReal rad = ctx->radius * ((j < 4) ? 0.5 : 1.0);
405: z = rad * PetscSinReal(theta);
406: coords[j][1] = z;
407: r = rad * PetscCosReal(theta);
408: coords[j][0] = r;
409: }
410: }
411: #endif
412: } else if (ctx->num_sections==3) {
413: numCells = 7;
414: numVerts = 12;
415: int cells2[][4] = { {0,1,5,4},
416: {1,2,6,5},
417: {2,3,7,6},
418: {4,5,9,8},
419: {5,6,10,9},
420: {6,7,11,10},
421: {8,9,10,11} };
422: for (i = 0; i < numCells; i++) for (j = 0; j < 4; j++) cells[i][j] = cells2[i][j];
423: PetscMalloc2(numVerts * 2, &flatCoords, numCells * 4, &flatCells);
424: {
425: PetscReal (*coords)[2] = (PetscReal (*) [2]) flatCoords;
426: for (j = 0; j < numVerts; j++) {
427: PetscReal z, r, theta = -PETSC_PI/2 + (j%4) * PETSC_PI/3;
428: PetscReal rad = (j >= 8) ? inner_radius1 : (j >= 4) ? inner_radius2 : ctx->radius;
429: z = rad * PetscSinReal(theta);
430: coords[j][1] = z;
431: r = rad * PetscCosReal(theta);
432: coords[j][0] = r;
433: }
434: }
435: } else if (ctx->num_sections==4) {
436: numCells = 10;
437: numVerts = 16;
438: int cells2[][4] = { {0,1,6,5},
439: {1,2,7,6},
440: {2,3,8,7},
441: {3,4,9,8},
442: {5,6,11,10},
443: {6,7,12,11},
444: {7,8,13,12},
445: {8,9,14,13},
446: {10,11,12,15},
447: {12,13,14,15}};
448: for (i = 0; i < numCells; i++) for (j = 0; j < 4; j++) cells[i][j] = cells2[i][j];
449: PetscMalloc2(numVerts * 2, &flatCoords, numCells * 4, &flatCells);
450: {
451: PetscReal (*coords)[2] = (PetscReal (*) [2]) flatCoords;
452: for (j = 0; j < numVerts-1; j++) {
453: PetscReal z, r, theta = -PETSC_PI/2 + (j%5) * PETSC_PI/4;
454: PetscReal rad = (j >= 10) ? inner_radius1 : (j >= 5) ? inner_radius2 : ctx->radius;
455: z = rad * PetscSinReal(theta);
456: coords[j][1] = z;
457: r = rad * PetscCosReal(theta);
458: coords[j][0] = r;
459: }
460: coords[numVerts-1][0] = coords[numVerts-1][1] = 0;
461: }
462: } else {
463: numCells = 0;
464: numVerts = 0;
465: }
466: for (j = 0, pcell = flatCells; j < numCells; j++, pcell += 4) {
467: pcell[0] = cells[j][0]; pcell[1] = cells[j][1];
468: pcell[2] = cells[j][2]; pcell[3] = cells[j][3];
469: }
470: DMPlexCreateFromCellListPetsc(comm,2,numCells,numVerts,4,ctx->interpolate,flatCells,2,flatCoords,dm);
471: PetscFree2(flatCoords,flatCells);
472: PetscObjectSetName((PetscObject) *dm, "semi-circle");
473: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Velocity space meshes does not support cubed sphere");
474: }
475: PetscObjectSetOptionsPrefix((PetscObject)*dm,prefix);
477: DMSetFromOptions(*dm); /* Plex refine */
479: { /* p4est? */
480: char convType[256];
481: PetscBool flg;
482: PetscOptionsBegin(PETSC_COMM_WORLD, prefix, "Mesh conversion options", "DMPLEX");
483: PetscOptionsFList("-dm_landau_type","Convert DMPlex to another format (should not be Plex!)","ex6f.c",DMList,DMPLEX,convType,256,&flg);
484: PetscOptionsEnd();
485: if (flg) {
486: DM dmforest;
487: DMConvert(*dm,convType,&dmforest);
488: if (dmforest) {
489: PetscBool isForest;
490: PetscObjectSetOptionsPrefix((PetscObject)dmforest,prefix);
491: DMIsForest(dmforest,&isForest);
492: if (isForest) {
493: if (ctx->sphere && ctx->inflate) {
494: DMForestSetBaseCoordinateMapping(dmforest,GeometryDMLandau,ctx);
495: }
496: DMDestroy(dm);
497: *dm = dmforest;
498: ctx->errorIndicator = ErrorIndicator_Simple; /* flag for Forest */
499: } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_USER, "Converted to non Forest?");
500: } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_USER, "Convert failed?");
501: }
502: }
503: PetscObjectSetName((PetscObject) *dm, "Mesh");
504: return(0);
505: }
507: static PetscErrorCode SetupDS(DM dm, PetscInt dim, LandauCtx *ctx)
508: {
509: PetscErrorCode ierr;
510: PetscInt ii;
512: for (ii=0;ii<ctx->num_species;ii++) {
513: char buf[256];
514: if (ii==0) PetscSNPrintf(buf, 256, "e");
515: else {PetscSNPrintf(buf, 256, "i%D", ii);}
516: /* Setup Discretization - FEM */
517: PetscFECreateDefault(PetscObjectComm((PetscObject) dm), dim, 1, ctx->simplex, NULL, PETSC_DECIDE, &ctx->fe[ii]);
518: PetscObjectSetName((PetscObject) ctx->fe[ii], buf);
519: DMSetField(dm, ii, NULL, (PetscObject) ctx->fe[ii]);
520: }
521: DMCreateDS(dm);
522: if (1) {
523: PetscInt ii;
524: PetscSection section;
525: DMGetSection(dm, §ion);
526: for (ii=0;ii<ctx->num_species;ii++){
527: char buf[256];
528: if (ii==0) PetscSNPrintf(buf, 256, "se");
529: else PetscSNPrintf(buf, 256, "si%D", ii);
530: PetscSectionSetComponentName(section, ii, 0, buf);
531: }
532: }
533: return(0);
534: }
536: /* Define a Maxwellian function for testing out the operator. */
538: /* Using cartesian velocity space coordinates, the particle */
539: /* density, [1/m^3], is defined according to */
541: /* $$ n=\int_{R^3} dv^3 \left(\frac{m}{2\pi T}\right)^{3/2}\exp [- mv^2/(2T)] $$ */
543: /* Using some constant, c, we normalize the velocity vector into a */
544: /* dimensionless variable according to v=c*x. Thus the density, $n$, becomes */
546: /* $$ n=\int_{R^3} dx^3 \left(\frac{mc^2}{2\pi T}\right)^{3/2}\exp [- mc^2/(2T)*x^2] $$ */
548: /* Defining $\theta=2T/mc^2$, we thus find that the probability density */
549: /* for finding the particle within the interval in a box dx^3 around x is */
551: /* f(x;\theta)=\left(\frac{1}{\pi\theta}\right)^{3/2} \exp [ -x^2/\theta ] */
553: typedef struct {
554: LandauCtx *ctx;
555: PetscReal kT_m;
556: PetscReal n;
557: PetscReal shift;
558: } MaxwellianCtx;
560: static PetscErrorCode maxwellian(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf_dummy, PetscScalar *u, void *actx)
561: {
562: MaxwellianCtx *mctx = (MaxwellianCtx*)actx;
563: LandauCtx *ctx = mctx->ctx;
564: PetscInt i;
565: PetscReal v2 = 0, theta = 2*mctx->kT_m/(ctx->v_0*ctx->v_0); /* theta = 2kT/mc^2 */
567: /* compute the exponents, v^2 */
568: for (i = 0; i < dim; ++i) v2 += x[i]*x[i];
569: /* evaluate the Maxwellian */
570: u[0] = mctx->n*PetscPowReal(PETSC_PI*theta,-1.5)*(PetscExpReal(-v2/theta));
571: if (mctx->shift!=0.) {
572: v2 = 0;
573: for (i = 0; i < dim-1; ++i) v2 += x[i]*x[i];
574: v2 += (x[dim-1]-mctx->shift)*(x[dim-1]-mctx->shift);
575: /* evaluate the shifted Maxwellian */
576: u[0] += mctx->n*PetscPowReal(PETSC_PI*theta,-1.5)*(PetscExpReal(-v2/theta));
577: }
578: return(0);
579: }
581: /*@
582: LandauAddMaxwellians - Add a Maxwellian distribution to a state
584: Collective on X
586: Input Parameters:
587: . dm - The mesh
588: + time - Current time
589: - temps - Temperatures of each species
590: . ns - Number density of each species
591: + actx - Landau context
593: Output Parameter:
594: . X - The state
596: Level: beginner
598: .keywords: mesh
599: .seealso: LandauCreateVelocitySpace()
600: @*/
601: PetscErrorCode LandauAddMaxwellians(DM dm, Vec X, PetscReal time, PetscReal temps[], PetscReal ns[], void *actx)
602: {
603: LandauCtx *ctx = (LandauCtx*)actx;
604: PetscErrorCode (*initu[LANDAU_MAX_SPECIES])(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar [], void *);
605: PetscErrorCode ierr,ii;
606: PetscInt dim;
607: MaxwellianCtx *mctxs[LANDAU_MAX_SPECIES], data[LANDAU_MAX_SPECIES];
610: DMGetDimension(dm, &dim);
611: if (!ctx) { DMGetApplicationContext(dm, &ctx); }
612: for (ii=0;ii<ctx->num_species;ii++) {
613: mctxs[ii] = &data[ii];
614: data[ii].ctx = ctx;
615: data[ii].kT_m = ctx->k*temps[ii]/ctx->masses[ii]; /* kT/m */
616: data[ii].n = ns[ii];
617: initu[ii] = maxwellian;
618: data[ii].shift = 0;
619: }
620: data[0].shift = ctx->electronShift;
621: /* need to make ADD_ALL_VALUES work - TODO */
622: DMProjectFunction(dm, time, initu, (void**)mctxs, INSERT_ALL_VALUES, X);
623: return(0);
624: }
626: /*
627: LandauSetInitialCondition - Addes Maxwellians with context
629: Collective on X
631: Input Parameters:
632: . dm - The mesh
633: + actx - Landau context with T and n
635: Output Parameter:
636: . X - The state
638: Level: beginner
640: .keywords: mesh
641: .seealso: LandauCreateVelocitySpace(), LandauAddMaxwellians()
642: */
643: static PetscErrorCode LandauSetInitialCondition(DM dm, Vec X, void *actx)
644: {
645: LandauCtx *ctx = (LandauCtx*)actx;
648: if (!ctx) { DMGetApplicationContext(dm, &ctx); }
649: VecZeroEntries(X);
650: LandauAddMaxwellians(dm, X, 0.0, ctx->thermal_temps, ctx->n, ctx);
651: return(0);
652: }
654: static PetscErrorCode adaptToleranceFEM(PetscFE fem, Vec sol, PetscReal refineTol[], PetscReal coarsenTol[], PetscInt type, LandauCtx *ctx, DM *newDM)
655: {
656: DM dm, plex, adaptedDM = NULL;
657: PetscDS prob;
658: PetscBool isForest;
659: PetscQuadrature quad;
660: PetscInt Nq, *Nb, cStart, cEnd, c, dim, qj, k;
661: DMLabel adaptLabel = NULL;
662: PetscErrorCode ierr;
665: VecGetDM(sol, &dm);
666: DMCreateDS(dm);
667: DMGetDS(dm, &prob);
668: DMGetDimension(dm, &dim);
669: DMIsForest(dm, &isForest);
670: DMConvert(dm, DMPLEX, &plex);
671: DMPlexGetHeightStratum(plex,0,&cStart,&cEnd);
672: DMLabelCreate(PETSC_COMM_SELF,"adapt",&adaptLabel);
673: PetscFEGetQuadrature(fem, &quad);
674: PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL);
675: if (Nq >LANDAU_MAX_NQ) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Order too high. Nq = %D > LANDAU_MAX_NQ (%D)",Nq,LANDAU_MAX_NQ);
676: PetscDSGetDimensions(prob, &Nb);
677: if (type==4) {
678: for (c = cStart; c < cEnd; c++) {
679: DMLabelSetValue(adaptLabel, c, DM_ADAPT_REFINE);
680: }
681: PetscInfo1(sol, "Phase:%s: Uniform refinement\n","adaptToleranceFEM");
682: } else if (type==2) {
683: PetscInt rCellIdx[8], eCellIdx[64], iCellIdx[64], eMaxIdx = -1, iMaxIdx = -1, nr = 0, nrmax = (dim==3 && !ctx->quarter3DDomain) ? 8 : 2;
684: PetscReal minRad = PETSC_INFINITY, r, eMinRad = PETSC_INFINITY, iMinRad = PETSC_INFINITY;
685: for (c = 0; c < 64; c++) { eCellIdx[c] = iCellIdx[c] = -1; }
686: for (c = cStart; c < cEnd; c++) {
687: PetscReal tt, v0[LANDAU_MAX_NQ*3], detJ[LANDAU_MAX_NQ];
688: DMPlexComputeCellGeometryFEM(plex, c, quad, v0, NULL, NULL, detJ);
689: for (qj = 0; qj < Nq; ++qj) {
690: tt = PetscSqr(v0[dim*qj+0]) + PetscSqr(v0[dim*qj+1]) + PetscSqr(((dim==3) ? v0[dim*qj+2] : 0));
691: r = PetscSqrtReal(tt);
692: if (r < minRad - PETSC_SQRT_MACHINE_EPSILON*10.) {
693: minRad = r;
694: nr = 0;
695: rCellIdx[nr++]= c;
696: PetscInfo4(sol, "\t\tPhase: adaptToleranceFEM Found first inner r=%e, cell %D, qp %D/%D\n", r, c, qj+1, Nq);
697: } else if ((r-minRad) < PETSC_SQRT_MACHINE_EPSILON*100. && nr < nrmax) {
698: for (k=0;k<nr;k++) if (c == rCellIdx[k]) break;
699: if (k==nr) {
700: rCellIdx[nr++]= c;
701: PetscInfo5(sol, "\t\t\tPhase: adaptToleranceFEM Found another inner r=%e, cell %D, qp %D/%D, d=%e\n", r, c, qj+1, Nq, r-minRad);
702: }
703: }
704: if (ctx->sphere) {
705: if ((tt=r-ctx->e_radius) > 0) {
706: PetscInfo2(sol, "\t\t\t %D cell r=%g\n",c,tt);
707: if (tt < eMinRad - PETSC_SQRT_MACHINE_EPSILON*100.) {
708: eMinRad = tt;
709: eMaxIdx = 0;
710: eCellIdx[eMaxIdx++] = c;
711: } else if (eMaxIdx > 0 && (tt-eMinRad) <= PETSC_SQRT_MACHINE_EPSILON && c != eCellIdx[eMaxIdx-1]) {
712: eCellIdx[eMaxIdx++] = c;
713: }
714: }
715: if ((tt=r-ctx->i_radius) > 0) {
716: if (tt < iMinRad - 1.e-5) {
717: iMinRad = tt;
718: iMaxIdx = 0;
719: iCellIdx[iMaxIdx++] = c;
720: } else if (iMaxIdx > 0 && (tt-iMinRad) <= PETSC_SQRT_MACHINE_EPSILON && c != iCellIdx[iMaxIdx-1]) {
721: iCellIdx[iMaxIdx++] = c;
722: }
723: }
724: }
725: }
726: }
727: for (k=0;k<nr;k++) {
728: DMLabelSetValue(adaptLabel, rCellIdx[k], DM_ADAPT_REFINE);
729: }
730: if (ctx->sphere) {
731: for (c = 0; c < eMaxIdx; c++) {
732: DMLabelSetValue(adaptLabel, eCellIdx[c], DM_ADAPT_REFINE);
733: PetscInfo3(sol, "\t\tPhase:%s: refine sphere e cell %D r=%g\n","adaptToleranceFEM",eCellIdx[c],eMinRad);
734: }
735: for (c = 0; c < iMaxIdx; c++) {
736: DMLabelSetValue(adaptLabel, iCellIdx[c], DM_ADAPT_REFINE);
737: PetscInfo3(sol, "\t\tPhase:%s: refine sphere i cell %D r=%g\n","adaptToleranceFEM",iCellIdx[c],iMinRad);
738: }
739: }
740: PetscInfo4(sol, "Phase:%s: Adaptive refine origin cells %D,%D r=%g\n","adaptToleranceFEM",rCellIdx[0],rCellIdx[1],minRad);
741: } else if (type==0 || type==1 || type==3) { /* refine along r=0 axis */
742: PetscScalar *coef = NULL;
743: Vec coords;
744: PetscInt csize,Nv,d,nz;
745: DM cdm;
746: PetscSection cs;
747: DMGetCoordinatesLocal(dm, &coords);
748: DMGetCoordinateDM(dm, &cdm);
749: DMGetLocalSection(cdm, &cs);
750: for (c = cStart; c < cEnd; c++) {
751: PetscInt doit = 0, outside = 0;
752: DMPlexVecGetClosure(cdm, cs, coords, c, &csize, &coef);
753: Nv = csize/dim;
754: for (nz = d = 0; d < Nv; d++) {
755: PetscReal z = PetscRealPart(coef[d*dim + (dim-1)]), x = PetscSqr(PetscRealPart(coef[d*dim + 0])) + ((dim==3) ? PetscSqr(PetscRealPart(coef[d*dim + 1])) : 0);
756: x = PetscSqrtReal(x);
757: if (x < PETSC_MACHINE_EPSILON*10. && PetscAbsReal(z)<PETSC_MACHINE_EPSILON*10.) doit = 1; /* refine origin */
758: else if (type==0 && (z < -PETSC_MACHINE_EPSILON*10. || z > ctx->re_radius+PETSC_MACHINE_EPSILON*10.)) outside++; /* first pass don't refine bottom */
759: else if (type==1 && (z > ctx->vperp0_radius1 || z < -ctx->vperp0_radius1)) outside++; /* don't refine outside electron refine radius */
760: else if (type==3 && (z > ctx->vperp0_radius2 || z < -ctx->vperp0_radius2)) outside++; /* don't refine outside ion refine radius */
761: if (x < PETSC_MACHINE_EPSILON*10.) nz++;
762: }
763: DMPlexVecRestoreClosure(cdm, cs, coords, c, &csize, &coef);
764: if (doit || (outside<Nv && nz)) {
765: DMLabelSetValue(adaptLabel, c, DM_ADAPT_REFINE);
766: }
767: }
768: PetscInfo1(sol, "Phase:%s: RE refinement\n","adaptToleranceFEM");
769: }
770: /* VecDestroy(&locX); */
771: DMDestroy(&plex);
772: DMAdaptLabel(dm, adaptLabel, &adaptedDM);
773: DMLabelDestroy(&adaptLabel);
774: *newDM = adaptedDM;
775: if (adaptedDM) {
776: if (isForest) {
777: DMForestSetAdaptivityForest(adaptedDM,NULL);
778: }
779: DMConvert(adaptedDM, DMPLEX, &plex);
780: DMPlexGetHeightStratum(plex,0,&cStart,&cEnd);
781: PetscInfo2(sol, "\tPhase: adaptToleranceFEM: %D cells, %d total quadrature points\n",cEnd-cStart,Nq*(cEnd-cStart));
782: DMDestroy(&plex);
783: }
784: return(0);
785: }
787: static PetscErrorCode adapt(DM *dm, LandauCtx *ctx, Vec *uu)
788: {
789: PetscErrorCode ierr;
790: PetscInt type, limits[5] = {ctx->numRERefine,ctx->nZRefine1,ctx->maxRefIts,ctx->nZRefine2,ctx->postAMRRefine};
791: PetscInt adaptIter;
794: for (type=0;type<5;type++) {
795: for (adaptIter = 0; adaptIter<limits[type];adaptIter++) {
796: DM dmNew = NULL;
797: adaptToleranceFEM(ctx->fe[0], *uu, ctx->refineTol, ctx->coarsenTol, type, ctx, &dmNew);
798: if (!dmNew) {
799: exit(113);
800: break;
801: } else {
802: DMDestroy(dm);
803: VecDestroy(uu);
804: DMCreateGlobalVector(dmNew,uu);
805: PetscObjectSetName((PetscObject) *uu, "u");
806: LandauSetInitialCondition(dmNew, *uu, ctx);
807: *dm = dmNew;
808: }
809: }
810: }
811: return(0);
812: }
814: static PetscErrorCode ProcessOptions(LandauCtx *ctx, const char prefix[])
815: {
816: PetscErrorCode ierr;
817: PetscBool flg, sph_flg;
818: PetscInt ii,nt,nm,nc;
819: DM dummy;
822: DMCreate(PETSC_COMM_WORLD,&dummy);
823: /* get options - initialize context */
824: ctx->normJ = 0;
825: ctx->verbose = 1;
826: ctx->interpolate = PETSC_TRUE;
827: ctx->simplex = PETSC_FALSE;
828: ctx->sphere = PETSC_FALSE;
829: ctx->inflate = PETSC_FALSE;
830: ctx->electronShift = 0;
831: ctx->errorIndicator = NULL;
832: ctx->radius = 5.; /* electron thermal radius (velocity) */
833: ctx->re_radius = 0.;
834: ctx->vperp0_radius1 = 0;
835: ctx->vperp0_radius2 = 0;
836: ctx->e_radius = .1;
837: ctx->i_radius = .01;
838: ctx->maxRefIts = 5;
839: ctx->postAMRRefine = 0;
840: ctx->nZRefine1 = 0;
841: ctx->nZRefine2 = 0;
842: ctx->numRERefine = 0;
843: ctx->num_sections = 3; /* 2, 3 or 4 */
844: /* species - [0] electrons, [1] one ion species eg, duetarium, [2] heavy impurity ion, ... */
845: ctx->charges[0] = -1; /* electron charge (MKS) */
846: ctx->masses[0] = 1/1835.5; /* temporary value in proton mass */
847: ctx->n[0] = 1;
848: ctx->thermal_temps[0] = 1;
849: /* constants, etc. */
850: ctx->epsilon0 = 8.8542e-12; /* permittivity of free space (MKS) F/m */
851: ctx->k = 1.38064852e-23; /* Boltzmann constant (MKS) J/K */
852: ctx->lnLam = 10; /* cross section ratio large - small angle collisions */
853: ctx->n_0 = 1.e20; /* typical plasma n, but could set it to 1 */
854: ctx->Ez = 0;
855: ctx->v_0 = 1; /* in electron thermal velocity */
856: ctx->subThreadBlockSize = 1; /* for device and maybe OMP */
857: ctx->quarter3DDomain = PETSC_FALSE;
858: PetscOptionsBegin(PETSC_COMM_WORLD, prefix, "Options for Fokker-Plank-Landau collision operator", "none");
859: {
860: char opstring[256];
861: #if defined(PETSC_HAVE_KOKKOS)
862: ctx->deviceType = LANDAU_KOKKOS;
863: PetscStrcpy(opstring,"kokkos");
864: #if !defined(PETSC_HAVE_CUDA)
865: ctx->subThreadBlockSize = 0;
866: #endif
867: #elif defined(PETSC_HAVE_CUDA)
868: ctx->deviceType = LANDAU_CUDA;
869: PetscStrcpy(opstring,"cuda");
870: #else
871: ctx->deviceType = LANDAU_CPU;
872: PetscStrcpy(opstring,"cpu");
873: ctx->subThreadBlockSize = 0;
874: #endif
875: PetscOptionsString("-dm_landau_device_type","Use kernels on 'cpu', 'cuda', or 'kokkos'","plexland.c",opstring,opstring,256,NULL);
876: PetscStrcmp("cpu",opstring,&flg);
877: if (flg) {
878: ctx->deviceType = LANDAU_CPU;
879: ctx->subThreadBlockSize = 0;
880: } else {
881: PetscStrcmp("cuda",opstring,&flg);
882: if (flg) ctx->deviceType = LANDAU_CUDA;
883: else {
884: PetscStrcmp("kokkos",opstring,&flg);
885: if (flg) ctx->deviceType = LANDAU_KOKKOS;
886: else SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"-dm_landau_device_type %s",opstring);
887: }
888: }
889: }
890: PetscOptionsReal("-dm_landau_electron_shift","Shift in thermal velocity of electrons","none",ctx->electronShift,&ctx->electronShift, NULL);
891: PetscOptionsBool("-dm_landau_sphere", "use sphere/semi-circle domain instead of rectangle", "plexland.c", ctx->sphere, &ctx->sphere, &sph_flg);
892: PetscOptionsBool("-dm_landau_inflate", "With sphere, inflate for curved edges (no AMR)", "plexland.c", ctx->inflate, &ctx->inflate, NULL);
893: /* PetscOptionsBool("-dm_landau_quarter_3d_domain", "Use symmetry in 3D to model 1/4 of domain", "plexland.c", ctx->quarter3DDomain, &ctx->quarter3DDomain, NULL); */
894: PetscOptionsInt("-dm_landau_amr_re_levels", "Number of levels to refine along v_perp=0, z>0", "plexland.c", ctx->numRERefine, &ctx->numRERefine, NULL);
895: PetscOptionsInt("-dm_landau_amr_z_refine1", "Number of levels to refine along v_perp=0", "plexland.c", ctx->nZRefine1, &ctx->nZRefine1, NULL);
896: PetscOptionsInt("-dm_landau_amr_z_refine2", "Number of levels to refine along v_perp=0", "plexland.c", ctx->nZRefine2, &ctx->nZRefine2, NULL);
897: PetscOptionsInt("-dm_landau_amr_levels_max", "Number of AMR levels of refinement around origin after r=0 refinements", "plexland.c", ctx->maxRefIts, &ctx->maxRefIts, NULL);
898: PetscOptionsInt("-dm_landau_amr_post_refine", "Number of levels to uniformly refine after AMR", "plexland.c", ctx->postAMRRefine, &ctx->postAMRRefine, NULL);
899: PetscOptionsInt("-dm_landau_verbose", "", "plexland.c", ctx->verbose, &ctx->verbose, NULL);
900: PetscOptionsReal("-dm_landau_re_radius","velocity range to refine on positive (z>0) r=0 axis for runaways","plexland.c",ctx->re_radius,&ctx->re_radius, &flg);
901: PetscOptionsReal("-dm_landau_z_radius1","velocity range to refine r=0 axis (for electrons)","plexland.c",ctx->vperp0_radius1,&ctx->vperp0_radius1, &flg);
902: PetscOptionsReal("-dm_landau_z_radius2","velocity range to refine r=0 axis (for ions) after origin AMR","plexland.c",ctx->vperp0_radius2,&ctx->vperp0_radius2, &flg);
903: PetscOptionsReal("-dm_landau_Ez","Initial parallel electric field in unites of Conner-Hastie criticle field","plexland.c",ctx->Ez,&ctx->Ez, NULL);
904: PetscOptionsReal("-dm_landau_n_0","Normalization constant for number density","plexland.c",ctx->n_0,&ctx->n_0, NULL);
905: PetscOptionsReal("-dm_landau_ln_lambda","Cross section parameter","plexland.c",ctx->lnLam,&ctx->lnLam, NULL);
906: PetscOptionsInt("-dm_landau_num_sections", "Number of tangential section in (2D) grid, 2, 3, of 4", "plexland.c", ctx->num_sections, &ctx->num_sections, NULL);
907: ctx->simplex = PETSC_FALSE;
908: /* get num species with tempurature*/
909: {
910: PetscReal arr[100];
911: nt = 100;
912: PetscOptionsRealArray("-dm_landau_thermal_temps", "Temperature of each species [e,i_0,i_1,...] in keV", "plexland.c", arr, &nt, &flg);
913: if (flg && nt > LANDAU_MAX_SPECIES) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"-thermal_temps ,t1,t2,.. number of species %D > MAX %D",nt,LANDAU_MAX_SPECIES);
914: }
915: nt = LANDAU_MAX_SPECIES;
916: for (ii=1;ii<LANDAU_MAX_SPECIES;ii++) {
917: ctx->thermal_temps[ii] = 1.;
918: ctx->charges[ii] = 1;
919: ctx->masses[ii] = 1;
920: ctx->n[ii] = (ii==1) ? 1 : 0;
921: }
922: PetscOptionsRealArray("-dm_landau_thermal_temps", "Temperature of each species [e,i_0,i_1,...] in keV (must be set to set number of species)", "plexland.c", ctx->thermal_temps, &nt, &flg);
923: if (flg) {
924: PetscInfo1(dummy, "num_species set to number of thermal temps provided (%D)\n",nt);
925: ctx->num_species = nt;
926: } else SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"-dm_landau_thermal_temps ,t1,t2,.. must be provided to set the number of species");
927: for (ii=0;ii<ctx->num_species;ii++) ctx->thermal_temps[ii] *= 1.1604525e7; /* convert to Kelvin */
928: nm = LANDAU_MAX_SPECIES-1;
929: PetscOptionsRealArray("-dm_landau_ion_masses", "Mass of each species in units of proton mass [i_0=2,i_1=40...]", "plexland.c", &ctx->masses[1], &nm, &flg);
930: if (flg && nm != ctx->num_species-1) {
931: SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"num ion masses %D != num species %D",nm,ctx->num_species-1);
932: }
933: nm = LANDAU_MAX_SPECIES;
934: PetscOptionsRealArray("-dm_landau_n", "Normalized (by -n_0) number density of each species", "plexland.c", ctx->n, &nm, &flg);
935: if (flg && nm != ctx->num_species) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"wrong num n: %D != num species %D",nm,ctx->num_species);
936: ctx->n_0 *= ctx->n[0]; /* normalized number density */
937: for (ii=1;ii<ctx->num_species;ii++) ctx->n[ii] = ctx->n[ii]/ctx->n[0];
938: ctx->n[0] = 1;
939: for (ii=0;ii<LANDAU_MAX_SPECIES;ii++) ctx->masses[ii] *= 1.6720e-27; /* scale by proton mass kg */
940: ctx->masses[0] = 9.10938356e-31; /* electron mass kg (should be about right already) */
941: ctx->m_0 = ctx->masses[0]; /* arbitrary reference mass, electrons */
942: PetscOptionsReal("-dm_landau_v_0","Velocity to normalize with in units of initial electrons thermal velocity (not recommended to change default)","plexland.c",ctx->v_0,&ctx->v_0, NULL);
943: ctx->v_0 *= PetscSqrtReal(ctx->k*ctx->thermal_temps[0]/(ctx->masses[0])); /* electron mean velocity in 1D (need 3D form in computing T from FE integral) */
944: nc = LANDAU_MAX_SPECIES-1;
945: PetscOptionsRealArray("-dm_landau_ion_charges", "Charge of each species in units of proton charge [i_0=2,i_1=18,...]", "plexland.c", &ctx->charges[1], &nc, &flg);
946: if (flg && nc != ctx->num_species-1) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"num charges %D != num species %D",nc,ctx->num_species-1);
947: for (ii=0;ii<LANDAU_MAX_SPECIES;ii++) ctx->charges[ii] *= 1.6022e-19; /* electron/proton charge (MKS) */
948: ctx->t_0 = 8*PETSC_PI*PetscSqr(ctx->epsilon0*ctx->m_0/PetscSqr(ctx->charges[0]))/ctx->lnLam/ctx->n_0*PetscPowReal(ctx->v_0,3); /* note, this t_0 makes nu[0,0]=1 */
949: /* geometry */
950: for (ii=0;ii<ctx->num_species;ii++) ctx->refineTol[ii] = PETSC_MAX_REAL;
951: for (ii=0;ii<ctx->num_species;ii++) ctx->coarsenTol[ii] = 0.;
952: ii = LANDAU_MAX_SPECIES;
953: PetscOptionsRealArray("-dm_landau_refine_tol","tolerance for refining cells in AMR","plexland.c",ctx->refineTol, &ii, &flg);
954: if (flg && ii != ctx->num_species) PetscInfo2(dummy, "Phase: Warning, #refine_tol %D != num_species %D\n",ii,ctx->num_species);
955: ii = LANDAU_MAX_SPECIES;
956: PetscOptionsRealArray("-dm_landau_coarsen_tol","tolerance for coarsening cells in AMR","plexland.c",ctx->coarsenTol, &ii, &flg);
957: if (flg && ii != ctx->num_species) PetscInfo2(dummy, "Phase: Warning, #coarsen_tol %D != num_species %D\n",ii,ctx->num_species);
958: PetscOptionsReal("-dm_landau_domain_radius","Phase space size in units of electron thermal velocity","plexland.c",ctx->radius,&ctx->radius, &flg);
959: if (flg && ctx->radius <= 0) { /* negative is ratio of c */
960: if (ctx->radius == 0) ctx->radius = 0.75;
961: else ctx->radius = -ctx->radius;
962: ctx->radius = ctx->radius*299792458.0/ctx->v_0;
963: PetscInfo1(dummy, "Change domain radius to %e\n",ctx->radius);
964: }
965: PetscOptionsReal("-dm_landau_i_radius","Ion thermal velocity, used for circular meshes","plexland.c",ctx->i_radius,&ctx->i_radius, &flg);
966: if (flg && !sph_flg) ctx->sphere = PETSC_TRUE; /* you gave me an ion radius but did not set sphere, user error really */
967: if (!flg) {
968: ctx->i_radius = 1.5*PetscSqrtReal(8*ctx->k*ctx->thermal_temps[1]/ctx->masses[1]/PETSC_PI)/ctx->v_0; /* normalized radius with thermal velocity of first ion */
969: }
970: PetscOptionsReal("-dm_landau_e_radius","Electron thermal velocity, used for circular meshes","plexland.c",ctx->e_radius,&ctx->e_radius, &flg);
971: if (flg && !sph_flg) ctx->sphere = PETSC_TRUE; /* you gave me an e radius but did not set sphere, user error really */
972: if (!flg) {
973: ctx->e_radius = 1.5*PetscSqrtReal(8*ctx->k*ctx->thermal_temps[0]/ctx->masses[0]/PETSC_PI)/ctx->v_0; /* normalized radius with thermal velocity of electrons */
974: }
975: if (ctx->sphere && (ctx->e_radius <= ctx->i_radius || ctx->radius <= ctx->e_radius)) SETERRQ3(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"bad radii: %g < %g < %g",ctx->i_radius,ctx->e_radius,ctx->radius);
976: PetscOptionsInt("-dm_landau_sub_thread_block_size", "Number of threads in CUDA integration point subblock", "plexland.c", ctx->subThreadBlockSize, &ctx->subThreadBlockSize, NULL);
977: #if !defined(PETSC_HAVE_KOKKOS)
978: if (ctx->subThreadBlockSize > LANDAU_MAX_SUB_THREAD_BLOCKS) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"num sub threads %D > MAX %d",ctx->subThreadBlockSize,LANDAU_MAX_SUB_THREAD_BLOCKS);
979: #endif
980: PetscOptionsEnd();
981: for (ii=ctx->num_species;ii<LANDAU_MAX_SPECIES;ii++) ctx->masses[ii] = ctx->thermal_temps[ii] = ctx->charges[ii] = 0;
982: if (ctx->verbose > 0) {
983: PetscPrintf(PETSC_COMM_WORLD, "masses: e=%10.3e; ions in proton mass units: %10.3e %10.3e ...\n",ctx->masses[0],ctx->masses[1]/1.6720e-27,ctx->num_species>2 ? ctx->masses[2]/1.6720e-27 : 0);
984: PetscPrintf(PETSC_COMM_WORLD, "charges: e=%10.3e; charges in elementary units: %10.3e %10.3e\n", ctx->charges[0],-ctx->charges[1]/ctx->charges[0],ctx->num_species>2 ? -ctx->charges[2]/ctx->charges[0] : 0);
985: PetscPrintf(PETSC_COMM_WORLD, "thermal T (K): e=%10.3e i=%10.3e imp=%10.3e. v_0=%10.3e n_0=%10.3e t_0=%10.3e domain=%10.3e\n",ctx->thermal_temps[0],ctx->thermal_temps[1],ctx->num_species>2 ? ctx->thermal_temps[2] : 0,ctx->v_0,ctx->n_0,ctx->t_0,ctx->radius);
986: }
987: DMDestroy(&dummy);
988: {
989: PetscMPIInt rank;
990: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
991: /* PetscLogStage setup_stage; */
992: PetscLogEventRegister("Landau Operator", DM_CLASSID, &ctx->events[0]); /* 0 */
993: PetscLogEventRegister(" Jac-vector", DM_CLASSID, &ctx->events[1]); /* 1 */
994: PetscLogEventRegister(" Jac-kern-init", DM_CLASSID, &ctx->events[3]); /* 3 */
995: PetscLogEventRegister(" Jac-kernel", DM_CLASSID, &ctx->events[4]); /* 4 */
996: PetscLogEventRegister(" Jac-kernel-post", DM_CLASSID, &ctx->events[5]); /* 5 */
997: PetscLogEventRegister(" Jac-assemble", DM_CLASSID, &ctx->events[6]); /* 6 */
998: PetscLogEventRegister(" Jac-end", DM_CLASSID, &ctx->events[7]); /* 7 */
999: PetscLogEventRegister(" Jac-geo-color", DM_CLASSID, &ctx->events[8]); /* 8 */
1000: PetscLogEventRegister(" Jac-cuda-sum", DM_CLASSID, &ctx->events[2]); /* 2 */
1001: PetscLogEventRegister("Landau Jacobian", DM_CLASSID, &ctx->events[9]); /* 9 */
1002: if (rank) { /* turn off output stuff for duplicate runs - do we need to add the prefix to all this? */
1003: PetscOptionsClearValue(NULL,"-snes_converged_reason");
1004: PetscOptionsClearValue(NULL,"-ksp_converged_reason");
1005: PetscOptionsClearValue(NULL,"-snes_monitor");
1006: PetscOptionsClearValue(NULL,"-ksp_monitor");
1007: PetscOptionsClearValue(NULL,"-ts_monitor");
1008: PetscOptionsClearValue(NULL,"-ts_adapt_monitor");
1009: PetscOptionsClearValue(NULL,"-dm_landau_amr_dm_view");
1010: PetscOptionsClearValue(NULL,"-dm_landau_amr_vec_view");
1011: PetscOptionsClearValue(NULL,"-dm_landau_pre_dm_view");
1012: PetscOptionsClearValue(NULL,"-dm_landau_pre_vec_view");
1013: PetscOptionsClearValue(NULL,"-info");
1014: }
1015: }
1016: return(0);
1017: }
1019: /*@C
1020: LandauCreateVelocitySpace - Create a DMPlex velocity space mesh
1022: Collective on comm
1024: Input Parameters:
1025: + comm - The MPI communicator
1026: . dim - velocity space dimension (2 for axisymmetric, 3 for full 3X + 3V solver)
1027: - prefix - prefix for options
1029: Output Parameter:
1030: . dm - The DM object representing the mesh
1031: + X - A vector (user destroys)
1032: - J - Optional matrix (object destroys)
1034: Level: beginner
1036: .keywords: mesh
1037: .seealso: DMPlexCreate(), LandauDestroyVelocitySpace()
1038: @*/
1039: PetscErrorCode LandauCreateVelocitySpace(MPI_Comm comm, PetscInt dim, const char prefix[], Vec *X, Mat *J, DM *dm)
1040: {
1041: PetscMPIInt size;
1043: LandauCtx *ctx;
1046: MPI_Comm_size(comm, &size);
1047: if (size!=1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Velocity space meshes should be serial (but should work in parallel)");
1048: if (dim!=2 && dim!=3) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Only 2D and 3D supported");
1049: ctx = (LandauCtx*)malloc(sizeof(LandauCtx));
1050: /* process options */
1051: ProcessOptions(ctx,prefix);
1052: /* Create Mesh */
1053: LandauDMCreateVMesh(comm, dim, prefix, ctx, dm);
1054: DMViewFromOptions(*dm,NULL,"-dm_landau_pre_dm_view");
1055: DMSetApplicationContext(*dm, ctx);
1056: /* create FEM */
1057: SetupDS(*dm,dim,ctx);
1058: /* set initial state */
1059: DMCreateGlobalVector(*dm,X);
1060: PetscObjectSetName((PetscObject) *X, "u");
1061: /* initial static refinement, no solve */
1062: LandauSetInitialCondition(*dm, *X, ctx);
1063: VecViewFromOptions(*X, NULL, "-dm_landau_pre_vec_view");
1064: /* forest refinement */
1065: if (ctx->errorIndicator) {
1066: /* AMR */
1067: adapt(dm,ctx,X);
1068: DMViewFromOptions(*dm,NULL,"-dm_landau_amr_dm_view");
1069: VecViewFromOptions(*X, NULL, "-dm_landau_amr_vec_view");
1070: }
1071: DMSetApplicationContext(*dm, ctx);
1072: ctx->dmv = *dm;
1073: DMCreateMatrix(ctx->dmv, &ctx->J);
1074: if (J) *J = ctx->J;
1075: return(0);
1076: }
1078: /*@
1079: LandauDestroyVelocitySpace - Destroy a DMPlex velocity space mesh
1081: Collective on dm
1083: Input/Output Parameters:
1084: . dm - the dm to destroy
1086: Level: beginner
1088: .keywords: mesh
1089: .seealso: LandauCreateVelocitySpace()
1090: @*/
1091: PetscErrorCode LandauDestroyVelocitySpace(DM *dm)
1092: {
1093: PetscErrorCode ierr,ii;
1094: LandauCtx *ctx;
1095: PetscContainer container = NULL;
1097: DMGetApplicationContext(*dm, &ctx);
1098: PetscObjectQuery((PetscObject)ctx->J,"coloring", (PetscObject*)&container);
1099: if (container) {
1100: PetscContainerDestroy(&container);
1101: }
1102: MatDestroy(&ctx->M);
1103: MatDestroy(&ctx->J);
1104: for (ii=0;ii<ctx->num_species;ii++) {
1105: PetscFEDestroy(&ctx->fe[ii]);
1106: }
1107: free(ctx);
1108: DMDestroy(dm);
1109: return(0);
1110: }
1112: /* < v, ru > */
1113: static void f0_s_den(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1114: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1115: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1116: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar *f0)
1117: {
1118: PetscInt ii = (PetscInt)PetscRealPart(constants[0]);
1119: f0[0] = u[ii];
1120: }
1122: /* < v, ru > */
1123: static void f0_s_mom(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1124: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1125: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1126: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar *f0)
1127: {
1128: PetscInt ii = (PetscInt)PetscRealPart(constants[0]), jj = (PetscInt)PetscRealPart(constants[1]);
1129: f0[0] = x[jj]*u[ii]; /* x momentum */
1130: }
1132: static void f0_s_v2(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1133: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1134: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1135: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar *f0)
1136: {
1137: PetscInt i, ii = (PetscInt)PetscRealPart(constants[0]);
1138: double tmp1 = 0.;
1139: for (i = 0; i < dim; ++i) tmp1 += x[i]*x[i];
1140: f0[0] = tmp1*u[ii];
1141: }
1143: /* < v, ru > */
1144: static void f0_s_rden(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1145: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1146: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1147: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar *f0)
1148: {
1149: PetscInt ii = (PetscInt)PetscRealPart(constants[0]);
1150: f0[0] = 2.*PETSC_PI*x[0]*u[ii];
1151: }
1153: /* < v, ru > */
1154: static void f0_s_rmom(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1155: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1156: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1157: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar *f0)
1158: {
1159: PetscInt ii = (PetscInt)PetscRealPart(constants[0]);
1160: f0[0] = 2.*PETSC_PI*x[0]*x[1]*u[ii];
1161: }
1163: static void f0_s_rv2(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1164: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1165: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1166: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar *f0)
1167: {
1168: PetscInt ii = (PetscInt)PetscRealPart(constants[0]);
1169: f0[0] = 2.*PETSC_PI*x[0]*(x[0]*x[0] + x[1]*x[1])*u[ii];
1170: }
1172: /*@
1173: LandauPrintNorms - collects moments and prints them
1175: Collective on dm
1177: Input Parameters:
1178: + X - the state
1179: - stepi - current step to print
1181: Level: beginner
1183: .keywords: mesh
1184: .seealso: LandauCreateVelocitySpace()
1185: @*/
1186: PetscErrorCode LandauPrintNorms(Vec X, PetscInt stepi)
1187: {
1189: LandauCtx *ctx;
1190: PetscDS prob;
1191: DM plex,dm;
1192: PetscInt cStart, cEnd, dim, ii;
1193: PetscScalar xmomentumtot=0, ymomentumtot=0, zmomentumtot=0, energytot=0, densitytot=0, tt[LANDAU_MAX_SPECIES];
1194: PetscScalar xmomentum[LANDAU_MAX_SPECIES], ymomentum[LANDAU_MAX_SPECIES], zmomentum[LANDAU_MAX_SPECIES], energy[LANDAU_MAX_SPECIES], density[LANDAU_MAX_SPECIES];
1197: VecGetDM(X, &dm);
1198: if (!dm) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "no DM");
1199: DMGetDimension(dm, &dim);
1200: DMGetApplicationContext(dm, &ctx);
1201: if (!ctx) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "no context");
1202: DMConvert(ctx->dmv, DMPLEX, &plex);
1203: DMCreateDS(plex);
1204: DMGetDS(plex, &prob);
1205: /* print momentum and energy */
1206: for (ii=0;ii<ctx->num_species;ii++) {
1207: PetscScalar user[2] = { (PetscScalar)ii, (PetscScalar)ctx->charges[ii]};
1208: PetscDSSetConstants(prob, 2, user);
1209: if (dim==2) { /* 2/3X + 3V (cylindrical coordinates) */
1210: PetscDSSetObjective(prob, 0, &f0_s_rden);
1211: DMPlexComputeIntegralFEM(plex,X,tt,ctx);
1212: density[ii] = tt[0]*ctx->n_0*ctx->charges[ii];
1213: PetscDSSetObjective(prob, 0, &f0_s_rmom);
1214: DMPlexComputeIntegralFEM(plex,X,tt,ctx);
1215: zmomentum[ii] = tt[0]*ctx->n_0*ctx->v_0*ctx->masses[ii];
1216: PetscDSSetObjective(prob, 0, &f0_s_rv2);
1217: DMPlexComputeIntegralFEM(plex,X,tt,ctx);
1218: energy[ii] = tt[0]*0.5*ctx->n_0*ctx->v_0*ctx->v_0*ctx->masses[ii];
1219: zmomentumtot += zmomentum[ii];
1220: energytot += energy[ii];
1221: densitytot += density[ii];
1222: PetscPrintf(PETSC_COMM_WORLD, "%3D) species-%D: charge density= %20.13e z-momentum= %20.13e energy= %20.13e",stepi,ii,PetscRealPart(density[ii]),PetscRealPart(zmomentum[ii]),PetscRealPart(energy[ii]));
1223: } else { /* 2/3X + 3V */
1224: PetscDSSetObjective(prob, 0, &f0_s_den);
1225: DMPlexComputeIntegralFEM(plex,X,tt,ctx);
1226: density[ii] = tt[0]*ctx->n_0*ctx->charges[ii];
1227: PetscDSSetObjective(prob, 0, &f0_s_mom);
1228: user[1] = 0;
1229: DMPlexComputeIntegralFEM(plex,X,tt,ctx);
1230: xmomentum[ii] = tt[0]*ctx->n_0*ctx->v_0*ctx->masses[ii];
1231: user[1] = 1;
1232: DMPlexComputeIntegralFEM(plex,X,tt,ctx);
1233: ymomentum[ii] = tt[0]*ctx->n_0*ctx->v_0*ctx->masses[ii];
1234: user[1] = 2;
1235: DMPlexComputeIntegralFEM(plex,X,tt,ctx);
1236: zmomentum[ii] = tt[0]*ctx->n_0*ctx->v_0*ctx->masses[ii];
1237: PetscDSSetObjective(prob, 0, &f0_s_v2);
1238: DMPlexComputeIntegralFEM(plex,X,tt,ctx);
1239: energy[ii] = 0.5*tt[0]*ctx->n_0*ctx->v_0*ctx->v_0*ctx->masses[ii];
1240: PetscPrintf(PETSC_COMM_WORLD, "%3D) species %D: density=%20.13e, x-momentum=%20.13e, y-momentum=%20.13e, z-momentum=%20.13e, energy=%21.13e",
1241: stepi,ii,PetscRealPart(density[ii]),PetscRealPart(xmomentum[ii]),PetscRealPart(ymomentum[ii]),PetscRealPart(zmomentum[ii]),PetscRealPart(energy[ii]));
1242: xmomentumtot += xmomentum[ii];
1243: ymomentumtot += ymomentum[ii];
1244: zmomentumtot += zmomentum[ii];
1245: energytot += energy[ii];
1246: densitytot += density[ii];
1247: }
1248: if (ctx->num_species>1) PetscPrintf(PETSC_COMM_WORLD, "\n");
1249: }
1250: /* totals */
1251: DMPlexGetHeightStratum(plex,0,&cStart,&cEnd);
1252: DMDestroy(&plex);
1253: if (ctx->num_species>1) {
1254: if (dim==2) {
1255: PetscPrintf(PETSC_COMM_WORLD, "\t%3D) Total: charge density=%21.13e, momentum=%21.13e, energy=%21.13e (m_i[0]/m_e = %g, %D cells)",
1256: stepi,PetscRealPart(densitytot),PetscRealPart(zmomentumtot),PetscRealPart(energytot),ctx->masses[1]/ctx->masses[0],cEnd-cStart);
1257: } else {
1258: PetscPrintf(PETSC_COMM_WORLD, "\t%3D) Total: charge density=%21.13e, x-momentum=%21.13e, y-momentum=%21.13e, z-momentum=%21.13e, energy=%21.13e (m_i[0]/m_e = %g, %D cells)",
1259: stepi,PetscRealPart(densitytot),PetscRealPart(xmomentumtot),PetscRealPart(ymomentumtot),PetscRealPart(zmomentumtot),PetscRealPart(energytot),ctx->masses[1]/ctx->masses[0],cEnd-cStart);
1260: }
1261: } else {
1262: PetscPrintf(PETSC_COMM_WORLD, " -- %D cells",cEnd-cStart);
1263: }
1264: if (ctx->verbose > 1) {PetscPrintf(PETSC_COMM_WORLD,", %D sub (vector) threads\n",ctx->subThreadBlockSize);}
1265: else {PetscPrintf(PETSC_COMM_WORLD,"\n");}
1266: return(0);
1267: }
1269: static PetscErrorCode destroy_coloring (void *is)
1270: {
1271: ISColoring tmp = (ISColoring)is;
1272: return ISColoringDestroy(&tmp);
1273: }
1275: /*@
1276: LandauCreateColoring - create a coloring and add to matrix (Landau context used just for 'print' flag, should be in DMPlex)
1278: Collective on JacP
1280: Input Parameters:
1281: + JacP - matrix to add coloring to
1282: - plex - The DM
1284: Output Parameter:
1285: . container - Container with coloring
1287: Level: beginner
1289: .keywords: mesh
1290: .seealso: LandauCreateVelocitySpace()
1291: @*/
1292: PetscErrorCode LandauCreateColoring(Mat JacP, DM plex, PetscContainer *container)
1293: {
1294: PetscErrorCode ierr;
1295: PetscInt dim,cell,i,ej,nc,Nv,totDim,numGCells,cStart,cEnd;
1296: ISColoring iscoloring = NULL;
1297: Mat G,Q;
1298: PetscScalar ones[128];
1299: MatColoring mc;
1300: IS *is;
1301: PetscInt csize,colour,j,k;
1302: const PetscInt *indices;
1303: PetscInt numComp[1];
1304: PetscInt numDof[4];
1305: PetscFE fe;
1306: DM colordm;
1307: PetscSection csection, section, globalSection;
1308: PetscDS prob;
1309: LandauCtx *ctx;
1312: DMGetApplicationContext(plex, &ctx);
1313: DMGetLocalSection(plex, §ion);
1314: DMGetGlobalSection(plex, &globalSection);
1315: DMGetDimension(plex, &dim);
1316: DMGetDS(plex, &prob);
1317: PetscDSGetTotalDimension(prob, &totDim);
1318: DMPlexGetHeightStratum(plex,0,&cStart,&cEnd);
1319: numGCells = cEnd - cStart;
1320: /* create cell centered DM */
1321: DMClone(plex, &colordm);
1322: PetscFECreateDefault(PetscObjectComm((PetscObject) plex), dim, 1, PETSC_FALSE, "color_", PETSC_DECIDE, &fe);
1323: PetscObjectSetName((PetscObject) fe, "color");
1324: DMSetField(colordm, 0, NULL, (PetscObject)fe);
1325: PetscFEDestroy(&fe);
1326: for (i = 0; i < (dim+1); ++i) numDof[i] = 0;
1327: numDof[dim] = 1;
1328: numComp[0] = 1;
1329: DMPlexCreateSection(colordm, NULL, numComp, numDof, 0, NULL, NULL, NULL, NULL, &csection);
1330: PetscSectionSetFieldName(csection, 0, "color");
1331: DMSetLocalSection(colordm, csection);
1332: DMViewFromOptions(colordm,NULL,"-color_dm_view");
1333: /* get vertex to element map Q and colroing graph G */
1334: MatGetSize(JacP,NULL,&Nv);
1335: MatCreateAIJ(PETSC_COMM_SELF,PETSC_DECIDE,PETSC_DECIDE,numGCells,Nv,totDim,NULL,0,NULL,&Q);
1336: for (i=0;i<128;i++) ones[i] = 1.0;
1337: for (cell = cStart, ej = 0 ; cell < cEnd; ++cell, ++ej) {
1338: PetscInt numindices,*indices;
1339: DMPlexGetClosureIndices(plex, section, globalSection, cell, PETSC_TRUE, &numindices, &indices, NULL, NULL);
1340: if (numindices>128) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "too many indices. %D > %D",numindices,128);
1341: MatSetValues(Q,1,&ej,numindices,indices,ones,ADD_VALUES);
1342: DMPlexRestoreClosureIndices(plex, section, globalSection, cell, PETSC_TRUE, &numindices, &indices, NULL, NULL);
1343: }
1344: MatAssemblyBegin(Q, MAT_FINAL_ASSEMBLY);
1345: MatAssemblyEnd(Q, MAT_FINAL_ASSEMBLY);
1346: MatMatTransposeMult(Q,Q,MAT_INITIAL_MATRIX,4.0,&G);
1347: PetscObjectSetName((PetscObject) Q, "Q");
1348: PetscObjectSetName((PetscObject) G, "coloring graph");
1349: MatViewFromOptions(G,NULL,"-coloring_mat_view");
1350: MatViewFromOptions(Q,NULL,"-coloring_mat_view");
1351: MatDestroy(&Q);
1352: /* coloring */
1353: MatColoringCreate(G,&mc);
1354: MatColoringSetDistance(mc,1);
1355: MatColoringSetType(mc,MATCOLORINGJP);
1356: MatColoringSetFromOptions(mc);
1357: MatColoringApply(mc,&iscoloring);
1358: MatColoringDestroy(&mc);
1359: /* view */
1360: ISColoringViewFromOptions(iscoloring,NULL,"-coloring_is_view");
1361: ISColoringGetIS(iscoloring,PETSC_USE_POINTER,&nc,&is);
1362: if (ctx && ctx->verbose > 5) {
1363: PetscViewer viewer;
1364: Vec color_vec, eidx_vec;
1365: DMGetGlobalVector(colordm, &color_vec);
1366: DMGetGlobalVector(colordm, &eidx_vec);
1367: for (colour=0; colour<nc; colour++) {
1368: ISGetLocalSize(is[colour],&csize);
1369: ISGetIndices(is[colour],&indices);
1370: for (j=0; j<csize; j++) {
1371: PetscScalar v = (PetscScalar)colour;
1372: k = indices[j];
1373: VecSetValues(color_vec,1,&k,&v,INSERT_VALUES);
1374: v = (PetscScalar)k;
1375: VecSetValues(eidx_vec,1,&k,&v,INSERT_VALUES);
1376: }
1377: ISRestoreIndices(is[colour],&indices);
1378: }
1379: /* view */
1380: PetscViewerVTKOpen(PETSC_COMM_WORLD, "color.vtu", FILE_MODE_WRITE, &viewer);
1381: PetscObjectSetName((PetscObject) color_vec, "color");
1382: VecView(color_vec, viewer);
1383: PetscViewerDestroy(&viewer);
1384: PetscViewerVTKOpen(PETSC_COMM_WORLD, "eidx.vtu", FILE_MODE_WRITE, &viewer);
1385: PetscObjectSetName((PetscObject) eidx_vec, "element-idx");
1386: VecView(eidx_vec, viewer);
1387: PetscViewerDestroy(&viewer);
1388: DMRestoreGlobalVector(colordm, &color_vec);
1389: DMRestoreGlobalVector(colordm, &eidx_vec);
1390: }
1391: PetscSectionDestroy(&csection);
1392: DMDestroy(&colordm);
1393: ISColoringRestoreIS(iscoloring,PETSC_USE_POINTER,&is);
1394: MatDestroy(&G);
1395: /* stash coloring */
1396: PetscContainerCreate(PETSC_COMM_SELF, container);
1397: PetscContainerSetPointer(*container,(void*)iscoloring);
1398: PetscContainerSetUserDestroy(*container, destroy_coloring);
1399: PetscObjectCompose((PetscObject)JacP,"coloring",(PetscObject)*container);
1400: if (ctx && ctx->verbose > 0) {
1401: PetscPrintf(PETSC_COMM_WORLD, "Made coloring with %D colors\n", nc);
1402: }
1403: return(0);
1404: }
1406: PetscErrorCode LandauAssembleOpenMP(PetscInt cStart, PetscInt cEnd, PetscInt totDim, DM plex, PetscSection section, PetscSection globalSection, Mat JacP, PetscScalar elemMats[], PetscContainer container)
1407: {
1408: PetscErrorCode ierr;
1409: IS *is;
1410: PetscInt nc,colour,j;
1411: const PetscInt *clr_idxs;
1412: ISColoring iscoloring;
1414: PetscContainerGetPointer(container,(void**)&iscoloring);
1415: ISColoringGetIS(iscoloring,PETSC_USE_POINTER,&nc,&is);
1416: for (colour=0; colour<nc; colour++) {
1417: PetscInt *idx_arr[1024]; /* need to make dynamic for general use */
1418: PetscScalar *new_el_mats[1024];
1419: PetscInt idx_size[1024],csize;
1420: ISGetLocalSize(is[colour],&csize);
1421: if (csize>1024) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "too many elements in color. %D > %D",csize,1024);
1422: ISGetIndices(is[colour],&clr_idxs);
1423: /* get indices and mats */
1424: for (j=0; j<csize; j++) {
1425: PetscInt cell = cStart + clr_idxs[j];
1426: PetscInt numindices,*indices;
1427: PetscScalar *elMat = &elemMats[clr_idxs[j]*totDim*totDim];
1428: PetscScalar *valuesOrig = elMat;
1429: DMPlexGetClosureIndices(plex, section, globalSection, cell, PETSC_TRUE, &numindices, &indices, NULL, (PetscScalar **) &elMat);
1430: idx_size[j] = numindices;
1431: PetscMalloc2(numindices,&idx_arr[j],numindices*numindices,&new_el_mats[j]);
1432: PetscMemcpy(idx_arr[j],indices,numindices*sizeof(PetscInt));
1433: PetscMemcpy(new_el_mats[j],elMat,numindices*numindices*sizeof(PetscScalar));
1434: DMPlexRestoreClosureIndices(plex, section, globalSection, cell, PETSC_TRUE, &numindices, &indices, NULL, (PetscScalar **) &elMat);
1435: if (elMat != valuesOrig) {DMRestoreWorkArray(plex, numindices*numindices, MPIU_SCALAR, &elMat);}
1436: }
1437: /* assemble matrix - pragmas break CI ? */
1438: //#pragma omp parallel default(JacP,idx_size,idx_arr,new_el_mats,colour,clr_idxs) private(j)
1439: //#pragma omp parallel for private(j)
1440: for (j=0; j<csize; j++) {
1441: PetscInt numindices = idx_size[j], *indices = idx_arr[j];
1442: PetscScalar *elMat = new_el_mats[j];
1443: MatSetValues(JacP,numindices,indices,numindices,indices,elMat,ADD_VALUES);
1444: }
1445: /* free */
1446: ISRestoreIndices(is[colour],&clr_idxs);
1447: for (j=0; j<csize; j++) {
1448: PetscFree2(idx_arr[j],new_el_mats[j]);
1449: }
1450: }
1451: ISColoringRestoreIS(iscoloring,PETSC_USE_POINTER,&is);
1452: return(0);
1453: }
1455: /* < v, u > */
1456: static void g0_1(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1457: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1458: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1459: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
1460: {
1461: g0[0] = 1.;
1462: }
1464: /* < v, u > */
1465: static void g0_r(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1466: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1467: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1468: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
1469: {
1470: g0[0] = 2.*PETSC_PI*x[0];
1471: }
1473: /*@
1474: LandauCreateMassMatrix - Create mass matrix for Landau
1476: Collective on dm
1478: Input Parameters:
1479: . dm - the DM object
1481: Output Parameters:
1482: . Amat - The mass matrix (optional), mass matrix is added to the DM context
1484: Level: beginner
1486: .keywords: mesh
1487: .seealso: LandauCreateVelocitySpace()
1488: @*/
1489: PetscErrorCode LandauCreateMassMatrix(DM dm, Mat *Amat)
1490: {
1491: DM massDM;
1492: PetscDS prob;
1493: PetscInt ii,dim,N1=1,N2;
1495: LandauCtx *ctx;
1496: Mat M;
1501: DMGetApplicationContext(dm, &ctx);
1502: if (!ctx) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "no context");
1503: DMGetDimension(dm, &dim);
1504: DMClone(dm, &massDM);
1505: DMCopyFields(dm, massDM);
1506: DMCreateDS(massDM);
1507: DMGetDS(massDM, &prob);
1508: for (ii=0;ii<ctx->num_species;ii++) {
1509: if (dim==3) {PetscDSSetJacobian(prob, ii, ii, g0_1, NULL, NULL, NULL);}
1510: else {PetscDSSetJacobian(prob, ii, ii, g0_r, NULL, NULL, NULL);}
1511: }
1512: DMViewFromOptions(massDM,NULL,"-dm_landau_mass_dm_view");
1513: DMCreateMatrix(massDM, &M);
1514: {
1515: Vec locX;
1516: DM plex;
1517: DMConvert(massDM, DMPLEX, &plex);
1518: DMGetLocalVector(massDM, &locX);
1519: /* Mass matrix is independent of the input, so no need to fill locX */
1520: DMPlexSNESComputeJacobianFEM(plex, locX, M, M, ctx);
1521: DMRestoreLocalVector(massDM, &locX);
1522: DMDestroy(&plex);
1523: }
1524: DMDestroy(&massDM);
1525: MatGetSize(ctx->J, &N1, NULL);
1526: MatGetSize(M, &N2, NULL);
1527: if (N1 != N2) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "Incorrect matrix sizes: |Jacobian| = %D, |Mass|=%D",N1,N2);
1528: MatViewFromOptions(M,NULL,"-dm_landau_mass_mat_view");
1529: ctx->M = M; /* this could be a noop, a = a */
1530: if (Amat) *Amat = M;
1531: return(0);
1532: }
1534: /*@
1535: LandauIFunction - TS residual calculation
1537: Collective on ts
1539: Input Parameters:
1540: + TS - The time stepping context
1541: . time_dummy - current time (not used)
1542: - X - Current state
1543: + X_t - Time derivative of current state
1544: . actx - Landau context
1546: Output Parameter:
1547: . F - The residual
1549: Level: beginner
1551: .keywords: mesh
1552: .seealso: LandauCreateVelocitySpace(), LandauIJacobian()
1553: @*/
1554: PetscErrorCode LandauIFunction(TS ts, PetscReal time_dummy, Vec X, Vec X_t, Vec F, void *actx)
1555: {
1557: LandauCtx *ctx=(LandauCtx*)actx;
1558: PetscReal unorm;
1559: PetscInt dim;
1560: DM dm;
1563: TSGetDM(ts,&dm);
1564: DMGetApplicationContext(dm, &ctx);
1565: if (!ctx) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "no context");
1566: VecNorm(X,NORM_2,&unorm);
1567: MPI_Barrier(PETSC_COMM_WORLD); // remove in real application
1568: PetscLogEventBegin(ctx->events[0],0,0,0,0);
1569: DMGetDimension(ctx->dmv, &dim);
1570: if (ctx->normJ != unorm) {
1571: ctx->normJ = unorm;
1572: PetscInfo1(ts, "Create Landau Jacobian t=%g\n",time_dummy);
1573: LandauFormJacobian_Internal(X,ctx->J,dim,(void*)ctx);
1574: ctx->aux_bool = PETSC_TRUE; /* debug: set flag that we made a new Jacobian */
1575: } else {
1576: ctx->aux_bool = PETSC_FALSE;
1577: PetscInfo1(ts, "Skip Landau Jacobian t=%g\n",(double)time_dummy);
1578: }
1579: /* mat vec for op */
1580: MatMult(ctx->J,X,F); /* C*f */
1581: /* add time term */
1582: if (X_t) {
1583: MatMultAdd(ctx->M,X_t,F,F);
1584: }
1585: PetscLogEventEnd(ctx->events[0],0,0,0,0);
1586: return(0);
1587: }
1589: /*@
1590: LandauIJacobian - TS Jacobian construction
1592: Collective on ts
1594: Input Parameters:
1595: + TS - The time stepping context
1596: . time_dummy - current time (not used)
1597: - X - Current state
1598: + U_tdummy - Time derivative of current state (not used)
1599: . shift - shift for du/dt term
1600: - actx - Landau context
1602: Output Parameter:
1603: . Amat - Jacobian
1604: + Pmat - same as Amat
1606: Level: beginner
1608: .keywords: mesh
1609: .seealso: LandauCreateVelocitySpace(), LandauIFunction()
1610: @*/
1611: PetscErrorCode LandauIJacobian(TS ts, PetscReal time_dummy, Vec X, Vec U_tdummy, PetscReal shift, Mat Amat, Mat Pmat, void *actx)
1612: {
1614: LandauCtx *ctx=(LandauCtx*)actx;
1615: PetscReal unorm;
1616: PetscInt dim;
1617: DM dm;
1620: TSGetDM(ts,&dm);
1621: DMGetApplicationContext(dm, &ctx);
1622: if (!ctx) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "no context");
1623: if (Amat!=Pmat || Amat!=ctx->J) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Amat!=Pmat || Amat!=ctx->J");
1624: DMGetDimension(ctx->dmv, &dim);
1625: /* get collision Jacobian into A */
1626: PetscLogEventBegin(ctx->events[9],0,0,0,0);
1627: VecNorm(X,NORM_2,&unorm);
1628: if (ctx->normJ!=unorm) {
1629: PetscInfo2(ts, "Create Landau Jacobian t=%g, shift=%g\n",(double)time_dummy,(double)shift);
1630: LandauFormJacobian_Internal(X,ctx->J,dim,(void*)ctx);
1631: ctx->normJ = unorm;
1632: ctx->aux_bool = PETSC_TRUE; /* debug: set flag that we made a new Jacobian */
1633: } else {
1634: ctx->aux_bool = PETSC_FALSE;
1635: PetscInfo3(ts, "Skip Landau Jacobian t=%g, shift=%g shift*|u|=%20.12e\n",(double)time_dummy,(double)shift,(double)shift*unorm);
1636: }
1637: /* add C */
1638: MatCopy(ctx->J,Pmat,SAME_NONZERO_PATTERN);
1639: /* add mass */
1640: MatAXPY(Pmat,shift,ctx->M,SAME_NONZERO_PATTERN);
1641: PetscLogEventEnd(ctx->events[9],0,0,0,0);
1642: return(0);
1643: }