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, void *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, 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: }

    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: }

    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:   PetscInt d;
    188:   for (d = 0; d < dim; ++d) g3[d * dim + d] = 1;
    189: }

    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]];

    209:   f0[0] += poissonBracket(dim, omegaDer, phiDer) - poissonBracket(dim, jzDer, psiDer);

    211:   if (u_t) f0[0] += u_t[OMEGA];
    212: }

    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]);

    219:   for (PetscInt d = 0; d < dim; ++d) f1[d] += mu * omegaDer[d];
    220: }

    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]);

    229:   f0[0] = -eta * u[uOff[JZ]];
    230:   f0[0] += poissonBracket(dim, psiDer, phiDer);

    232:   if (u_t) f0[0] += u_t[PSI];
    233:   // printf("psiDer = %20.15e %20.15e psi = %20.15e
    234: }

    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]);

    241:   f0[0] = -eta * u[uOff[JZ]];

    243:   if (u_t) f0[0] += u_t[PSI];
    244:   // printf("psiDer = %20.15e %20.15e psi = %20.15e
    245: }

    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: }

    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]];

    256:   for (PetscInt d = 0; d < dim; ++d) f1[d] = phiDer[d];
    257: }

    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]);

    264:   g0[0] = -eta;
    265: }

    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: }

    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]];

    277:   for (PetscInt d = 0; d < dim; ++d) f1[d] = psiDer[d];
    278: }

    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: }

    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: }

    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;

    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: }

    357: static PetscErrorCode CreateBCLabel(DM dm, const char name[])
    358: {
    359:   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;

    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: }

    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: }

    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: }

    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: }

    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: }

    463: static PetscErrorCode SetupProblem(PetscDS prob, DM dm, AppCtx *ctx)
    464: {
    465:   PetscInt f;

    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));

    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

    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));

    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));

    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));

    513:     ctx->initialFuncs[JZ]  = initialSolution_Jz;
    514:     ctx->initialFuncs[PSI] = initialSolution_Psi;

    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: }

    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);

    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]));

    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: }

    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;

    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: }

    624: /*TEST

    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_

    639: TEST*/