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_INT_MAX; 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 PetscFunctionReturn(PETSC_SUCCESS); 100 } 101 102 // | 0 1 | matrix to apply bracket 103 // |-1 0 | 104 static PetscReal s_K[2][2] = { 105 {0, 1}, 106 {-1, 0} 107 }; 108 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 } 116 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 } 158 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, j; 163 const PetscScalar *pOmegaDer = &u_x[uOff_x[OMEGA]]; 164 for (i = 0; i < dim; ++i) 165 for (j = 0; j < dim; ++j) g1[j] += pOmegaDer[i] * s_K[i][j]; 166 } 167 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, j; 172 const PetscScalar *pPsiDer = &u_x[uOff_x[PSI]]; 173 for (i = 0; i < dim; ++i) 174 for (j = 0; j < dim; ++j) g1[j] += pPsiDer[i] * s_K[i][j]; 175 } 176 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 } 183 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 PetscInt d; 188 for (d = 0; d < dim; ++d) g3[d * dim + d] = 1; 189 } 190 191 /* residual point methods */ 192 static PetscScalar poissonBracket(PetscInt dim, const PetscScalar df[], const PetscScalar dg[]) 193 { 194 PetscScalar ret = df[0] * dg[1] - df[1] * dg[0]; 195 if (dim == 3) { 196 ret += df[1] * dg[2] - df[2] * dg[1]; 197 ret += df[2] * dg[0] - df[0] * dg[2]; 198 } 199 return ret; 200 } 201 // 202 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[]) 203 { 204 const PetscScalar *omegaDer = &u_x[uOff_x[OMEGA]]; 205 const PetscScalar *psiDer = &u_x[uOff_x[PSI]]; 206 const PetscScalar *phiDer = &u_x[uOff_x[PHI]]; 207 const PetscScalar *jzDer = &u_x[uOff_x[JZ]]; 208 209 f0[0] += poissonBracket(dim, omegaDer, phiDer) - poissonBracket(dim, jzDer, psiDer); 210 211 if (u_t) f0[0] += u_t[OMEGA]; 212 } 213 214 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[]) 215 { 216 const PetscScalar *omegaDer = &u_x[uOff_x[OMEGA]]; 217 PetscReal mu = PetscRealPart(constants[MU_CONST]); 218 219 for (PetscInt d = 0; d < dim; ++d) f1[d] += mu * omegaDer[d]; 220 } 221 222 // d/dt + {psi,phi} - eta j_z 223 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[]) 224 { 225 const PetscScalar *psiDer = &u_x[uOff_x[PSI]]; 226 const PetscScalar *phiDer = &u_x[uOff_x[PHI]]; 227 PetscReal eta = PetscRealPart(constants[ETA_CONST]); 228 229 f0[0] = -eta * u[uOff[JZ]]; 230 f0[0] += poissonBracket(dim, psiDer, phiDer); 231 232 if (u_t) f0[0] += u_t[PSI]; 233 // printf("psiDer = %20.15e %20.15e psi = %20.15e 234 } 235 236 // d/dt - eta j_z 237 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[]) 238 { 239 PetscReal eta = PetscRealPart(constants[ETA_CONST]); 240 241 f0[0] = -eta * u[uOff[JZ]]; 242 243 if (u_t) f0[0] += u_t[PSI]; 244 // printf("psiDer = %20.15e %20.15e psi = %20.15e 245 } 246 247 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[]) 248 { 249 f0[0] += u[uOff[OMEGA]]; 250 } 251 252 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[]) 253 { 254 const PetscScalar *phiDer = &u_x[uOff_x[PHI]]; 255 256 for (PetscInt d = 0; d < dim; ++d) f1[d] = phiDer[d]; 257 } 258 259 /* - eta M */ 260 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[]) 261 { 262 PetscReal eta = PetscRealPart(constants[ETA_CONST]); 263 264 g0[0] = -eta; 265 } 266 267 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[]) 268 { 269 f0[0] = u[uOff[JZ]]; 270 } 271 272 /* -del^2 psi = (grad v, grad psi) */ 273 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[]) 274 { 275 const PetscScalar *psiDer = &u_x[uOff_x[PSI]]; 276 277 for (PetscInt d = 0; d < dim; ++d) f1[d] = psiDer[d]; 278 } 279 280 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) 281 { 282 const PetscScalar *psiDer = &u_x[uOff_x[PSI]]; 283 PetscScalar b2 = 0; 284 for (int i = 0; i < dim; ++i) b2 += psiDer[i] * psiDer[i]; 285 f0[0] = b2; 286 } 287 288 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) 289 { 290 const PetscScalar *phiDer = &u_x[uOff_x[PHI]]; 291 PetscScalar v2 = 0; 292 for (int i = 0; i < dim; ++i) v2 += phiDer[i] * phiDer[i]; 293 f0[0] = v2; 294 } 295 296 static PetscErrorCode Monitor(TS ts, PetscInt stepi, PetscReal time, Vec X, void *actx) 297 { 298 AppCtx *ctx = (AppCtx *)actx; /* user-defined application context */ 299 SNES snes; 300 SNESConvergedReason reason; 301 TSConvergedReason tsreason; 302 303 PetscFunctionBegin; 304 // PetscCall(TSGetApplicationContext(ts, &ctx)); 305 if (ctx->debug < 1) PetscFunctionReturn(PETSC_SUCCESS); 306 PetscCall(TSGetSNES(ts, &snes)); 307 PetscCall(SNESGetConvergedReason(snes, &reason)); 308 if (reason < 0) { 309 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "\t\t ***************** Monitor: SNES diverged with reason %d.\n", (int)reason)); 310 PetscFunctionReturn(PETSC_SUCCESS); 311 } 312 if (stepi > ctx->plotStep && ctx->plotting) { 313 ctx->plotting = PETSC_FALSE; /* was doing diagnostics, now done */ 314 ctx->plotIdx++; 315 } 316 PetscCall(TSGetTime(ts, &time)); 317 PetscCall(TSGetConvergedReason(ts, &tsreason)); 318 if (((time - ctx->plotStartTime) / ctx->plotDt >= (PetscReal)ctx->plotIdx && time >= ctx->plotStartTime) || (tsreason == TS_CONVERGED_TIME || tsreason == TS_CONVERGED_ITS) || ctx->plotIdx == 0) { 319 DM dm, plex; 320 Vec X; 321 PetscReal val; 322 PetscScalar tt[12]; // FE integral seems to need a large array 323 PetscDS prob; 324 if (!ctx->plotting) { /* first step of possible backtracks */ 325 ctx->plotting = PETSC_TRUE; 326 } else { 327 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\t\t ?????? ------\n")); 328 ctx->plotting = PETSC_TRUE; 329 } 330 ctx->plotStep = stepi; 331 PetscCall(TSGetSolution(ts, &X)); 332 PetscCall(VecGetDM(X, &dm)); 333 PetscCall(DMGetOutputSequenceNumber(dm, NULL, &val)); 334 PetscCall(DMSetOutputSequenceNumber(dm, ctx->plotIdx, val)); 335 PetscCall(VecViewFromOptions(X, NULL, "-vec_view_mhd")); 336 if (ctx->debug > 2) { 337 Vec R; 338 PetscCall(SNESGetFunction(snes, &R, NULL, NULL)); 339 PetscCall(VecViewFromOptions(R, NULL, "-vec_view_res")); 340 } 341 // compute energy 342 PetscCall(DMGetDS(dm, &prob)); 343 PetscCall(DMConvert(dm, DMPLEX, &plex)); 344 PetscCall(PetscDSSetObjective(prob, 0, &f0_mhd_v_energy2)); 345 PetscCall(DMPlexComputeIntegralFEM(plex, X, &tt[0], ctx)); 346 val = PetscRealPart(tt[0]); 347 PetscCall(PetscDSSetObjective(prob, 0, &f0_mhd_B_energy2)); 348 PetscCall(DMPlexComputeIntegralFEM(plex, X, &tt[0], ctx)); 349 val = PetscSqrtReal(val) * 0.5 + PetscSqrtReal(PetscRealPart(tt[0])) * 0.5; 350 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)); 351 /* clean up */ 352 PetscCall(DMDestroy(&plex)); 353 } 354 PetscFunctionReturn(PETSC_SUCCESS); 355 } 356 357 static PetscErrorCode CreateBCLabel(DM dm, const char name[]) 358 { 359 DMLabel label; 360 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 377 PetscFunctionBeginUser; 378 PetscCall(PetscStrlen(filename, &len)); 379 if (len) { 380 PetscCall(DMPlexCreateFromFile(comm, filename, "", PETSC_TRUE, dm)); 381 } else { 382 PetscCall(DMCreate(comm, dm)); 383 PetscCall(DMSetType(*dm, DMPLEX)); 384 } 385 PetscCallMPI(MPI_Comm_size(comm, &size)); 386 while (nface * nface < size) nface *= 2; // 2D 387 if (nface < 2) nface = 2; 388 PetscCall(PetscSNPrintf(buff, sizeof(buff), "-dm_plex_box_faces %d,%d -petscpartitioner_type simple", (int)nface, (int)nface)); 389 PetscCall(PetscOptionsInsertString(NULL, buff)); 390 PetscCall(PetscOptionsInsertString(NULL, "-dm_plex_box_lower -2,-2 -dm_plex_box_upper 2,2")); 391 PetscCall(DMSetFromOptions(*dm)); 392 PetscCall(DMPlexDistributeSetDefault(*dm, PETSC_FALSE)); 393 PetscCall(DMGetDimension(*dm, &ctx->dim)); 394 { 395 char convType[256]; 396 PetscBool flg; 397 PetscOptionsBegin(comm, "", "Mesh conversion options", "DMPLEX"); 398 PetscCall(PetscOptionsFList("-dm_plex_convert_type", "Convert DMPlex to another format", "mhd", DMList, DMPLEX, convType, 256, &flg)); 399 PetscOptionsEnd(); 400 if (flg) { 401 DM dmConv; 402 PetscCall(DMConvert(*dm, convType, &dmConv)); 403 if (dmConv) { 404 PetscCall(DMDestroy(dm)); 405 *dm = dmConv; 406 } 407 } 408 } 409 PetscCall(DMLocalizeCoordinates(*dm)); /* needed for periodic */ 410 { 411 PetscBool hasLabel; 412 PetscCall(DMHasLabel(*dm, "marker", &hasLabel)); 413 if (!hasLabel) PetscCall(CreateBCLabel(*dm, "marker")); 414 } 415 PetscCall(PetscObjectSetName((PetscObject)*dm, "Mesh")); 416 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view_mhd")); 417 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view_res")); 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 persistent? */ 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 unlimited -ts_max_steps 1 -ts_max_time -ts_monitor -ts_type arkimex 637 filter: grep -v DM_ 638 639 TEST*/ 640