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"; 3 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*/ 9 10 #include <assert.h> 11 #include <petscdmplex.h> 12 #include <petscds.h> 13 #include <petscts.h> 14 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; 39 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, void *ctx); 50 PetscReal mu, eta; 51 PetscReal perturb; 52 TestType testType; 53 ModelType modelType; 54 PetscInt Nf; 55 } AppCtx; 56 57 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 58 { 59 PetscInt ii; 60 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_MAX_INT; 73 options->plotting = PETSC_FALSE; 74 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; 88 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 100 PetscFunctionReturn(PETSC_SUCCESS); 101 } 102 103 // | 0 1 | matrix to apply bracket 104 // |-1 0 | 105 static PetscReal s_K[2][2] = { 106 {0, 1}, 107 {-1, 0} 108 }; 109 110 /* 111 dt - "g0" are mass terms 112 */ 113 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[]) 114 { 115 g0[0] = u_tShift; 116 } 117 118 /* 119 Identity, Mass 120 */ 121 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[]) 122 { 123 g0[0] = 1; 124 } 125 /* 'right' Poisson bracket -<.,phi>, linearized variable is left (column), data 126 * variable right */ 127 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[]) 128 { 129 PetscInt i, j; 130 const PetscScalar *pphiDer = &u_x[uOff_x[PHI]]; // get derivative of the 'right' ("dg") and apply to 131 // live var "df" 132 for (i = 0; i < dim; ++i) 133 for (j = 0; j < dim; ++j) 134 // indexing with inner, j, index generates the left live variable [dy,-] 135 // by convension, put j index on right, with i destination: [ d/dy, 136 // -d/dx]' 137 g1[i] += s_K[i][j] * pphiDer[j]; 138 } 139 /* 'left' bracket -{jz,.}, "n" for negative, linearized variable right (column), 140 * data variable left */ 141 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[]) 142 { 143 PetscInt i, j; 144 const PetscScalar *jzDer = &u_x[uOff_x[JZ]]; // get derivative of the 'left' ("df") and apply to live 145 // var "dg" 146 for (i = 0; i < dim; ++i) 147 for (j = 0; j < dim; ++j) 148 // live right: Der[i] * K: Der[j] --> j: [d/dy, -d/dx]' 149 g1[j] += -jzDer[i] * s_K[i][j]; 150 } 151 /* 'left' Poisson bracket -< . , psi> */ 152 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[]) 153 { 154 PetscInt i, j; 155 const PetscScalar *psiDer = &u_x[uOff_x[PSI]]; 156 for (i = 0; i < dim; ++i) 157 for (j = 0; j < dim; ++j) g1[i] += -s_K[i][j] * psiDer[j]; 158 } 159 160 /* < Omega , . > */ 161 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[]) 162 { 163 PetscInt i, j; 164 const PetscScalar *pOmegaDer = &u_x[uOff_x[OMEGA]]; 165 for (i = 0; i < dim; ++i) 166 for (j = 0; j < dim; ++j) g1[j] += pOmegaDer[i] * s_K[i][j]; 167 } 168 169 /* < psi , . > */ 170 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[]) 171 { 172 PetscInt i, j; 173 const PetscScalar *pPsiDer = &u_x[uOff_x[PSI]]; 174 for (i = 0; i < dim; ++i) 175 for (j = 0; j < dim; ++j) g1[j] += pPsiDer[i] * s_K[i][j]; 176 } 177 178 // -Lapacians (resistivity), negative sign goes away from IBP 179 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[]) 180 { 181 PetscReal mu = PetscRealPart(constants[MU_CONST]); 182 for (PetscInt d = 0; d < dim; ++d) g3[d * dim + d] = mu; 183 } 184 185 // Auxilary variable = -del^2 x, negative sign goes away from IBP 186 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[]) 187 { 188 PetscInt d; 189 for (d = 0; d < dim; ++d) g3[d * dim + d] = 1; 190 } 191 192 /* residual point methods */ 193 static PetscScalar poissonBracket(PetscInt dim, const PetscScalar df[], const PetscScalar dg[]) 194 { 195 PetscScalar ret = df[0] * dg[1] - df[1] * dg[0]; 196 if (dim == 3) { 197 ret += df[1] * dg[2] - df[2] * dg[1]; 198 ret += df[2] * dg[0] - df[0] * dg[2]; 199 } 200 return ret; 201 } 202 // 203 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[]) 204 { 205 const PetscScalar *omegaDer = &u_x[uOff_x[OMEGA]]; 206 const PetscScalar *psiDer = &u_x[uOff_x[PSI]]; 207 const PetscScalar *phiDer = &u_x[uOff_x[PHI]]; 208 const PetscScalar *jzDer = &u_x[uOff_x[JZ]]; 209 210 f0[0] += poissonBracket(dim, omegaDer, phiDer) - poissonBracket(dim, jzDer, psiDer); 211 212 if (u_t) f0[0] += u_t[OMEGA]; 213 } 214 215 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[]) 216 { 217 const PetscScalar *omegaDer = &u_x[uOff_x[OMEGA]]; 218 PetscReal mu = PetscRealPart(constants[MU_CONST]); 219 220 for (PetscInt d = 0; d < dim; ++d) f1[d] += mu * omegaDer[d]; 221 } 222 223 // d/dt + {psi,phi} - eta j_z 224 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[]) 225 { 226 const PetscScalar *psiDer = &u_x[uOff_x[PSI]]; 227 const PetscScalar *phiDer = &u_x[uOff_x[PHI]]; 228 PetscReal eta = PetscRealPart(constants[ETA_CONST]); 229 230 f0[0] = -eta * u[uOff[JZ]]; 231 f0[0] += poissonBracket(dim, psiDer, phiDer); 232 233 if (u_t) f0[0] += u_t[PSI]; 234 // printf("psiDer = %20.15e %20.15e psi = %20.15e 235 } 236 237 // d/dt - eta j_z 238 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[]) 239 { 240 PetscReal eta = PetscRealPart(constants[ETA_CONST]); 241 242 f0[0] = -eta * u[uOff[JZ]]; 243 244 if (u_t) f0[0] += u_t[PSI]; 245 // printf("psiDer = %20.15e %20.15e psi = %20.15e 246 } 247 248 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[]) 249 { 250 f0[0] += u[uOff[OMEGA]]; 251 } 252 253 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[]) 254 { 255 const PetscScalar *phiDer = &u_x[uOff_x[PHI]]; 256 257 for (PetscInt d = 0; d < dim; ++d) f1[d] = phiDer[d]; 258 } 259 260 /* - eta M */ 261 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[]) 262 { 263 PetscReal eta = PetscRealPart(constants[ETA_CONST]); 264 265 g0[0] = -eta; 266 } 267 268 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[]) 269 { 270 f0[0] = u[uOff[JZ]]; 271 } 272 273 /* -del^2 psi = (grad v, grad psi) */ 274 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[]) 275 { 276 const PetscScalar *psiDer = &u_x[uOff_x[PSI]]; 277 278 for (PetscInt d = 0; d < dim; ++d) f1[d] = psiDer[d]; 279 } 280 281 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) 282 { 283 const PetscScalar *psiDer = &u_x[uOff_x[PSI]]; 284 PetscScalar b2 = 0; 285 for (int i = 0; i < dim; ++i) b2 += psiDer[i] * psiDer[i]; 286 f0[0] = b2; 287 } 288 289 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) 290 { 291 const PetscScalar *phiDer = &u_x[uOff_x[PHI]]; 292 PetscScalar v2 = 0; 293 for (int i = 0; i < dim; ++i) v2 += phiDer[i] * phiDer[i]; 294 f0[0] = v2; 295 } 296 297 static PetscErrorCode Monitor(TS ts, PetscInt stepi, PetscReal time, Vec X, void *actx) 298 { 299 AppCtx *ctx = (AppCtx *)actx; /* user-defined application context */ 300 SNES snes; 301 SNESConvergedReason reason; 302 TSConvergedReason tsreason; 303 304 PetscFunctionBegin; 305 // PetscCall(TSGetApplicationContext(ts, &ctx)); 306 if (ctx->debug < 1) PetscFunctionReturn(PETSC_SUCCESS); 307 PetscCall(TSGetSNES(ts, &snes)); 308 PetscCall(SNESGetConvergedReason(snes, &reason)); 309 if (reason < 0) { 310 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "\t\t ***************** Monitor: SNES diverged with reason %d.\n", (int)reason)); 311 PetscFunctionReturn(PETSC_SUCCESS); 312 } 313 if (stepi > ctx->plotStep && ctx->plotting) { 314 ctx->plotting = PETSC_FALSE; /* was doing diagnostics, now done */ 315 ctx->plotIdx++; 316 } 317 PetscCall(TSGetTime(ts, &time)); 318 PetscCall(TSGetConvergedReason(ts, &tsreason)); 319 if (((time - ctx->plotStartTime) / ctx->plotDt >= (PetscReal)ctx->plotIdx && time >= ctx->plotStartTime) || (tsreason == TS_CONVERGED_TIME || tsreason == TS_CONVERGED_ITS) || ctx->plotIdx == 0) { 320 DM dm, plex; 321 Vec X; 322 PetscReal val; 323 PetscScalar tt[12]; // FE integral seems to need a large array 324 PetscDS prob; 325 if (!ctx->plotting) { /* first step of possible backtracks */ 326 ctx->plotting = PETSC_TRUE; 327 } else { 328 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\t\t ?????? ------\n")); 329 ctx->plotting = PETSC_TRUE; 330 } 331 ctx->plotStep = stepi; 332 PetscCall(TSGetSolution(ts, &X)); 333 PetscCall(VecGetDM(X, &dm)); 334 PetscCall(DMGetOutputSequenceNumber(dm, NULL, &val)); 335 PetscCall(DMSetOutputSequenceNumber(dm, ctx->plotIdx, val)); 336 PetscCall(VecViewFromOptions(X, NULL, "-vec_view_mhd")); 337 if (ctx->debug > 2) { 338 Vec R; 339 PetscCall(SNESGetFunction(snes, &R, NULL, NULL)); 340 PetscCall(VecViewFromOptions(R, NULL, "-vec_view_res")); 341 } 342 // compute energy 343 PetscCall(DMGetDS(dm, &prob)); 344 PetscCall(DMConvert(dm, DMPLEX, &plex)); 345 PetscCall(PetscDSSetObjective(prob, 0, &f0_mhd_v_energy2)); 346 PetscCall(DMPlexComputeIntegralFEM(plex, X, &tt[0], ctx)); 347 val = PetscRealPart(tt[0]); 348 PetscCall(PetscDSSetObjective(prob, 0, &f0_mhd_B_energy2)); 349 PetscCall(DMPlexComputeIntegralFEM(plex, X, &tt[0], ctx)); 350 val = PetscSqrtReal(val) * 0.5 + PetscSqrtReal(PetscRealPart(tt[0])) * 0.5; 351 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)); 352 /* clean up */ 353 PetscCall(DMDestroy(&plex)); 354 } 355 PetscFunctionReturn(PETSC_SUCCESS); 356 } 357 358 static PetscErrorCode CreateBCLabel(DM dm, const char name[]) 359 { 360 DMLabel label; 361 PetscFunctionBeginUser; 362 PetscCall(DMCreateLabel(dm, name)); 363 PetscCall(DMGetLabel(dm, name, &label)); 364 PetscCall(DMPlexMarkBoundaryFaces(dm, PETSC_DETERMINE, label)); 365 PetscCall(DMPlexLabelComplete(dm, label)); 366 PetscFunctionReturn(PETSC_SUCCESS); 367 } 368 // Create mesh, dim is set here 369 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *ctx, DM *dm) 370 { 371 const char *filename = ctx->filename; 372 size_t len; 373 char buff[256]; 374 PetscMPIInt size; 375 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 418 PetscFunctionReturn(PETSC_SUCCESS); 419 } 420 421 static PetscErrorCode initialSolution_Omega(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) 422 { 423 u[0] = 0.0; 424 return PETSC_SUCCESS; 425 } 426 427 static PetscErrorCode initialSolution_Psi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *a_ctx) 428 { 429 AppCtx *ctx = (AppCtx *)a_ctx; 430 PetscReal r = 0, theta, cos_theta; 431 // k = sp.jn_zeros(1, 1)[0] 432 const PetscReal k = 3.8317059702075125; 433 for (PetscInt i = 0; i < dim; i++) r += x[i] * x[i]; 434 r = PetscSqrtReal(r); 435 // r = sqrt(dot(x,x)) 436 theta = PetscAtan2Real(x[1], x[0]); 437 cos_theta = PetscCosReal(theta); 438 // f = conditional(gt(r, 1.0), outer_f, inner_f) 439 if (r < 1.0) { 440 // inner_f = 441 // (2/(Constant(k)*bessel_J(0,Constant(k))))*bessel_J(1,Constant(k)*r)*cos_theta 442 u[0] = 2.0 / (k * j0(k)) * j1(k * r) * cos_theta; 443 } else { 444 // outer_f = (1/r - r)*cos_theta 445 u[0] = (r - 1.0 / r) * cos_theta; 446 } 447 u[0] += ctx->perturb * ((double)rand() / (double)RAND_MAX - 0.5); 448 return PETSC_SUCCESS; 449 } 450 451 static PetscErrorCode initialSolution_Phi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) 452 { 453 u[0] = 0.0; 454 return PETSC_SUCCESS; 455 } 456 457 static PetscErrorCode initialSolution_Jz(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) 458 { 459 u[0] = 0.0; 460 return PETSC_SUCCESS; 461 } 462 463 static PetscErrorCode SetupProblem(PetscDS prob, DM dm, AppCtx *ctx) 464 { 465 PetscInt f; 466 467 PetscFunctionBeginUser; 468 // for both 2 & 4 field (j_z is same) 469 PetscCall(PetscDSSetJacobian(prob, JZ, JZ, g0_1, NULL, NULL, NULL)); 470 PetscCall(PetscDSSetJacobian(prob, JZ, PSI, NULL, NULL, NULL, g3_n1)); 471 PetscCall(PetscDSSetResidual(prob, JZ, f0_jz, f1_jz)); 472 473 PetscCall(PetscDSSetJacobian(prob, PSI, JZ, g0_neta, NULL, NULL, NULL)); 474 if (ctx->modelType == ONE_FILD) { 475 PetscCall(PetscDSSetJacobian(prob, PSI, PSI, g0_dt, NULL, NULL, 476 NULL)); // remove phi term 477 478 PetscCall(PetscDSSetResidual(prob, PSI, f0_psi_2f, NULL)); 479 } else { 480 PetscCall(PetscDSSetJacobian(prob, PSI, PSI, g0_dt, g1_phi_right, NULL, NULL)); 481 PetscCall(PetscDSSetJacobian(prob, PSI, PHI, NULL, g1_psi_left, NULL, NULL)); 482 PetscCall(PetscDSSetResidual(prob, PSI, f0_psi_4f, NULL)); 483 484 PetscCall(PetscDSSetJacobian(prob, PHI, PHI, NULL, NULL, NULL, g3_n1)); 485 PetscCall(PetscDSSetJacobian(prob, PHI, OMEGA, g0_1, NULL, NULL, NULL)); 486 PetscCall(PetscDSSetResidual(prob, PHI, f0_phi, f1_phi)); 487 488 PetscCall(PetscDSSetJacobian(prob, OMEGA, OMEGA, g0_dt, g1_phi_right, NULL, g3_nmu)); 489 PetscCall(PetscDSSetJacobian(prob, OMEGA, PSI, NULL, g1_njz_left, NULL, NULL)); 490 PetscCall(PetscDSSetJacobian(prob, OMEGA, PHI, NULL, g1_omega_left, NULL, NULL)); 491 PetscCall(PetscDSSetJacobian(prob, OMEGA, JZ, NULL, g1_npsi_right, NULL, NULL)); 492 PetscCall(PetscDSSetResidual(prob, OMEGA, f0_Omega, f1_Omega)); 493 } 494 /* Setup constants - is this persistant? */ 495 { 496 PetscScalar scales[NUM_CONSTS]; // +1 adding in testType for use in the f 497 // and g functions 498 /* These could be set from the command line */ 499 scales[MU_CONST] = ctx->mu; 500 scales[ETA_CONST] = ctx->eta; 501 // scales[TEST_CONST] = (PetscReal)ctx->testType; -- how to make work with complex 502 PetscCall(PetscDSSetConstants(prob, NUM_CONSTS, scales)); 503 } 504 for (f = 0; f < ctx->Nf; ++f) { 505 ctx->initialFuncs[f] = NULL; 506 PetscCall(PetscDSSetImplicit(prob, f, PETSC_TRUE)); 507 } 508 if (ctx->testType == TEST_TILT) { 509 const PetscInt id = 1; 510 DMLabel label; 511 PetscCall(DMGetLabel(dm, "marker", &label)); 512 513 ctx->initialFuncs[JZ] = initialSolution_Jz; 514 ctx->initialFuncs[PSI] = initialSolution_Psi; 515 516 PetscCall(PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "Jz for tilt test", label, 1, &id, JZ, 0, NULL, (void (*)(void))initialSolution_Jz, NULL, ctx, NULL)); 517 PetscCall(PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "Psi for tilt test", label, 1, &id, PSI, 0, NULL, (void (*)(void))initialSolution_Psi, NULL, ctx, NULL)); 518 if (ctx->modelType == TWO_FILD) { 519 ctx->initialFuncs[OMEGA] = initialSolution_Omega; 520 ctx->initialFuncs[PHI] = initialSolution_Phi; 521 PetscCall(PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "Omega for tilt test", label, 1, &id, OMEGA, 0, NULL, (void (*)(void))initialSolution_Omega, NULL, ctx, NULL)); 522 PetscCall(PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "Phi for tilt test", label, 1, &id, PHI, 0, NULL, (void (*)(void))initialSolution_Phi, NULL, ctx, NULL)); 523 } 524 } else { 525 PetscCheck(0, PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_WRONG, "Unsupported test type: %s (%d)", testTypes[PetscMin(ctx->testType, NUM_TEST_TYPES)], (int)ctx->testType); 526 } 527 PetscCall(PetscDSSetContext(prob, 0, ctx)); 528 PetscCall(PetscDSSetFromOptions(prob)); 529 PetscFunctionReturn(PETSC_SUCCESS); 530 } 531 532 static PetscErrorCode SetupDiscretization(DM dm, AppCtx *ctx) 533 { 534 DM cdm; 535 const PetscInt dim = ctx->dim; 536 PetscFE fe[NUM_COMP]; 537 PetscDS prob; 538 PetscInt Nf = ctx->Nf, f; 539 PetscBool cell_simplex = PETSC_TRUE; 540 MPI_Comm comm = PetscObjectComm((PetscObject)dm); 541 542 PetscFunctionBeginUser; 543 /* Create finite element */ 544 PetscCall(PetscFECreateDefault(comm, dim, 1, cell_simplex, NULL, -1, &fe[JZ])); 545 PetscCall(PetscObjectSetName((PetscObject)fe[JZ], "j_z")); 546 PetscCall(DMSetField(dm, JZ, NULL, (PetscObject)fe[JZ])); 547 PetscCall(PetscFECreateDefault(comm, dim, 1, cell_simplex, NULL, -1, &fe[PSI])); 548 PetscCall(PetscObjectSetName((PetscObject)fe[PSI], "psi")); 549 PetscCall(DMSetField(dm, PSI, NULL, (PetscObject)fe[PSI])); 550 if (ctx->modelType == TWO_FILD) { 551 PetscCall(PetscFECreateDefault(comm, dim, 1, cell_simplex, NULL, -1, &fe[OMEGA])); 552 PetscCall(PetscObjectSetName((PetscObject)fe[OMEGA], "Omega")); 553 PetscCall(DMSetField(dm, OMEGA, NULL, (PetscObject)fe[OMEGA])); 554 555 PetscCall(PetscFECreateDefault(comm, dim, 1, cell_simplex, NULL, -1, &fe[PHI])); 556 PetscCall(PetscObjectSetName((PetscObject)fe[PHI], "phi")); 557 PetscCall(DMSetField(dm, PHI, NULL, (PetscObject)fe[PHI])); 558 } 559 /* Set discretization and boundary conditions for each mesh */ 560 PetscCall(DMCreateDS(dm)); 561 PetscCall(DMGetDS(dm, &prob)); 562 for (f = 0; f < Nf; ++f) PetscCall(PetscDSSetDiscretization(prob, f, (PetscObject)fe[f])); 563 PetscCall(SetupProblem(prob, dm, ctx)); 564 cdm = dm; 565 while (cdm) { 566 PetscCall(DMCopyDisc(dm, cdm)); 567 if (dm != cdm) PetscCall(PetscObjectSetName((PetscObject)cdm, "Coarse")); 568 PetscCall(DMGetCoarseDM(cdm, &cdm)); 569 } 570 for (f = 0; f < Nf; ++f) PetscCall(PetscFEDestroy(&fe[f])); 571 PetscFunctionReturn(PETSC_SUCCESS); 572 } 573 574 int main(int argc, char **argv) 575 { 576 DM dm; 577 TS ts; 578 Vec u, r; 579 AppCtx ctx; 580 PetscReal t = 0.0; 581 AppCtx *ctxarr[] = {&ctx, &ctx, &ctx, &ctx}; // each variable could have a different context 582 PetscMPIInt rank; 583 584 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 585 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 586 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &ctx)); // dim is not set 587 /* create mesh and problem */ 588 PetscCall(CreateMesh(PETSC_COMM_WORLD, &ctx, &dm)); 589 PetscCall(DMView(dm, PETSC_VIEWER_STDOUT_WORLD)); 590 PetscCall(DMSetApplicationContext(dm, &ctx)); 591 PetscCall(PetscMalloc1(ctx.Nf, &ctx.initialFuncs)); 592 PetscCall(SetupDiscretization(dm, &ctx)); 593 PetscCall(DMCreateGlobalVector(dm, &u)); 594 PetscCall(PetscObjectSetName((PetscObject)u, "u")); 595 PetscCall(VecDuplicate(u, &r)); 596 PetscCall(PetscObjectSetName((PetscObject)r, "r")); 597 /* create TS */ 598 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 599 PetscCall(TSSetDM(ts, dm)); 600 PetscCall(TSSetApplicationContext(ts, &ctx)); 601 PetscCall(DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &ctx)); 602 PetscCall(DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &ctx)); 603 PetscCall(DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &ctx)); 604 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); 605 PetscCall(TSSetMaxTime(ts, 15.0)); 606 PetscCall(TSSetFromOptions(ts)); 607 PetscCall(TSMonitorSet(ts, Monitor, &ctx, NULL)); 608 /* make solution */ 609 PetscCall(DMProjectFunction(dm, t, ctx.initialFuncs, (void **)ctxarr, INSERT_ALL_VALUES, u)); 610 ctx.perturb = 0.0; 611 PetscCall(TSSetSolution(ts, u)); 612 // solve 613 PetscCall(TSSolve(ts, u)); 614 // cleanup 615 PetscCall(VecDestroy(&u)); 616 PetscCall(VecDestroy(&r)); 617 PetscCall(TSDestroy(&ts)); 618 PetscCall(DMDestroy(&dm)); 619 PetscCall(PetscFree(ctx.initialFuncs)); 620 PetscCall(PetscFinalize()); 621 return 0; 622 } 623 624 /*TEST 625 626 test: 627 suffix: 0 628 requires: triangle !complex 629 nsize: 4 630 args: -dm_plex_box_lower -2,-2 -dm_plex_box_upper 2,2 -dm_plex_simplex 1 -dm_refine_hierarchy 2 \ 631 -eta 0.0001 -ksp_converged_reason -ksp_max_it 50 -ksp_rtol 1e-3 -ksp_type fgmres -mg_coarse_ksp_rtol 1e-1 \ 632 -mg_coarse_ksp_type fgmres -mg_coarse_mg_levels_ksp_type gmres -mg_coarse_pc_type gamg -mg_levels_ksp_max_it 4 \ 633 -mg_levels_ksp_type gmres -mg_levels_pc_type jacobi -mu 0.005 -pc_mg_type full -pc_type mg \ 634 -petscpartitioner_type simple -petscspace_degree 2 -snes_converged_reason -snes_max_it 10 -snes_monitor \ 635 -snes_rtol 1.e-9 -snes_stol 1.e-9 -ts_adapt_dt_max 0.01 -ts_adapt_monitor -ts_arkimex_type 1bee \ 636 -ts_dt 0.001 -ts_max_reject 10 -ts_max_snes_failures -1 -ts_max_steps 1 -ts_max_time -ts_monitor -ts_type arkimex 637 filter: grep -v DM_ 638 639 TEST*/ 640