Actual source code: ex48.c
1: static char help[] = "Magnetohydrodynamics (MHD) with Poisson brackets and "
2: "stream functions, solver testbed for M3D-C1. Used in https://arxiv.org/abs/2302.10242";
4: /*F
5: The strong form of a two field model for vorticity $\Omega$ and magnetic flux
6: $\psi$, using auxiliary variables potential $\phi$ and (negative) current
7: density $j_z$ \cite{Jardin04,Strauss98}.See http://arxiv.org/abs/ for more details
8: F*/
10: #include <assert.h>
11: #include <petscdmplex.h>
12: #include <petscds.h>
13: #include <petscts.h>
15: typedef enum _testidx {
16: TEST_TILT,
17: NUM_TEST_TYPES
18: } TestType;
19: const char *testTypes[NUM_TEST_TYPES + 1] = {"tilt", "unknown"};
20: typedef enum _modelidx {
21: TWO_FILD,
22: ONE_FILD,
23: NUM_MODELS
24: } ModelType;
25: const char *modelTypes[NUM_MODELS + 1] = {"two-field", "one-field", "unknown"};
26: typedef enum _fieldidx {
27: JZ,
28: PSI,
29: PHI,
30: OMEGA,
31: NUM_COMP
32: } FieldIdx; // add more
33: typedef enum _const_idx {
34: MU_CONST,
35: ETA_CONST,
36: TEST_CONST,
37: NUM_CONSTS
38: } ConstIdx;
40: typedef struct {
41: PetscInt debug; /* The debugging level */
42: PetscReal plotDt;
43: PetscReal plotStartTime;
44: PetscInt plotIdx;
45: PetscInt plotStep;
46: PetscBool plotting;
47: PetscInt dim; /* The topological mesh dimension */
48: char filename[PETSC_MAX_PATH_LEN]; /* The optional ExodusII file */
49: PetscErrorCode (**initialFuncs)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, PetscCtx ctx);
50: PetscReal mu, eta;
51: PetscReal perturb;
52: TestType testType;
53: ModelType modelType;
54: PetscInt Nf;
55: } AppCtx;
57: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
58: {
59: PetscInt ii;
61: PetscFunctionBeginUser;
62: options->debug = 1;
63: options->filename[0] = '\0';
64: options->testType = TEST_TILT;
65: options->modelType = TWO_FILD;
66: options->mu = 0.005;
67: options->eta = 0.001;
68: options->perturb = 0;
69: options->plotDt = 0.1;
70: options->plotStartTime = 0.0;
71: options->plotIdx = 0;
72: options->plotStep = PETSC_INT_MAX;
73: options->plotting = PETSC_FALSE;
75: PetscOptionsBegin(comm, "", "MHD Problem Options", "DMPLEX");
76: PetscCall(PetscOptionsInt("-debug", "The debugging level", "mhd.c", options->debug, &options->debug, NULL));
77: ii = (PetscInt)options->testType;
78: options->testType = TEST_TILT;
79: ii = options->testType;
80: PetscCall(PetscOptionsEList("-test_type", "The test type: 'tilt' Tilt instability", "mhd.c", testTypes, NUM_TEST_TYPES, testTypes[options->testType], &ii, NULL));
81: options->testType = (TestType)ii;
82: ii = (PetscInt)options->modelType;
83: options->modelType = TWO_FILD;
84: ii = options->modelType;
85: PetscCall(PetscOptionsEList("-model_type", "The model type: 'two', 'one' field", "mhd.c", modelTypes, NUM_MODELS, modelTypes[options->modelType], &ii, NULL));
86: options->modelType = (ModelType)ii;
87: options->Nf = options->modelType == TWO_FILD ? 4 : 2;
89: PetscCall(PetscOptionsReal("-mu", "Magnetic resistivity", "mhd.c", options->mu, &options->mu, NULL));
90: PetscCall(PetscOptionsReal("-eta", "Viscosity", "mhd.c", options->eta, &options->eta, NULL));
91: PetscCall(PetscOptionsReal("-plot_dt", "Plot frequency in time", "mhd.c", options->plotDt, &options->plotDt, NULL));
92: PetscCall(PetscOptionsReal("-plot_start_time", "Time to delay start of plotting", "mhd.c", options->plotStartTime, &options->plotStartTime, NULL));
93: PetscCall(PetscOptionsReal("-perturbation", "Random perturbation of initial psi scale", "mhd.c", options->perturb, &options->perturb, NULL));
94: PetscCall(PetscPrintf(comm, "Test Type = %s\n", testTypes[options->testType]));
95: PetscCall(PetscPrintf(comm, "Model Type = %s\n", modelTypes[options->modelType]));
96: PetscCall(PetscPrintf(comm, "eta = %g\n", (double)options->eta));
97: PetscCall(PetscPrintf(comm, "mu = %g\n", (double)options->mu));
98: PetscOptionsEnd();
99: PetscFunctionReturn(PETSC_SUCCESS);
100: }
102: // | 0 1 | matrix to apply bracket
103: // |-1 0 |
104: static PetscReal s_K[2][2] = {
105: {0, 1},
106: {-1, 0}
107: };
109: /*
110: dt - "g0" are mass terms
111: */
112: static void g0_dt(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
113: {
114: g0[0] = u_tShift;
115: }
117: /*
118: Identity, Mass
119: */
120: static void g0_1(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
121: {
122: g0[0] = 1;
123: }
124: /* 'right' Poisson bracket -<.,phi>, linearized variable is left (column), data
125: * variable right */
126: static void g1_phi_right(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[])
127: {
128: PetscInt i, j;
129: const PetscScalar *pphiDer = &u_x[uOff_x[PHI]]; // get derivative of the 'right' ("dg") and apply to
130: // live var "df"
131: for (i = 0; i < dim; ++i)
132: for (j = 0; j < dim; ++j)
133: // indexing with inner, j, index generates the left live variable [dy,-]
134: // by convention, put j index on right, with i destination: [ d/dy,
135: // -d/dx]'
136: g1[i] += s_K[i][j] * pphiDer[j];
137: }
138: /* 'left' bracket -{jz,.}, "n" for negative, linearized variable right (column),
139: * data variable left */
140: static void g1_njz_left(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[])
141: {
142: PetscInt i, j;
143: const PetscScalar *jzDer = &u_x[uOff_x[JZ]]; // get derivative of the 'left' ("df") and apply to live
144: // var "dg"
145: for (i = 0; i < dim; ++i)
146: for (j = 0; j < dim; ++j)
147: // live right: Der[i] * K: Der[j] --> j: [d/dy, -d/dx]'
148: g1[j] += -jzDer[i] * s_K[i][j];
149: }
150: /* 'left' Poisson bracket -< . , psi> */
151: static void g1_npsi_right(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[])
152: {
153: PetscInt i, j;
154: const PetscScalar *psiDer = &u_x[uOff_x[PSI]];
155: for (i = 0; i < dim; ++i)
156: for (j = 0; j < dim; ++j) g1[i] += -s_K[i][j] * psiDer[j];
157: }
159: /* < Omega , . > */
160: static void g1_omega_left(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[])
161: {
162: PetscInt i;
163: const PetscScalar *pOmegaDer = &u_x[uOff_x[OMEGA]];
164: for (i = 0; i < dim; ++i)
165: for (PetscInt j = 0; j < dim; ++j) g1[j] += pOmegaDer[i] * s_K[i][j];
166: }
168: /* < psi , . > */
169: static void g1_psi_left(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[])
170: {
171: PetscInt i;
172: const PetscScalar *pPsiDer = &u_x[uOff_x[PSI]];
173: for (i = 0; i < dim; ++i)
174: for (PetscInt j = 0; j < dim; ++j) g1[j] += pPsiDer[i] * s_K[i][j];
175: }
177: // -Lapacians (resistivity), negative sign goes away from IBP
178: static void g3_nmu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
179: {
180: PetscReal mu = PetscRealPart(constants[MU_CONST]);
181: for (PetscInt d = 0; d < dim; ++d) g3[d * dim + d] = mu;
182: }
184: // Auxiliary variable = -del^2 x, negative sign goes away from IBP
185: static void g3_n1(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
186: {
187: for (PetscInt d = 0; d < dim; ++d) g3[d * dim + d] = 1;
188: }
190: /* residual point methods */
191: static PetscScalar poissonBracket(PetscInt dim, const PetscScalar df[], const PetscScalar dg[])
192: {
193: PetscScalar ret = df[0] * dg[1] - df[1] * dg[0];
194: if (dim == 3) {
195: ret += df[1] * dg[2] - df[2] * dg[1];
196: ret += df[2] * dg[0] - df[0] * dg[2];
197: }
198: return ret;
199: }
200: //
201: static void f0_Omega(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
202: {
203: const PetscScalar *omegaDer = &u_x[uOff_x[OMEGA]];
204: const PetscScalar *psiDer = &u_x[uOff_x[PSI]];
205: const PetscScalar *phiDer = &u_x[uOff_x[PHI]];
206: const PetscScalar *jzDer = &u_x[uOff_x[JZ]];
208: f0[0] += poissonBracket(dim, omegaDer, phiDer) - poissonBracket(dim, jzDer, psiDer);
210: if (u_t) f0[0] += u_t[OMEGA];
211: }
213: static void f1_Omega(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
214: {
215: const PetscScalar *omegaDer = &u_x[uOff_x[OMEGA]];
216: PetscReal mu = PetscRealPart(constants[MU_CONST]);
218: for (PetscInt d = 0; d < dim; ++d) f1[d] += mu * omegaDer[d];
219: }
221: // d/dt + {psi,phi} - eta j_z
222: static void f0_psi_4f(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
223: {
224: const PetscScalar *psiDer = &u_x[uOff_x[PSI]];
225: const PetscScalar *phiDer = &u_x[uOff_x[PHI]];
226: PetscReal eta = PetscRealPart(constants[ETA_CONST]);
228: f0[0] = -eta * u[uOff[JZ]];
229: f0[0] += poissonBracket(dim, psiDer, phiDer);
231: if (u_t) f0[0] += u_t[PSI];
232: // printf("psiDer = %20.15e %20.15e psi = %20.15e
233: }
235: // d/dt - eta j_z
236: static void f0_psi_2f(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
237: {
238: PetscReal eta = PetscRealPart(constants[ETA_CONST]);
240: f0[0] = -eta * u[uOff[JZ]];
242: if (u_t) f0[0] += u_t[PSI];
243: // printf("psiDer = %20.15e %20.15e psi = %20.15e
244: }
246: static void f0_phi(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
247: {
248: f0[0] += u[uOff[OMEGA]];
249: }
251: static void f1_phi(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
252: {
253: const PetscScalar *phiDer = &u_x[uOff_x[PHI]];
255: for (PetscInt d = 0; d < dim; ++d) f1[d] = phiDer[d];
256: }
258: /* - eta M */
259: static void g0_neta(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
260: {
261: PetscReal eta = PetscRealPart(constants[ETA_CONST]);
263: g0[0] = -eta;
264: }
266: static void f0_jz(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
267: {
268: f0[0] = u[uOff[JZ]];
269: }
271: /* -del^2 psi = (grad v, grad psi) */
272: static void f1_jz(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
273: {
274: const PetscScalar *psiDer = &u_x[uOff_x[PSI]];
276: for (PetscInt d = 0; d < dim; ++d) f1[d] = psiDer[d];
277: }
279: static void f0_mhd_B_energy2(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar *f0)
280: {
281: const PetscScalar *psiDer = &u_x[uOff_x[PSI]];
282: PetscScalar b2 = 0;
283: for (int i = 0; i < dim; ++i) b2 += psiDer[i] * psiDer[i];
284: f0[0] = b2;
285: }
287: static void f0_mhd_v_energy2(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar *f0)
288: {
289: const PetscScalar *phiDer = &u_x[uOff_x[PHI]];
290: PetscScalar v2 = 0;
291: for (int i = 0; i < dim; ++i) v2 += phiDer[i] * phiDer[i];
292: f0[0] = v2;
293: }
295: static PetscErrorCode Monitor(TS ts, PetscInt stepi, PetscReal time, Vec X, void *actx)
296: {
297: AppCtx *ctx = (AppCtx *)actx; /* user-defined application context */
298: SNES snes;
299: SNESConvergedReason reason;
300: TSConvergedReason tsreason;
302: PetscFunctionBegin;
303: // PetscCall(TSGetApplicationContext(ts, &ctx));
304: if (ctx->debug < 1) PetscFunctionReturn(PETSC_SUCCESS);
305: PetscCall(TSGetSNES(ts, &snes));
306: PetscCall(SNESGetConvergedReason(snes, &reason));
307: if (reason < 0) {
308: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "\t\t ***************** Monitor: SNES diverged with reason %d.\n", (int)reason));
309: PetscFunctionReturn(PETSC_SUCCESS);
310: }
311: if (stepi > ctx->plotStep && ctx->plotting) {
312: ctx->plotting = PETSC_FALSE; /* was doing diagnostics, now done */
313: ctx->plotIdx++;
314: }
315: PetscCall(TSGetTime(ts, &time));
316: PetscCall(TSGetConvergedReason(ts, &tsreason));
317: if (((time - ctx->plotStartTime) / ctx->plotDt >= (PetscReal)ctx->plotIdx && time >= ctx->plotStartTime) || (tsreason == TS_CONVERGED_TIME || tsreason == TS_CONVERGED_ITS) || ctx->plotIdx == 0) {
318: DM dm, plex;
319: Vec X;
320: PetscReal val;
321: PetscScalar tt[12]; // FE integral seems to need a large array
322: PetscDS prob;
323: if (!ctx->plotting) { /* first step of possible backtracks */
324: ctx->plotting = PETSC_TRUE;
325: } else {
326: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\t\t ?????? ------\n"));
327: ctx->plotting = PETSC_TRUE;
328: }
329: ctx->plotStep = stepi;
330: PetscCall(TSGetSolution(ts, &X));
331: PetscCall(VecGetDM(X, &dm));
332: PetscCall(DMGetOutputSequenceNumber(dm, NULL, &val));
333: PetscCall(DMSetOutputSequenceNumber(dm, ctx->plotIdx, val));
334: PetscCall(VecViewFromOptions(X, NULL, "-vec_view_mhd"));
335: if (ctx->debug > 2) {
336: Vec R;
337: PetscCall(SNESGetFunction(snes, &R, NULL, NULL));
338: PetscCall(VecViewFromOptions(R, NULL, "-vec_view_res"));
339: }
340: // compute energy
341: PetscCall(DMGetDS(dm, &prob));
342: PetscCall(DMConvert(dm, DMPLEX, &plex));
343: PetscCall(PetscDSSetObjective(prob, 0, &f0_mhd_v_energy2));
344: PetscCall(DMPlexComputeIntegralFEM(plex, X, &tt[0], ctx));
345: val = PetscRealPart(tt[0]);
346: PetscCall(PetscDSSetObjective(prob, 0, &f0_mhd_B_energy2));
347: PetscCall(DMPlexComputeIntegralFEM(plex, X, &tt[0], ctx));
348: val = PetscSqrtReal(val) * 0.5 + PetscSqrtReal(PetscRealPart(tt[0])) * 0.5;
349: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "MHD %4d) time = %9.5g, Eergy= %20.13e (plot ID %d)\n", (int)ctx->plotIdx, (double)time, (double)val, (int)ctx->plotIdx));
350: /* clean up */
351: PetscCall(DMDestroy(&plex));
352: }
353: PetscFunctionReturn(PETSC_SUCCESS);
354: }
356: static PetscErrorCode CreateBCLabel(DM dm, const char name[])
357: {
358: DMLabel label;
360: PetscFunctionBeginUser;
361: PetscCall(DMCreateLabel(dm, name));
362: PetscCall(DMGetLabel(dm, name, &label));
363: PetscCall(DMPlexMarkBoundaryFaces(dm, PETSC_DETERMINE, label));
364: PetscCall(DMPlexLabelComplete(dm, label));
365: PetscFunctionReturn(PETSC_SUCCESS);
366: }
367: // Create mesh, dim is set here
368: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *ctx, DM *dm)
369: {
370: const char *filename = ctx->filename;
371: size_t len;
372: char buff[256];
373: PetscMPIInt size;
374: PetscInt nface = 1;
376: PetscFunctionBeginUser;
377: PetscCall(PetscStrlen(filename, &len));
378: if (len) {
379: PetscCall(DMPlexCreateFromFile(comm, filename, "", PETSC_TRUE, dm));
380: } else {
381: PetscCall(DMCreate(comm, dm));
382: PetscCall(DMSetType(*dm, DMPLEX));
383: }
384: PetscCallMPI(MPI_Comm_size(comm, &size));
385: while (nface * nface < size) nface *= 2; // 2D
386: if (nface < 2) nface = 2;
387: PetscCall(PetscSNPrintf(buff, sizeof(buff), "-dm_plex_box_faces %d,%d -petscpartitioner_type simple", (int)nface, (int)nface));
388: PetscCall(PetscOptionsInsertString(NULL, buff));
389: PetscCall(PetscOptionsInsertString(NULL, "-dm_plex_box_lower -2,-2 -dm_plex_box_upper 2,2"));
390: PetscCall(DMSetFromOptions(*dm));
391: PetscCall(DMPlexDistributeSetDefault(*dm, PETSC_FALSE));
392: PetscCall(DMGetDimension(*dm, &ctx->dim));
393: {
394: char convType[256];
395: PetscBool flg;
396: PetscOptionsBegin(comm, "", "Mesh conversion options", "DMPLEX");
397: PetscCall(PetscOptionsFList("-dm_plex_convert_type", "Convert DMPlex to another format", "mhd", DMList, DMPLEX, convType, 256, &flg));
398: PetscOptionsEnd();
399: if (flg) {
400: DM dmConv;
401: PetscCall(DMConvert(*dm, convType, &dmConv));
402: if (dmConv) {
403: PetscCall(DMDestroy(dm));
404: *dm = dmConv;
405: }
406: }
407: }
408: PetscCall(DMLocalizeCoordinates(*dm)); /* needed for periodic */
409: {
410: PetscBool hasLabel;
411: PetscCall(DMHasLabel(*dm, "marker", &hasLabel));
412: if (!hasLabel) PetscCall(CreateBCLabel(*dm, "marker"));
413: }
414: PetscCall(PetscObjectSetName((PetscObject)*dm, "Mesh"));
415: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view_mhd"));
416: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view_res"));
417: PetscFunctionReturn(PETSC_SUCCESS);
418: }
420: static PetscErrorCode initialSolution_Omega(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, PetscCtx ctx)
421: {
422: u[0] = 0.0;
423: return PETSC_SUCCESS;
424: }
426: static PetscErrorCode initialSolution_Psi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *a_ctx)
427: {
428: AppCtx *ctx = (AppCtx *)a_ctx;
429: PetscReal r = 0, theta, cos_theta;
430: // k = sp.jn_zeros(1, 1)[0]
431: const PetscReal k = 3.8317059702075125;
432: for (PetscInt i = 0; i < dim; i++) r += x[i] * x[i];
433: r = PetscSqrtReal(r);
434: // r = sqrt(dot(x,x))
435: theta = PetscAtan2Real(x[1], x[0]);
436: cos_theta = PetscCosReal(theta);
437: // f = conditional(gt(r, 1.0), outer_f, inner_f)
438: if (r < 1.0) {
439: // inner_f =
440: // (2/(Constant(k)*bessel_J(0,Constant(k))))*bessel_J(1,Constant(k)*r)*cos_theta
441: u[0] = 2.0 / (k * j0(k)) * j1(k * r) * cos_theta;
442: } else {
443: // outer_f = (1/r - r)*cos_theta
444: u[0] = (r - 1.0 / r) * cos_theta;
445: }
446: u[0] += ctx->perturb * ((double)rand() / (double)RAND_MAX - 0.5);
447: return PETSC_SUCCESS;
448: }
450: static PetscErrorCode initialSolution_Phi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, PetscCtx ctx)
451: {
452: u[0] = 0.0;
453: return PETSC_SUCCESS;
454: }
456: static PetscErrorCode initialSolution_Jz(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, PetscCtx ctx)
457: {
458: u[0] = 0.0;
459: return PETSC_SUCCESS;
460: }
462: static PetscErrorCode SetupProblem(PetscDS prob, DM dm, AppCtx *ctx)
463: {
464: PetscInt f;
466: PetscFunctionBeginUser;
467: // for both 2 & 4 field (j_z is same)
468: PetscCall(PetscDSSetJacobian(prob, JZ, JZ, g0_1, NULL, NULL, NULL));
469: PetscCall(PetscDSSetJacobian(prob, JZ, PSI, NULL, NULL, NULL, g3_n1));
470: PetscCall(PetscDSSetResidual(prob, JZ, f0_jz, f1_jz));
472: PetscCall(PetscDSSetJacobian(prob, PSI, JZ, g0_neta, NULL, NULL, NULL));
473: if (ctx->modelType == ONE_FILD) {
474: PetscCall(PetscDSSetJacobian(prob, PSI, PSI, g0_dt, NULL, NULL,
475: NULL)); // remove phi term
477: PetscCall(PetscDSSetResidual(prob, PSI, f0_psi_2f, NULL));
478: } else {
479: PetscCall(PetscDSSetJacobian(prob, PSI, PSI, g0_dt, g1_phi_right, NULL, NULL));
480: PetscCall(PetscDSSetJacobian(prob, PSI, PHI, NULL, g1_psi_left, NULL, NULL));
481: PetscCall(PetscDSSetResidual(prob, PSI, f0_psi_4f, NULL));
483: PetscCall(PetscDSSetJacobian(prob, PHI, PHI, NULL, NULL, NULL, g3_n1));
484: PetscCall(PetscDSSetJacobian(prob, PHI, OMEGA, g0_1, NULL, NULL, NULL));
485: PetscCall(PetscDSSetResidual(prob, PHI, f0_phi, f1_phi));
487: PetscCall(PetscDSSetJacobian(prob, OMEGA, OMEGA, g0_dt, g1_phi_right, NULL, g3_nmu));
488: PetscCall(PetscDSSetJacobian(prob, OMEGA, PSI, NULL, g1_njz_left, NULL, NULL));
489: PetscCall(PetscDSSetJacobian(prob, OMEGA, PHI, NULL, g1_omega_left, NULL, NULL));
490: PetscCall(PetscDSSetJacobian(prob, OMEGA, JZ, NULL, g1_npsi_right, NULL, NULL));
491: PetscCall(PetscDSSetResidual(prob, OMEGA, f0_Omega, f1_Omega));
492: }
493: /* Setup constants - is this persistent? */
494: {
495: PetscScalar scales[NUM_CONSTS]; // +1 adding in testType for use in the f
496: // and g functions
497: /* These could be set from the command line */
498: scales[MU_CONST] = ctx->mu;
499: scales[ETA_CONST] = ctx->eta;
500: // scales[TEST_CONST] = (PetscReal)ctx->testType; -- how to make work with complex
501: PetscCall(PetscDSSetConstants(prob, NUM_CONSTS, scales));
502: }
503: for (f = 0; f < ctx->Nf; ++f) {
504: ctx->initialFuncs[f] = NULL;
505: PetscCall(PetscDSSetImplicit(prob, f, PETSC_TRUE));
506: }
507: if (ctx->testType == TEST_TILT) {
508: const PetscInt id = 1;
509: DMLabel label;
510: PetscCall(DMGetLabel(dm, "marker", &label));
512: ctx->initialFuncs[JZ] = initialSolution_Jz;
513: ctx->initialFuncs[PSI] = initialSolution_Psi;
515: PetscCall(PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "Jz for tilt test", label, 1, &id, JZ, 0, NULL, (PetscVoidFn *)initialSolution_Jz, NULL, ctx, NULL));
516: PetscCall(PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "Psi for tilt test", label, 1, &id, PSI, 0, NULL, (PetscVoidFn *)initialSolution_Psi, NULL, ctx, NULL));
517: if (ctx->modelType == TWO_FILD) {
518: ctx->initialFuncs[OMEGA] = initialSolution_Omega;
519: ctx->initialFuncs[PHI] = initialSolution_Phi;
520: PetscCall(PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "Omega for tilt test", label, 1, &id, OMEGA, 0, NULL, (PetscVoidFn *)initialSolution_Omega, NULL, ctx, NULL));
521: PetscCall(PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "Phi for tilt test", label, 1, &id, PHI, 0, NULL, (PetscVoidFn *)initialSolution_Phi, NULL, ctx, NULL));
522: }
523: } else {
524: PetscCheck(0, PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_WRONG, "Unsupported test type: %s (%d)", testTypes[PetscMin(ctx->testType, NUM_TEST_TYPES)], (int)ctx->testType);
525: }
526: PetscCall(PetscDSSetContext(prob, 0, ctx));
527: PetscCall(PetscDSSetFromOptions(prob));
528: PetscFunctionReturn(PETSC_SUCCESS);
529: }
531: static PetscErrorCode SetupDiscretization(DM dm, AppCtx *ctx)
532: {
533: DM cdm;
534: const PetscInt dim = ctx->dim;
535: PetscFE fe[NUM_COMP];
536: PetscDS prob;
537: PetscInt Nf = ctx->Nf, f;
538: PetscBool cell_simplex = PETSC_TRUE;
539: MPI_Comm comm = PetscObjectComm((PetscObject)dm);
541: PetscFunctionBeginUser;
542: /* Create finite element */
543: PetscCall(PetscFECreateDefault(comm, dim, 1, cell_simplex, NULL, -1, &fe[JZ]));
544: PetscCall(PetscObjectSetName((PetscObject)fe[JZ], "j_z"));
545: PetscCall(DMSetField(dm, JZ, NULL, (PetscObject)fe[JZ]));
546: PetscCall(PetscFECreateDefault(comm, dim, 1, cell_simplex, NULL, -1, &fe[PSI]));
547: PetscCall(PetscObjectSetName((PetscObject)fe[PSI], "psi"));
548: PetscCall(DMSetField(dm, PSI, NULL, (PetscObject)fe[PSI]));
549: if (ctx->modelType == TWO_FILD) {
550: PetscCall(PetscFECreateDefault(comm, dim, 1, cell_simplex, NULL, -1, &fe[OMEGA]));
551: PetscCall(PetscObjectSetName((PetscObject)fe[OMEGA], "Omega"));
552: PetscCall(DMSetField(dm, OMEGA, NULL, (PetscObject)fe[OMEGA]));
554: PetscCall(PetscFECreateDefault(comm, dim, 1, cell_simplex, NULL, -1, &fe[PHI]));
555: PetscCall(PetscObjectSetName((PetscObject)fe[PHI], "phi"));
556: PetscCall(DMSetField(dm, PHI, NULL, (PetscObject)fe[PHI]));
557: }
558: /* Set discretization and boundary conditions for each mesh */
559: PetscCall(DMCreateDS(dm));
560: PetscCall(DMGetDS(dm, &prob));
561: for (f = 0; f < Nf; ++f) PetscCall(PetscDSSetDiscretization(prob, f, (PetscObject)fe[f]));
562: PetscCall(SetupProblem(prob, dm, ctx));
563: cdm = dm;
564: while (cdm) {
565: PetscCall(DMCopyDisc(dm, cdm));
566: if (dm != cdm) PetscCall(PetscObjectSetName((PetscObject)cdm, "Coarse"));
567: PetscCall(DMGetCoarseDM(cdm, &cdm));
568: }
569: for (f = 0; f < Nf; ++f) PetscCall(PetscFEDestroy(&fe[f]));
570: PetscFunctionReturn(PETSC_SUCCESS);
571: }
573: int main(int argc, char **argv)
574: {
575: DM dm;
576: TS ts;
577: Vec u, r;
578: AppCtx ctx;
579: PetscReal t = 0.0;
580: AppCtx *ctxarr[] = {&ctx, &ctx, &ctx, &ctx}; // each variable could have a different context
581: PetscMPIInt rank;
583: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
584: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
585: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &ctx)); // dim is not set
586: /* create mesh and problem */
587: PetscCall(CreateMesh(PETSC_COMM_WORLD, &ctx, &dm));
588: PetscCall(DMView(dm, PETSC_VIEWER_STDOUT_WORLD));
589: PetscCall(DMSetApplicationContext(dm, &ctx));
590: PetscCall(PetscMalloc1(ctx.Nf, &ctx.initialFuncs));
591: PetscCall(SetupDiscretization(dm, &ctx));
592: PetscCall(DMCreateGlobalVector(dm, &u));
593: PetscCall(PetscObjectSetName((PetscObject)u, "u"));
594: PetscCall(VecDuplicate(u, &r));
595: PetscCall(PetscObjectSetName((PetscObject)r, "r"));
596: /* create TS */
597: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
598: PetscCall(TSSetDM(ts, dm));
599: PetscCall(TSSetApplicationContext(ts, &ctx));
600: PetscCall(DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &ctx));
601: PetscCall(DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &ctx));
602: PetscCall(DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &ctx));
603: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
604: PetscCall(TSSetMaxTime(ts, 15.0));
605: PetscCall(TSSetFromOptions(ts));
606: PetscCall(TSMonitorSet(ts, Monitor, &ctx, NULL));
607: /* make solution */
608: PetscCall(DMProjectFunction(dm, t, ctx.initialFuncs, (void **)ctxarr, INSERT_ALL_VALUES, u));
609: ctx.perturb = 0.0;
610: PetscCall(TSSetSolution(ts, u));
611: // solve
612: PetscCall(TSSolve(ts, u));
613: // cleanup
614: PetscCall(VecDestroy(&u));
615: PetscCall(VecDestroy(&r));
616: PetscCall(TSDestroy(&ts));
617: PetscCall(DMDestroy(&dm));
618: PetscCall(PetscFree(ctx.initialFuncs));
619: PetscCall(PetscFinalize());
620: return 0;
621: }
623: /*TEST
625: test:
626: suffix: 0
627: requires: triangle !complex
628: nsize: 4
629: args: -dm_plex_box_lower -2,-2 -dm_plex_box_upper 2,2 -dm_plex_simplex 1 -dm_refine_hierarchy 2 \
630: -eta 0.0001 -ksp_converged_reason -ksp_max_it 50 -ksp_rtol 1e-3 -ksp_type fgmres -mg_coarse_ksp_rtol 1e-1 \
631: -mg_coarse_ksp_type fgmres -mg_coarse_mg_levels_ksp_type gmres -mg_coarse_pc_type gamg -mg_levels_ksp_max_it 4 \
632: -mg_levels_ksp_type gmres -mg_levels_pc_type jacobi -mu 0.005 -pc_mg_type full -pc_type mg \
633: -petscpartitioner_type simple -petscspace_degree 2 -snes_converged_reason -snes_max_it 10 -snes_monitor \
634: -snes_rtol 1.e-9 -snes_stol 1.e-9 -ts_adapt_dt_max 0.01 -ts_adapt_monitor -ts_arkimex_type 1bee \
635: -ts_time_step 0.001 -ts_max_step_rejections 10 -ts_max_snes_failures unlimited -ts_max_steps 1 -ts_max_time -ts_monitor -ts_type arkimex
636: filter: grep -v DM_
638: TEST*/