Actual source code: ex48.c
petsc-3.11.1 2019-04-12
1: static char help[] = " Evolution of magnetic islands.\n\
2: The aim of this model is to self-consistently study the interaction between the tearing mode and small scale drift-wave turbulence.\n\n\n" ;
This is a three field model for the density $\tilde n$, vorticity $\tilde\Omega$, and magnetic flux $\tilde\psi$, using auxiliary variables potential $\tilde\phi$ and current $j_z$.
\begin{equation}
\begin{aligned}
\partial_t \tilde n &= \left\{ \tilde n, \tilde\phi \right\} + \beta \left\{ j_z, \tilde\psi \right\} + \left\{ \ln n_0, \tilde\phi \right\} + \mu \nabla^2_\perp \tilde n \\
\partial_t \tilde\Omega &= \left\{ \tilde\Omega, \tilde\phi \right\} + \beta \left\{ j_z, \tilde\psi \right\} + \mu \nabla^2_\perp \tilde\Omega \\
\partial_t \tilde\psi &= \left\{ \psi_0 + \tilde\psi, \tilde\phi - \tilde n \right\} - \left\{ \ln n_0, \tilde\psi \right\} + \frac{\eta}{\beta} \nabla^2_\perp \tilde\psi \\
\nabla^2_\perp\tilde\phi &= \tilde\Omega \\
j_z &= -\nabla^2_\perp \left( \tilde\psi + \psi_0 \right)\\
\end{aligned}
\end{equation}
17: #include <petscdmplex.h>
18: #include <petscts.h>
19: #include <petscds.h>
20: #include <assert.h>
22: typedef struct {
23: PetscInt debug; /* The debugging level */
24: PetscBool plotRef; /* Plot the reference fields */
25: /* Domain and mesh definition */
26: PetscInt dim; /* The topological mesh dimension */
27: char filename[2048]; /* The optional ExodusII file */
28: PetscBool cell_simplex; /* Simplicial mesh */
29: DMBoundaryType boundary_types[3];
30: PetscInt cells[3];
31: PetscInt refine;
32: /* geometry */
33: PetscReal domain_lo[3], domain_hi[3];
34: DMBoundaryType periodicity[3]; /* The domain periodicity */
35: PetscReal b0[3]; /* not used */
36: /* Problem definition */
37: PetscErrorCode ( **initialFuncs)( PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
38: PetscReal mu, eta, beta;
39: PetscReal a,b,Jo,Jop,m,ke,kx,ky,DeltaPrime,eps;
40: /* solver */
41: PetscBool implicit;
42: } AppCtx;
44: static AppCtx *s_ctx;
46: static PetscScalar poissonBracket( PetscInt dim, const PetscScalar df[], const PetscScalar dg[])
47: {
48: PetscScalar ret = df[0]*dg[1] - df[1]*dg[0];
49: return ret;
50: }
52: enum field_idx {DENSITY,OMEGA,PSI,PHI,JZ};
54: static void f0_n( PetscInt dim, PetscInt Nf, PetscInt NfAux,
55: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
56: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
57: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
58: {
59: const PetscScalar *pnDer = &u_x[uOff_x[DENSITY]];
60: const PetscScalar *ppsiDer = &u_x[uOff_x[PSI]];
61: const PetscScalar *pphiDer = &u_x[uOff_x[PHI]];
62: const PetscScalar *jzDer = &u_x[uOff_x[JZ]];
63: const PetscScalar *logRefDenDer = &a_x[aOff_x[DENSITY]];
64: f0[0] += - poissonBracket( dim,pnDer, pphiDer) - s_ctx->beta*poissonBracket( dim,jzDer, ppsiDer) - poissonBracket( dim,logRefDenDer, pphiDer);
65: if ( u_t) f0[0] += u_t[DENSITY];
66: }
68: static void f1_n( PetscInt dim, PetscInt Nf, PetscInt NfAux,
69: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
70: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
71: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
72: {
73: const PetscScalar *pnDer = &u_x[uOff_x[DENSITY]];
74: PetscInt d;
76: for ( d = 0; d < dim-1; ++d) f1[d] = -s_ctx->mu*pnDer[d];
77: }
79: static void f0_Omega( PetscInt dim, PetscInt Nf, PetscInt NfAux,
80: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
81: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
82: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
83: {
84: const PetscScalar *pOmegaDer = &u_x[uOff_x[OMEGA]];
85: const PetscScalar *ppsiDer = &u_x[uOff_x[PSI]];
86: const PetscScalar *pphiDer = &u_x[uOff_x[PHI]];
87: const PetscScalar *jzDer = &u_x[uOff_x[JZ]];
89: f0[0] += - poissonBracket( dim,pOmegaDer, pphiDer) - s_ctx->beta*poissonBracket( dim,jzDer, ppsiDer);
90: if ( u_t) f0[0] += u_t[OMEGA];
91: }
93: static void f1_Omega( PetscInt dim, PetscInt Nf, PetscInt NfAux,
94: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
95: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
96: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
97: {
98: const PetscScalar *pOmegaDer = &u_x[uOff_x[OMEGA]];
99: PetscInt d;
101: for ( d = 0; d < dim-1; ++d) f1[d] = -s_ctx->mu*pOmegaDer[d];
102: }
104: static void f0_psi( PetscInt dim, PetscInt Nf, PetscInt NfAux,
105: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
106: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
107: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
108: {
109: const PetscScalar *pnDer = &u_x[uOff_x[DENSITY]];
110: const PetscScalar *ppsiDer = &u_x[uOff_x[PSI]];
111: const PetscScalar *pphiDer = &u_x[uOff_x[PHI]];
112: const PetscScalar *refPsiDer = &a_x[aOff_x[PSI]];
113: const PetscScalar *logRefDenDer= &a_x[aOff_x[DENSITY]];
114: PetscScalar psiDer[3];
115: PetscScalar phi_n_Der[3];
116: PetscInt d;
117: if ( dim < 2) {MPI_Abort ( MPI_COMM_WORLD,1); return ;} /* this is needed so that the clang static analyzer does not generate a warning about variables used by not set */
118: for ( d = 0; d < dim; ++d) {
119: psiDer[d] = refPsiDer[d] + ppsiDer[d];
120: phi_n_Der[d] = pphiDer[d] - pnDer[d];
121: }
122: f0[0] = - poissonBracket( dim,psiDer, phi_n_Der) + poissonBracket( dim,logRefDenDer, ppsiDer);
123: if ( u_t) f0[0] += u_t[PSI];
124: }
126: static void f1_psi( PetscInt dim, PetscInt Nf, PetscInt NfAux,
127: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
128: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
129: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
130: {
131: const PetscScalar *ppsi = &u_x[uOff_x[PSI]];
132: PetscInt d;
134: for ( d = 0; d < dim-1; ++d) f1[d] = -( s_ctx->eta/s_ctx->beta)*ppsi[d];
135: }
137: static void f0_phi( PetscInt dim, PetscInt Nf, PetscInt NfAux,
138: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
139: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
140: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
141: {
142: f0[0] = -u[uOff[OMEGA]];
143: }
145: static void f1_phi( PetscInt dim, PetscInt Nf, PetscInt NfAux,
146: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
147: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
148: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
149: {
150: const PetscScalar *pphi = &u_x[uOff_x[PHI]];
151: PetscInt d;
153: for ( d = 0; d < dim-1; ++d) f1[d] = pphi[d];
154: }
156: static void f0_jz( PetscInt dim, PetscInt Nf, PetscInt NfAux,
157: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
158: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
159: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
160: {
161: f0[0] = u[uOff[JZ]];
162: }
164: static void f1_jz( PetscInt dim, PetscInt Nf, PetscInt NfAux,
165: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
166: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
167: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
168: {
169: const PetscScalar *ppsi = &u_x[uOff_x[PSI]];
170: const PetscScalar *refPsiDer = &a_x[aOff_x[PSI]]; /* aOff_x[PSI] = = 2*PSI */
171: PetscInt d;
173: for ( d = 0; d < dim-1; ++d) f1[d] = ppsi[d] + refPsiDer[d];
174: }
176: static PetscErrorCode ProcessOptions( MPI_Comm comm, AppCtx *options)
177: {
178: PetscBool flg;
180: PetscInt ii, bd;
182: options->debug = 1;
183: options->plotRef = PETSC_FALSE ;
184: options->dim = 2;
185: options->filename[0] = '\0';
186: options->cell_simplex = PETSC_FALSE ;
187: options->implicit = PETSC_FALSE ;
188: options->refine = 2;
189: options->domain_lo[0] = 0.0;
190: options->domain_lo[1] = 0.0;
191: options->domain_lo[2] = 0.0;
192: options->domain_hi[0] = 2.0;
193: options->domain_hi[1] = 2.0*PETSC_PI;
194: options->domain_hi[2] = 2.0;
195: options->periodicity[0] = DM_BOUNDARY_NONE ;
196: options->periodicity[1] = DM_BOUNDARY_NONE ;
197: options->periodicity[2] = DM_BOUNDARY_NONE ;
198: options->mu = 0;
199: options->eta = 0;
200: options->beta = 1;
201: options->a = 1;
202: options->b = PETSC_PI;
203: options->Jop = 0;
204: options->m = 1;
205: options->eps = 1.e-6;
207: for ( ii = 0; ii < options->dim; ++ii) options->cells[ii] = 4;
208: PetscOptionsBegin ( comm, " " , " Poisson Problem Options" , " DMPLEX " );
209: PetscOptionsInt ( " -debug" , " The debugging level" , " ex48.c" , options->debug, &options->debug, NULL);
210: PetscOptionsBool ( " -plot_ref" , " Plot the reference fields" , " ex48.c" , options->plotRef, &options->plotRef, NULL);
211: PetscOptionsInt ( " -dim" , " The topological mesh dimension" , " ex48.c" , options->dim, &options->dim, NULL);
212: if ( options->dim < 2 || options->dim > 3) SETERRQ1 ( PETSC_COMM_WORLD ,PETSC_ERR_ARG_OUTOFRANGE," Dim %D must be 2 or 3" ,options->dim);
213: PetscOptionsInt ( " -dm_refine" , " Hack to get refinement level for cylinder" , " ex48.c" , options->refine, &options->refine, NULL);
214: PetscOptionsReal ( " -mu" , " mu" , " ex48.c" , options->mu, &options->mu, NULL);
215: PetscOptionsReal ( " -eta" , " eta" , " ex48.c" , options->eta, &options->eta, NULL);
216: PetscOptionsReal ( " -beta" , " beta" , " ex48.c" , options->beta, &options->beta, NULL);
217: PetscOptionsReal ( " -Jop" , " Jop" , " ex48.c" , options->Jop, &options->Jop, NULL);
218: PetscOptionsReal ( " -m" , " m" , " ex48.c" , options->m, &options->m, NULL);
219: PetscOptionsReal ( " -eps" , " eps" , " ex48.c" , options->eps, &options->eps, NULL);
220: PetscOptionsString ( " -f" , " Exodus.II filename to read" , " ex48.c" , options->filename, options->filename, sizeof ( options->filename), &flg);
221: PetscOptionsBool ( " -cell_simplex" , " Simplicial ( true) or tensor ( false) mesh" , " ex48.c" , options->cell_simplex, &options->cell_simplex, NULL);
222: PetscOptionsBool ( " -implicit" , " Use implicit time integrator" , " ex48.c" , options->implicit, &options->implicit, NULL);
223: ii = options->dim;
224: PetscOptionsRealArray ( " -domain_hi" , " Domain size" , " ex48.c" , options->domain_hi, &ii, NULL);
225: ii = options->dim;
226: PetscOptionsRealArray ( " -domain_lo" , " Domain size" , " ex48.c" , options->domain_lo, &ii, NULL);
227: ii = options->dim;
228: bd = options->periodicity[0];
229: PetscOptionsEList ( " -x_periodicity" , " The x-boundary periodicity" , " ex48.c" , DMBoundaryTypes, 5, DMBoundaryTypes[options->periodicity[0]], &bd, NULL);
230: options->periodicity[0] = ( DMBoundaryType ) bd;
231: bd = options->periodicity[1];
232: PetscOptionsEList ( " -y_periodicity" , " The y-boundary periodicity" , " ex48.c" , DMBoundaryTypes, 5, DMBoundaryTypes[options->periodicity[1]], &bd, NULL);
233: options->periodicity[1] = ( DMBoundaryType ) bd;
234: bd = options->periodicity[2];
235: PetscOptionsEList ( " -z_periodicity" , " The z-boundary periodicity" , " ex48.c" , DMBoundaryTypes, 5, DMBoundaryTypes[options->periodicity[2]], &bd, NULL);
236: options->periodicity[2] = ( DMBoundaryType ) bd;
237: ii = options->dim;
238: PetscOptionsIntArray ( " -cells" , " Number of cells in each dimension" , " ex48.c" , options->cells, &ii, NULL);
239: PetscOptionsEnd ( );
240: options->a = ( options->domain_hi[0]-options->domain_lo[0])/2.0;
241: options->b = ( options->domain_hi[1]-options->domain_lo[1])/2.0;
242: for ( ii = 0; ii < options->dim; ++ii) {
243: if ( options->domain_hi[ii] <= options->domain_lo[ii]) SETERRQ3 ( comm,PETSC_ERR_ARG_WRONG," Domain %D lo= %g hi= %g" ,ii,options->domain_lo[ii],options->domain_hi[ii]);
244: }
245: options->ke = PetscSqrtScalar( options->Jop);
246: if ( options->Jop= = 0.0) {
247: options->Jo = 1.0/PetscPowScalar( options->a,2);
248: } else {
249: options->Jo = options->Jop*PetscCosReal( options->ke*options->a)/( 1.0-PetscCosReal( options->ke*options->a));
250: }
251: options->ky = PETSC_PI*options->m/options->b;
252: if ( PetscPowReal( options->ky, 2) < options->Jop) {
253: options->kx = PetscSqrtScalar( options->Jop-PetscPowScalar( options->ky,2));
254: options->DeltaPrime = -2.0*options->kx*options->a*PetscCosReal( options->kx*options->a)/PetscSinReal( options->kx*options->a);
255: } else if ( PetscPowReal( options->ky, 2) > options->Jop) {
256: options->kx = PetscSqrtScalar( PetscPowScalar( options->ky,2)-options->Jop);
257: options->DeltaPrime = -2.0*options->kx*options->a*PetscCoshReal( options->kx*options->a)/PetscSinhReal( options->kx*options->a);
258: } else { /*they're equal ( or there's a NaN), lim( x*cot( x))_x->0= 1*/
259: options->kx = 0;
260: options->DeltaPrime = -2.0;
261: }
262: PetscPrintf ( comm, " DeltaPrime= %g\n" ,options->DeltaPrime);
264: return ( 0);
265: }
269: static void f_n( PetscInt dim, PetscInt Nf, PetscInt NfAux,
270: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
271: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
272: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar *f0)
273: {
274: const PetscScalar *pn = &u[uOff[DENSITY]];
275: *f0 = *pn;
276: }
280: static PetscErrorCode PostStep( TS ts)
281: {
282: PetscErrorCode ierr;
283: DM dm;
284: AppCtx *ctx;
285: PetscInt stepi,num;
286: Vec X;
288: TSGetApplicationContext ( ts, &ctx); assert( ctx);
289: if ( ctx->debug<1) return ( 0);
290: TSGetSolution ( ts, &X);
291: VecGetDM ( X, &dm);
292: TSGetStepNumber ( ts, &stepi);
293: DMGetOutputSequenceNumber ( dm, &num, NULL);
294: if ( num < 0) {DMSetOutputSequenceNumber ( dm, 0, 0.0);}
295: PetscObjectSetName ( ( PetscObject ) X, " u" );
296: VecViewFromOptions( X, NULL, " -vec_view" );
297: /* print integrals */
298: {
299: PetscDS prob;
300: DM plex;
301: PetscScalar den, tt[5];
302: DMConvert ( dm, DMPLEX , &plex);
303: DMGetDS ( plex, &prob);
304: PetscDSSetObjective( prob, 0, &f_n);
305: DMPlexComputeIntegralFEM ( plex,X,tt,ctx);
306: den = tt[0];
307: DMDestroy ( &plex);
308: PetscPrintf ( PetscObjectComm ( ( PetscObject )dm), " %D) total perturbed mass = %g\n" , stepi, ( double) PetscRealPart ( den));
309: }
310: return ( 0);
311: }
313: static PetscErrorCode CreateBCLabel( DM dm, const char name[])
314: {
315: DMLabel label;
318: DMCreateLabel ( dm, name);
319: DMGetLabel ( dm, name, &label);
320: DMPlexMarkBoundaryFaces ( dm, 1, label);
321: DMPlexLabelComplete ( dm, label);
322: return ( 0);
323: }
325: static PetscErrorCode CreateMesh( MPI_Comm comm, AppCtx *ctx, DM *dm)
326: {
327: PetscInt dim = ctx->dim;
328: const char *filename = ctx->filename;
329: size_t len;
330: PetscMPIInt numProcs;
334: MPI_Comm_size ( comm, &numProcs);
335: PetscStrlen ( filename, &len);
336: if ( len) {
337: DMPlexCreateFromFile ( comm, filename, PETSC_TRUE , dm);
338: } else {
339: PetscInt d;
341: /* create DM */
342: if ( ctx->cell_simplex && dim = = 3) SETERRQ ( comm, PETSC_ERR_ARG_WRONG, " Cannot mesh a cylinder with simplices" );
343: if ( dim= = 2) {
344: PetscInt refineRatio, totCells = 1;
345: if ( ctx->cell_simplex) SETERRQ ( comm, PETSC_ERR_ARG_WRONG, " Cannot mesh 2D with simplices" );
346: refineRatio = PetscMax ( ( PetscInt ) ( PetscPowReal( numProcs, 1.0/dim) + 0.1) - 1, 1);
347: for ( d = 0; d < dim; ++d) {
348: if ( ctx->cells[d] < refineRatio) ctx->cells[d] = refineRatio;
349: if ( ctx->periodicity[d]= = DM_BOUNDARY_PERIODIC && ctx->cells[d]*refineRatio <= 2) refineRatio = 2;
350: }
351: for ( d = 0; d < dim; ++d) {
352: ctx->cells[d] *= refineRatio;
353: totCells *= ctx->cells[d];
354: }
355: if ( totCells % numProcs) SETERRQ2 ( comm,PETSC_ERR_ARG_WRONG," Total cells %D not divisible by processes %D" , totCells, numProcs);
356: DMPlexCreateBoxMesh ( comm, dim, PETSC_FALSE , ctx->cells, ctx->domain_lo, ctx->domain_hi, ctx->periodicity, PETSC_TRUE , dm);
357: } else {
358: if ( ctx->periodicity[0]= = DM_BOUNDARY_PERIODIC || ctx->periodicity[1]= = DM_BOUNDARY_PERIODIC ) SETERRQ ( comm, PETSC_ERR_ARG_WRONG, " Cannot do periodic in x or y in a cylinder" );
359: /* we stole dm_refine so clear it */
360: PetscOptionsClearValue ( NULL," -dm_refine" );
361: DMPlexCreateHexCylinderMesh ( comm, ctx->refine, ctx->periodicity[2], dm);
362: }
363: }
364: {
365: DM distributedMesh = NULL;
366: /* Distribute mesh over processes */
367: DMPlexDistribute ( *dm, 0, NULL, &distributedMesh);
368: if ( distributedMesh) {
369: DMDestroy ( dm);
370: *dm = distributedMesh;
371: }
372: }
373: {
374: PetscBool hasLabel;
375: DMHasLabel ( *dm, " marker" , &hasLabel);
376: if ( !hasLabel) {CreateBCLabel( *dm, " marker" );}
377: }
378: {
379: char convType[256];
380: PetscBool flg;
381: PetscOptionsBegin ( comm, " " , " Mesh conversion options" , " DMPLEX " );
382: PetscOptionsFList ( " -dm_plex_convert_type" ," Convert DMPlex to another format" ," ex48" ,DMList,DMPLEX ,convType,256,&flg);
383: PetscOptionsEnd ( );
384: if ( flg) {
385: DM dmConv;
386: DMConvert ( *dm,convType,&dmConv);
387: if ( dmConv) {
388: DMDestroy ( dm);
389: *dm = dmConv;
390: }
391: }
392: }
393: PetscObjectSetName ( ( PetscObject ) *dm, " Mesh" );
394: DMSetFromOptions ( *dm);
395: DMLocalizeCoordinates ( *dm); /* needed for periodic */
396: return ( 0);
397: }
399: static PetscErrorCode log_n_0( PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
400: {
401: AppCtx *lctx = ( AppCtx*)ctx;
402: assert( ctx);
403: u[0] = ( lctx->domain_hi-lctx->domain_lo)+x[0];
404: return 0;
405: }
407: static PetscErrorCode Omega_0( PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
408: {
409: u[0] = 0.0;
410: return 0;
411: }
413: static PetscErrorCode psi_0( PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
414: {
415: AppCtx *lctx = ( AppCtx*)ctx;
416: assert( ctx);
417: /* This sets up a symmetrix By flux aroound the mid point in x, which represents a current density flux along z. The stability
418: is analytically known and reported in ProcessOptions. */
419: if ( lctx->ke!= 0.0) {
420: u[0] = ( PetscCosReal( lctx->ke*( x[0]-lctx->a))-PetscCosReal( lctx->ke*lctx->a))/( 1.0-PetscCosReal( lctx->ke*lctx->a));
421: } else {
422: u[0] = 1.0-PetscPowScalar( ( x[0]-lctx->a)/lctx->a,2);
423: }
424: return 0;
425: }
427: static PetscErrorCode initialSolution_n( PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
428: {
429: u[0] = 0.0;
430: return 0;
431: }
433: static PetscErrorCode initialSolution_Omega( PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
434: {
435: u[0] = 0.0;
436: return 0;
437: }
439: static PetscErrorCode initialSolution_psi( PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *a_ctx)
440: {
441: AppCtx *ctx = ( AppCtx*)a_ctx;
442: PetscScalar r = ctx->eps*( PetscScalar ) ( rand( )) / ( PetscScalar ) ( RAND_MAX);
443: assert( ctx);
444: if ( x[0] = = ctx->domain_lo[0] || x[0] = = ctx->domain_hi[0]) r = 0;
445: u[0] = r;
446: /* PetscPrintf ( PETSC_COMM_WORLD , " rand psi %lf\n" ,u[0]); */
447: return 0;
448: }
450: static PetscErrorCode initialSolution_phi( PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
451: {
452: u[0] = 0.0;
453: return 0;
454: }
456: static PetscErrorCode initialSolution_jz( PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
457: {
458: u[0] = 0.0;
459: return 0;
460: }
462: static PetscErrorCode SetupProblem( DM dm, AppCtx *ctx)
463: {
464: PetscDS prob;
465: const PetscInt id = 1;
466: PetscErrorCode ierr, f;
469: DMGetDS ( dm, &prob);
470: PetscDSSetResidual ( prob, 0, f0_n, f1_n);
471: PetscDSSetResidual ( prob, 1, f0_Omega, f1_Omega);
472: PetscDSSetResidual ( prob, 2, f0_psi, f1_psi);
473: PetscDSSetResidual ( prob, 3, f0_phi, f1_phi);
474: PetscDSSetResidual ( prob, 4, f0_jz, f1_jz);
475: ctx->initialFuncs[0] = initialSolution_n;
476: ctx->initialFuncs[1] = initialSolution_Omega;
477: ctx->initialFuncs[2] = initialSolution_psi;
478: ctx->initialFuncs[3] = initialSolution_phi;
479: ctx->initialFuncs[4] = initialSolution_jz;
480: for ( f = 0; f < 5; ++f) {
481: PetscDSSetImplicit ( prob, f, ctx->implicit);
482: PetscDSAddBoundary ( prob, DM_BC_ESSENTIAL , " wall" , " marker" , f, 0, NULL, ( void ( *)( void)) ctx->initialFuncs[f], 1, &id, ctx);
483: }
484: PetscDSSetContext( prob, 0, ctx);
485: return ( 0);
486: }
488: static PetscErrorCode SetupEquilibriumFields( DM dm, DM dmAux, AppCtx *ctx)
489: {
490: PetscErrorCode ( *eqFuncs[3])( PetscInt , PetscReal , const PetscReal [], PetscInt , PetscScalar [], void *) = {log_n_0, Omega_0, psi_0};
491: Vec eq;
493: AppCtx *ctxarr[3];
495: ctxarr[0] = ctxarr[1] = ctxarr[2] = ctx; /* each variable could have a different context */
497: DMCreateLocalVector ( dmAux, &eq);
498: DMProjectFunctionLocal( dmAux, 0.0, eqFuncs, ( void **)ctxarr, INSERT_ALL_VALUES , eq);
499: PetscObjectCompose ( ( PetscObject ) dm, " A" , ( PetscObject ) eq);
500: if ( ctx->plotRef) { /* plot reference functions */
501: PetscViewer viewer = NULL;
502: PetscBool isHDF5,isVTK;
503: char buf[256];
504: Vec global;
505: DMCreateGlobalVector ( dmAux,&global);
506: VecSet ( global,.0); /* BCs! */
507: DMLocalToGlobalBegin ( dmAux,eq,INSERT_VALUES ,global);
508: DMLocalToGlobalEnd ( dmAux,eq,INSERT_VALUES ,global);
509: PetscViewerCreate ( PetscObjectComm ( ( PetscObject )dmAux),&viewer);
510: #ifdef PETSC_HAVE_HDF5
511: PetscViewerSetType ( viewer,PETSCVIEWERHDF5 );
512: #else
513: PetscViewerSetType ( viewer,PETSCVIEWERVTK );
514: #endif
515: PetscViewerSetFromOptions ( viewer);
516: PetscObjectTypeCompare ( ( PetscObject )viewer,PETSCVIEWERHDF5 ,&isHDF5);
517: PetscObjectTypeCompare ( ( PetscObject )viewer,PETSCVIEWERVTK ,&isVTK);
518: if ( isHDF5) {
519: PetscSNPrintf ( buf, 256, " uEquilibrium-%dD.h5" , ctx->dim);
520: } else if ( isVTK) {
521: PetscSNPrintf ( buf, 256, " uEquilibrium-%dD.vtu" , ctx->dim);
522: PetscViewerPushFormat ( viewer,PETSC_VIEWER_VTK_VTU );
523: }
524: PetscViewerFileSetMode ( viewer,FILE_MODE_WRITE );
525: PetscViewerFileSetName ( viewer,buf);
526: if ( isHDF5) {DMView ( dmAux,viewer);}
527: /* view equilibrium fields, this will overwrite fine grids with coarse grids! */
528: PetscObjectSetName ( ( PetscObject ) global, " u0" );
529: VecView ( global,viewer);
530: PetscViewerDestroy ( &viewer);
531: VecDestroy ( &global);
532: }
533: VecDestroy ( &eq);
534: return ( 0);
535: }
537: static PetscErrorCode SetupAuxDM( DM dm, PetscInt NfAux, PetscFE feAux[], AppCtx *user)
538: {
539: DM dmAux, coordDM;
540: PetscInt f;
544: /* MUST call DMGetCoordinateDM ( ) in order to get p4est setup if present */
545: DMGetCoordinateDM ( dm, &coordDM);
546: if ( !feAux) return ( 0);
547: DMClone ( dm, &dmAux);
548: PetscObjectCompose ( ( PetscObject ) dm, " dmAux" , ( PetscObject ) dmAux);
549: DMSetCoordinateDM ( dmAux, coordDM);
550: for ( f = 0; f < NfAux; ++f) {DMSetField ( dmAux, f, NULL, ( PetscObject ) feAux[f]);}
551: DMCreateDS ( dmAux);
552: SetupEquilibriumFields( dm, dmAux, user);
553: DMDestroy ( &dmAux);
554: return ( 0);
555: }
557: static PetscErrorCode SetupDiscretization( DM dm, AppCtx *ctx)
558: {
559: DM cdm = dm;
560: const PetscInt dim = ctx->dim;
561: PetscQuadrature q;
562: PetscFE fe[5], feAux[3];
563: PetscInt Nf = 5, NfAux = 3, f;
564: PetscBool cell_simplex = ctx->cell_simplex;
565: MPI_Comm comm;
566: PetscErrorCode ierr;
569: /* Create finite element */
570: PetscObjectGetComm ( ( PetscObject ) dm, &comm);
571: PetscFECreateDefault ( comm, dim, 1, cell_simplex, NULL, -1, &fe[0]);
572: PetscObjectSetName ( ( PetscObject ) fe[0], " density" );
573: PetscFECreateDefault ( comm, dim, 1, cell_simplex, NULL, -1, &fe[1]);
574: PetscObjectSetName ( ( PetscObject ) fe[1], " vorticity" );
575: PetscFECreateDefault ( comm, dim, 1, cell_simplex, NULL, -1, &fe[2]);
576: PetscObjectSetName ( ( PetscObject ) fe[2], " flux" );
577: PetscFECreateDefault ( comm, dim, 1, cell_simplex, NULL, -1, &fe[3]);
578: PetscObjectSetName ( ( PetscObject ) fe[3], " potential" );
579: PetscFECreateDefault ( comm, dim, 1, cell_simplex, NULL, -1, &fe[4]);
580: PetscObjectSetName ( ( PetscObject ) fe[4], " current" );
582: PetscFECreateDefault ( comm, dim, 1, cell_simplex, NULL, -1, &feAux[0]);
583: PetscFEGetQuadrature ( fe[0], &q);
584: PetscFESetQuadrature ( feAux[0], q);
585: PetscObjectSetName ( ( PetscObject ) feAux[0], " n_0" );
586: PetscFECreateDefault ( comm, dim, 1, cell_simplex, NULL, -1, &feAux[1]);
587: PetscFEGetQuadrature ( fe[1], &q);
588: PetscFESetQuadrature ( feAux[1], q);
589: PetscObjectSetName ( ( PetscObject ) feAux[1], " vorticity_0" );
590: PetscFECreateDefault ( comm, dim, 1, cell_simplex, NULL, -1, &feAux[2]);
591: PetscFEGetQuadrature ( fe[2], &q);
592: PetscFESetQuadrature ( feAux[2], q);
593: PetscObjectSetName ( ( PetscObject ) feAux[2], " flux_0" );
594: /* Set discretization and boundary conditions for each mesh */
595: for ( f = 0; f < Nf; ++f) {DMSetField ( dm, f, NULL, ( PetscObject ) fe[f]);}
596: DMCreateDS ( dm);
597: SetupProblem( dm, ctx);
598: while ( cdm) {
599: DMCopyDisc ( dm, cdm);
600: SetupAuxDM( dm, NfAux, feAux, ctx);
601: {
602: PetscBool hasLabel;
604: DMHasLabel ( cdm, " marker" , &hasLabel);
605: if ( !hasLabel) {CreateBCLabel( cdm, " marker" );}
606: }
607: DMGetCoarseDM ( cdm, &cdm);
608: }
609: for ( f = 0; f < Nf; ++f) {PetscFEDestroy ( &fe[f]);}
610: for ( f = 0; f < NfAux; ++f) {PetscFEDestroy ( &feAux[f]);}
611: return ( 0);
612: }
614: int main( int argc, char **argv)
615: {
616: DM dm;
617: TS ts;
618: Vec u, r;
619: AppCtx ctx;
620: PetscReal t = 0.0;
621: PetscReal L2error = 0.0;
623: AppCtx *ctxarr[5];
625: ctxarr[0] = ctxarr[1] = ctxarr[2] = ctxarr[3] = ctxarr[4] = &ctx; /* each variable could have a different context */
626: s_ctx = &ctx;
627: PetscInitialize ( &argc, &argv, NULL,help);if ( ierr) return ierr;
628: ProcessOptions( PETSC_COMM_WORLD , &ctx);
629: /* create mesh and problem */
630: CreateMesh( PETSC_COMM_WORLD , &ctx, &dm);
631: DMSetApplicationContext ( dm, &ctx);
632: PetscMalloc1 ( 5, &ctx.initialFuncs);
633: SetupDiscretization( dm, &ctx);
634: DMCreateGlobalVector ( dm, &u);
635: PetscObjectSetName ( ( PetscObject ) u, " u" );
636: VecDuplicate ( u, &r);
637: PetscObjectSetName ( ( PetscObject ) r, " r" );
638: /* create TS */
639: TSCreate ( PETSC_COMM_WORLD , &ts);
640: TSSetDM ( ts, dm);
641: TSSetApplicationContext ( ts, &ctx);
642: DMTSSetBoundaryLocal ( dm, DMPlexTSComputeBoundary , &ctx);
643: if ( ctx.implicit) {
644: DMTSSetIFunctionLocal ( dm, DMPlexTSComputeIFunctionFEM , &ctx);
645: DMTSSetIJacobianLocal ( dm, DMPlexTSComputeIJacobianFEM , &ctx);
646: } else {
647: DMTSSetRHSFunctionLocal ( dm, DMPlexTSComputeRHSFunctionFVM , &ctx);
648: }
649: TSSetExactFinalTime ( ts, TS_EXACTFINALTIME_STEPOVER );
650: TSSetFromOptions ( ts);
651: TSSetPostStep ( ts, PostStep);
652: /* make solution & solve */
653: DMProjectFunction ( dm, t, ctx.initialFuncs, ( void **)ctxarr, INSERT_ALL_VALUES , u);
654: TSSetSolution ( ts,u);
655: DMViewFromOptions( dm, NULL, " -dm_view" );
656: PostStep( ts); /* print the initial state */
657: TSSolve ( ts, u);
658: TSGetTime ( ts, &t);
659: DMComputeL2Diff ( dm, t, ctx.initialFuncs, ( void **)ctxarr, u, &L2error);
660: if ( L2error < 1.0e-11 ) {PetscPrintf ( PETSC_COMM_WORLD , " L_2 Error: < 1.0e-11 \n" );}
661: else {PetscPrintf ( PETSC_COMM_WORLD , " L_2 Error: %g\n" , L2error);}
662: #if 0
663: {
664: PetscReal res = 0.0;
665: /* Check discretization error */
666: VecViewFromOptions( u, NULL, " -initial_guess_view" );
667: DMComputeL2Diff ( dm, 0.0, ctx.exactFuncs, NULL, u, &error);
668: if ( error < 1.0e-11 ) {PetscPrintf ( PETSC_COMM_WORLD , " L_2 Error: < 1.0e-11 \n" );}
669: else {PetscPrintf ( PETSC_COMM_WORLD , " L_2 Error: %g\n" , error);}
670: /* Check residual */
671: SNESComputeFunction ( snes, u, r);
672: VecChop ( r, 1.0e-10);
673: VecViewFromOptions( r, NULL, " -initial_residual_view" );
674: VecNorm ( r, NORM_2 , &res);
675: PetscPrintf ( PETSC_COMM_WORLD , " L_2 Residual: %g\n" , res);
676: /* Check Jacobian */
677: {
678: Mat A;
679: Vec b;
681: SNESGetJacobian ( snes, &A, NULL, NULL, NULL);
682: SNESComputeJacobian ( snes, u, A, A);
683: VecDuplicate ( u, &b);
684: VecSet ( r, 0.0);
685: SNESComputeFunction ( snes, r, b);
686: MatMult ( A, u, r);
687: VecAXPY ( r, 1.0, b);
688: VecDestroy ( &b);
689: PetscPrintf ( PETSC_COMM_WORLD , " Au - b = Au + F( 0)\n" );
690: VecChop ( r, 1.0e-10);
691: VecViewFromOptions( r, NULL, " -linear_residual_view" );
692: VecNorm ( r, NORM_2 , &res);
693: PetscPrintf ( PETSC_COMM_WORLD , " Linear L_2 Residual: %g\n" , res);
694: }
695: }
696: #endif
697: VecDestroy ( &u);
698: VecDestroy ( &r);
699: TSDestroy ( &ts);
700: DMDestroy ( &dm);
701: PetscFree ( ctx.initialFuncs);
702: PetscFinalize ( );
703: return ierr;
704: }
706: /*TEST
708: test:
709: suffix: 0
710: args: -debug 1 -dim 2 -dm_refine 1 -x_periodicity PERIODIC -ts_max_steps 1 -ts_max_time 10. -ts_dt 1.0
711: test:
712: suffix: 1
713: args: -debug 1 -dim 3 -dm_refine 1 -z_periodicity PERIODIC -ts_max_steps 1 -ts_max_time 10. -ts_dt 1.0 -domain_lo -2,-1,-1 -domain_hi 2,1,1
715: TEST*/