1 static char help[] = "Evolution of magnetic islands.\n\ 2 The aim of this model is to self-consistently study the interaction between the tearing mode and small scale drift-wave turbulence.\n\n\n"; 3 4 /*F 5 This is a three field model for the density $\tilde n$, vorticity $\tilde\Omega$, and magnetic flux $\tilde\psi$, using auxiliary variables potential $\tilde\phi$ and current $j_z$. 6 \begin{equation} 7 \begin{aligned} 8 \partial_t \tilde n &= \left\{ \tilde n, \tilde\phi \right\} + \beta \left\{ j_z, \tilde\psi \right\} + \left\{ \ln n_0, \tilde\phi \right\} + \mu \nabla^2_\perp \tilde n \\ 9 \partial_t \tilde\Omega &= \left\{ \tilde\Omega, \tilde\phi \right\} + \beta \left\{ j_z, \tilde\psi \right\} + \mu \nabla^2_\perp \tilde\Omega \\ 10 \partial_t \tilde\psi &= \left\{ \psi_0 + \tilde\psi, \tilde\phi - \tilde n \right\} - \left\{ \ln n_0, \tilde\psi \right\} + \frac{\eta}{\beta} \nabla^2_\perp \tilde\psi \\ 11 \nabla^2_\perp\tilde\phi &= \tilde\Omega \\ 12 j_z &= -\nabla^2_\perp \left(\tilde\psi + \psi_0 \right)\\ 13 \end{aligned} 14 \end{equation} 15 F*/ 16 17 #include <petscdmplex.h> 18 #include <petscts.h> 19 #include <petscds.h> 20 21 typedef struct { 22 PetscInt debug; /* The debugging level */ 23 PetscBool plotRef; /* Plot the reference fields */ 24 PetscReal lower[3], upper[3]; 25 /* Problem definition */ 26 PetscErrorCode (**initialFuncs)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx); 27 PetscReal mu, eta, beta; 28 PetscReal a, b, Jo, Jop, m, ke, kx, ky, DeltaPrime, eps; 29 /* solver */ 30 PetscBool implicit; 31 } AppCtx; 32 33 static AppCtx *s_ctx; 34 35 static PetscScalar poissonBracket(PetscInt dim, const PetscScalar df[], const PetscScalar dg[]) { 36 PetscScalar ret = df[0] * dg[1] - df[1] * dg[0]; 37 return ret; 38 } 39 40 enum field_idx { 41 DENSITY, 42 OMEGA, 43 PSI, 44 PHI, 45 JZ 46 }; 47 48 static void f0_n(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[]) { 49 const PetscScalar *pnDer = &u_x[uOff_x[DENSITY]]; 50 const PetscScalar *ppsiDer = &u_x[uOff_x[PSI]]; 51 const PetscScalar *pphiDer = &u_x[uOff_x[PHI]]; 52 const PetscScalar *jzDer = &u_x[uOff_x[JZ]]; 53 const PetscScalar *logRefDenDer = &a_x[aOff_x[DENSITY]]; 54 f0[0] += -poissonBracket(dim, pnDer, pphiDer) - s_ctx->beta * poissonBracket(dim, jzDer, ppsiDer) - poissonBracket(dim, logRefDenDer, pphiDer); 55 if (u_t) f0[0] += u_t[DENSITY]; 56 } 57 58 static void f1_n(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[]) { 59 const PetscScalar *pnDer = &u_x[uOff_x[DENSITY]]; 60 PetscInt d; 61 62 for (d = 0; d < dim - 1; ++d) f1[d] = -s_ctx->mu * pnDer[d]; 63 } 64 65 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[]) { 66 const PetscScalar *pOmegaDer = &u_x[uOff_x[OMEGA]]; 67 const PetscScalar *ppsiDer = &u_x[uOff_x[PSI]]; 68 const PetscScalar *pphiDer = &u_x[uOff_x[PHI]]; 69 const PetscScalar *jzDer = &u_x[uOff_x[JZ]]; 70 71 f0[0] += -poissonBracket(dim, pOmegaDer, pphiDer) - s_ctx->beta * poissonBracket(dim, jzDer, ppsiDer); 72 if (u_t) f0[0] += u_t[OMEGA]; 73 } 74 75 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[]) { 76 const PetscScalar *pOmegaDer = &u_x[uOff_x[OMEGA]]; 77 PetscInt d; 78 79 for (d = 0; d < dim - 1; ++d) f1[d] = -s_ctx->mu * pOmegaDer[d]; 80 } 81 82 static void f0_psi(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[]) { 83 const PetscScalar *pnDer = &u_x[uOff_x[DENSITY]]; 84 const PetscScalar *ppsiDer = &u_x[uOff_x[PSI]]; 85 const PetscScalar *pphiDer = &u_x[uOff_x[PHI]]; 86 const PetscScalar *refPsiDer = &a_x[aOff_x[PSI]]; 87 const PetscScalar *logRefDenDer = &a_x[aOff_x[DENSITY]]; 88 PetscScalar psiDer[3]; 89 PetscScalar phi_n_Der[3]; 90 PetscInt d; 91 if (dim < 2) { 92 MPI_Abort(MPI_COMM_WORLD, 1); 93 return; 94 } /* this is needed so that the clang static analyzer does not generate a warning about variables used by not set */ 95 for (d = 0; d < dim; ++d) { 96 psiDer[d] = refPsiDer[d] + ppsiDer[d]; 97 phi_n_Der[d] = pphiDer[d] - pnDer[d]; 98 } 99 f0[0] = -poissonBracket(dim, psiDer, phi_n_Der) + poissonBracket(dim, logRefDenDer, ppsiDer); 100 if (u_t) f0[0] += u_t[PSI]; 101 } 102 103 static void f1_psi(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[]) { 104 const PetscScalar *ppsi = &u_x[uOff_x[PSI]]; 105 PetscInt d; 106 107 for (d = 0; d < dim - 1; ++d) f1[d] = -(s_ctx->eta / s_ctx->beta) * ppsi[d]; 108 } 109 110 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[]) { 111 f0[0] = -u[uOff[OMEGA]]; 112 } 113 114 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[]) { 115 const PetscScalar *pphi = &u_x[uOff_x[PHI]]; 116 PetscInt d; 117 118 for (d = 0; d < dim - 1; ++d) f1[d] = pphi[d]; 119 } 120 121 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[]) { 122 f0[0] = u[uOff[JZ]]; 123 } 124 125 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[]) { 126 const PetscScalar *ppsi = &u_x[uOff_x[PSI]]; 127 const PetscScalar *refPsiDer = &a_x[aOff_x[PSI]]; /* aOff_x[PSI] == 2*PSI */ 128 PetscInt d; 129 130 for (d = 0; d < dim - 1; ++d) f1[d] = ppsi[d] + refPsiDer[d]; 131 } 132 133 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) { 134 PetscFunctionBeginUser; 135 options->debug = 1; 136 options->plotRef = PETSC_FALSE; 137 options->implicit = PETSC_FALSE; 138 options->mu = 0; 139 options->eta = 0; 140 options->beta = 1; 141 options->a = 1; 142 options->b = PETSC_PI; 143 options->Jop = 0; 144 options->m = 1; 145 options->eps = 1.e-6; 146 147 PetscOptionsBegin(comm, "", "Poisson Problem Options", "DMPLEX"); 148 PetscCall(PetscOptionsInt("-debug", "The debugging level", "ex48.c", options->debug, &options->debug, NULL)); 149 PetscCall(PetscOptionsBool("-plot_ref", "Plot the reference fields", "ex48.c", options->plotRef, &options->plotRef, NULL)); 150 PetscCall(PetscOptionsReal("-mu", "mu", "ex48.c", options->mu, &options->mu, NULL)); 151 PetscCall(PetscOptionsReal("-eta", "eta", "ex48.c", options->eta, &options->eta, NULL)); 152 PetscCall(PetscOptionsReal("-beta", "beta", "ex48.c", options->beta, &options->beta, NULL)); 153 PetscCall(PetscOptionsReal("-Jop", "Jop", "ex48.c", options->Jop, &options->Jop, NULL)); 154 PetscCall(PetscOptionsReal("-m", "m", "ex48.c", options->m, &options->m, NULL)); 155 PetscCall(PetscOptionsReal("-eps", "eps", "ex48.c", options->eps, &options->eps, NULL)); 156 PetscCall(PetscOptionsBool("-implicit", "Use implicit time integrator", "ex48.c", options->implicit, &options->implicit, NULL)); 157 PetscOptionsEnd(); 158 options->ke = PetscSqrtScalar(options->Jop); 159 if (options->Jop == 0.0) { 160 options->Jo = 1.0 / PetscPowScalar(options->a, 2); 161 } else { 162 options->Jo = options->Jop * PetscCosReal(options->ke * options->a) / (1.0 - PetscCosReal(options->ke * options->a)); 163 } 164 options->ky = PETSC_PI * options->m / options->b; 165 if (PetscPowReal(options->ky, 2) < options->Jop) { 166 options->kx = PetscSqrtScalar(options->Jop - PetscPowScalar(options->ky, 2)); 167 options->DeltaPrime = -2.0 * options->kx * options->a * PetscCosReal(options->kx * options->a) / PetscSinReal(options->kx * options->a); 168 } else if (PetscPowReal(options->ky, 2) > options->Jop) { 169 options->kx = PetscSqrtScalar(PetscPowScalar(options->ky, 2) - options->Jop); 170 options->DeltaPrime = -2.0 * options->kx * options->a * PetscCoshReal(options->kx * options->a) / PetscSinhReal(options->kx * options->a); 171 } else { /*they're equal (or there's a NaN), lim(x*cot(x))_x->0=1*/ 172 options->kx = 0; 173 options->DeltaPrime = -2.0; 174 } 175 PetscCall(PetscPrintf(comm, "DeltaPrime=%g\n", (double)options->DeltaPrime)); 176 177 PetscFunctionReturn(0); 178 } 179 180 static void f_n(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) { 181 const PetscScalar *pn = &u[uOff[DENSITY]]; 182 *f0 = *pn; 183 } 184 185 static PetscErrorCode PostStep(TS ts) { 186 DM dm; 187 AppCtx *ctx; 188 PetscInt stepi, num; 189 Vec X; 190 191 PetscFunctionBeginUser; 192 PetscCall(TSGetApplicationContext(ts, &ctx)); 193 if (ctx->debug < 1) PetscFunctionReturn(0); 194 PetscCall(TSGetSolution(ts, &X)); 195 PetscCall(VecGetDM(X, &dm)); 196 PetscCall(TSGetStepNumber(ts, &stepi)); 197 PetscCall(DMGetOutputSequenceNumber(dm, &num, NULL)); 198 if (num < 0) PetscCall(DMSetOutputSequenceNumber(dm, 0, 0.0)); 199 PetscCall(PetscObjectSetName((PetscObject)X, "u")); 200 PetscCall(VecViewFromOptions(X, NULL, "-vec_view")); 201 /* print integrals */ 202 { 203 PetscDS prob; 204 DM plex; 205 PetscScalar den, tt[5]; 206 PetscCall(DMConvert(dm, DMPLEX, &plex)); 207 PetscCall(DMGetDS(plex, &prob)); 208 PetscCall(PetscDSSetObjective(prob, 0, &f_n)); 209 PetscCall(DMPlexComputeIntegralFEM(plex, X, tt, ctx)); 210 den = tt[0]; 211 PetscCall(DMDestroy(&plex)); 212 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "%" PetscInt_FMT ") total perturbed mass = %g\n", stepi, (double)PetscRealPart(den))); 213 } 214 PetscFunctionReturn(0); 215 } 216 217 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *ctx, DM *dm) { 218 PetscFunctionBeginUser; 219 PetscCall(DMCreate(comm, dm)); 220 PetscCall(DMSetType(*dm, DMPLEX)); 221 PetscCall(DMSetFromOptions(*dm)); 222 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); 223 224 PetscCall(DMGetBoundingBox(*dm, ctx->lower, ctx->upper)); 225 ctx->a = (ctx->upper[0] - ctx->lower[0]) / 2.0; 226 ctx->b = (ctx->upper[1] - ctx->lower[1]) / 2.0; 227 PetscFunctionReturn(0); 228 } 229 230 static PetscErrorCode log_n_0(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) { 231 AppCtx *lctx = (AppCtx *)ctx; 232 u[0] = 2. * lctx->a + x[0]; 233 return 0; 234 } 235 236 static PetscErrorCode Omega_0(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) { 237 u[0] = 0.0; 238 return 0; 239 } 240 241 static PetscErrorCode psi_0(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) { 242 AppCtx *lctx = (AppCtx *)ctx; 243 /* This sets up a symmetrix By flux aroound the mid point in x, which represents a current density flux along z. The stability 244 is analytically known and reported in ProcessOptions. */ 245 if (lctx->ke != 0.0) { 246 u[0] = (PetscCosReal(lctx->ke * (x[0] - lctx->a)) - PetscCosReal(lctx->ke * lctx->a)) / (1.0 - PetscCosReal(lctx->ke * lctx->a)); 247 } else { 248 u[0] = 1.0 - PetscPowScalar((x[0] - lctx->a) / lctx->a, 2); 249 } 250 return 0; 251 } 252 253 static PetscErrorCode initialSolution_n(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) { 254 u[0] = 0.0; 255 return 0; 256 } 257 258 static PetscErrorCode initialSolution_Omega(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) { 259 u[0] = 0.0; 260 return 0; 261 } 262 263 static PetscErrorCode initialSolution_psi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *a_ctx) { 264 AppCtx *ctx = (AppCtx *)a_ctx; 265 PetscScalar r = ctx->eps * (PetscScalar)(rand()) / (PetscScalar)(RAND_MAX); 266 if (x[0] == ctx->lower[0] || x[0] == ctx->upper[0]) r = 0; 267 u[0] = r; 268 /* PetscPrintf(PETSC_COMM_WORLD, "rand psi %lf\n",u[0]); */ 269 return 0; 270 } 271 272 static PetscErrorCode initialSolution_phi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) { 273 u[0] = 0.0; 274 return 0; 275 } 276 277 static PetscErrorCode initialSolution_jz(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) { 278 u[0] = 0.0; 279 return 0; 280 } 281 282 static PetscErrorCode SetupProblem(DM dm, AppCtx *ctx) { 283 PetscDS ds; 284 DMLabel label; 285 const PetscInt id = 1; 286 287 PetscFunctionBeginUser; 288 PetscCall(DMGetLabel(dm, "marker", &label)); 289 PetscCall(DMGetDS(dm, &ds)); 290 PetscCall(PetscDSSetResidual(ds, 0, f0_n, f1_n)); 291 PetscCall(PetscDSSetResidual(ds, 1, f0_Omega, f1_Omega)); 292 PetscCall(PetscDSSetResidual(ds, 2, f0_psi, f1_psi)); 293 PetscCall(PetscDSSetResidual(ds, 3, f0_phi, f1_phi)); 294 PetscCall(PetscDSSetResidual(ds, 4, f0_jz, f1_jz)); 295 ctx->initialFuncs[0] = initialSolution_n; 296 ctx->initialFuncs[1] = initialSolution_Omega; 297 ctx->initialFuncs[2] = initialSolution_psi; 298 ctx->initialFuncs[3] = initialSolution_phi; 299 ctx->initialFuncs[4] = initialSolution_jz; 300 for (PetscInt f = 0; f < 5; ++f) { 301 PetscCall(PetscDSSetImplicit(ds, f, ctx->implicit)); 302 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, f, 0, NULL, (void (*)(void))ctx->initialFuncs[f], NULL, ctx, NULL)); 303 } 304 PetscCall(PetscDSSetContext(ds, 0, ctx)); 305 PetscFunctionReturn(0); 306 } 307 308 static PetscErrorCode SetupEquilibriumFields(DM dm, DM dmAux, AppCtx *ctx) { 309 PetscErrorCode (*eqFuncs[3])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar[], void *) = {log_n_0, Omega_0, psi_0}; 310 Vec eq; 311 AppCtx *ctxarr[3]; 312 313 ctxarr[0] = ctxarr[1] = ctxarr[2] = ctx; /* each variable could have a different context */ 314 PetscFunctionBeginUser; 315 PetscCall(DMCreateLocalVector(dmAux, &eq)); 316 PetscCall(DMProjectFunctionLocal(dmAux, 0.0, eqFuncs, (void **)ctxarr, INSERT_ALL_VALUES, eq)); 317 PetscCall(DMSetAuxiliaryVec(dm, NULL, 0, 0, eq)); 318 if (ctx->plotRef) { /* plot reference functions */ 319 PetscViewer viewer = NULL; 320 PetscBool isHDF5, isVTK; 321 char buf[256]; 322 Vec global; 323 PetscInt dim; 324 325 PetscCall(DMGetDimension(dm, &dim)); 326 PetscCall(DMCreateGlobalVector(dmAux, &global)); 327 PetscCall(VecSet(global, .0)); /* BCs! */ 328 PetscCall(DMLocalToGlobalBegin(dmAux, eq, INSERT_VALUES, global)); 329 PetscCall(DMLocalToGlobalEnd(dmAux, eq, INSERT_VALUES, global)); 330 PetscCall(PetscViewerCreate(PetscObjectComm((PetscObject)dmAux), &viewer)); 331 #ifdef PETSC_HAVE_HDF5 332 PetscCall(PetscViewerSetType(viewer, PETSCVIEWERHDF5)); 333 #else 334 PetscCall(PetscViewerSetType(viewer, PETSCVIEWERVTK)); 335 #endif 336 PetscCall(PetscViewerSetFromOptions(viewer)); 337 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &isHDF5)); 338 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isVTK)); 339 if (isHDF5) { 340 PetscCall(PetscSNPrintf(buf, 256, "uEquilibrium-%" PetscInt_FMT "D.h5", dim)); 341 } else if (isVTK) { 342 PetscCall(PetscSNPrintf(buf, 256, "uEquilibrium-%" PetscInt_FMT "D.vtu", dim)); 343 PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_VTK_VTU)); 344 } 345 PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_WRITE)); 346 PetscCall(PetscViewerFileSetName(viewer, buf)); 347 if (isHDF5) PetscCall(DMView(dmAux, viewer)); 348 /* view equilibrium fields, this will overwrite fine grids with coarse grids! */ 349 PetscCall(PetscObjectSetName((PetscObject)global, "u0")); 350 PetscCall(VecView(global, viewer)); 351 PetscCall(PetscViewerDestroy(&viewer)); 352 PetscCall(VecDestroy(&global)); 353 } 354 PetscCall(VecDestroy(&eq)); 355 PetscFunctionReturn(0); 356 } 357 358 static PetscErrorCode SetupAuxDM(DM dm, PetscInt NfAux, PetscFE feAux[], AppCtx *user) { 359 DM dmAux, coordDM; 360 PetscInt f; 361 362 PetscFunctionBeginUser; 363 /* MUST call DMGetCoordinateDM() in order to get p4est setup if present */ 364 PetscCall(DMGetCoordinateDM(dm, &coordDM)); 365 if (!feAux) PetscFunctionReturn(0); 366 PetscCall(DMClone(dm, &dmAux)); 367 PetscCall(DMSetCoordinateDM(dmAux, coordDM)); 368 for (f = 0; f < NfAux; ++f) PetscCall(DMSetField(dmAux, f, NULL, (PetscObject)feAux[f])); 369 PetscCall(DMCreateDS(dmAux)); 370 PetscCall(SetupEquilibriumFields(dm, dmAux, user)); 371 PetscCall(DMDestroy(&dmAux)); 372 PetscFunctionReturn(0); 373 } 374 375 static PetscErrorCode SetupDiscretization(DM dm, AppCtx *ctx) { 376 DM cdm = dm; 377 PetscFE fe[5], feAux[3]; 378 PetscInt dim, Nf = 5, NfAux = 3, f; 379 PetscBool simplex; 380 MPI_Comm comm; 381 382 PetscFunctionBeginUser; 383 /* Create finite element */ 384 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 385 PetscCall(DMGetDimension(dm, &dim)); 386 PetscCall(DMPlexIsSimplex(dm, &simplex)); 387 PetscCall(PetscFECreateDefault(comm, dim, 1, simplex, NULL, -1, &fe[0])); 388 PetscCall(PetscObjectSetName((PetscObject)fe[0], "density")); 389 PetscCall(PetscFECreateDefault(comm, dim, 1, simplex, NULL, -1, &fe[1])); 390 PetscCall(PetscObjectSetName((PetscObject)fe[1], "vorticity")); 391 PetscCall(PetscFECopyQuadrature(fe[0], fe[1])); 392 PetscCall(PetscFECreateDefault(comm, dim, 1, simplex, NULL, -1, &fe[2])); 393 PetscCall(PetscObjectSetName((PetscObject)fe[2], "flux")); 394 PetscCall(PetscFECopyQuadrature(fe[0], fe[2])); 395 PetscCall(PetscFECreateDefault(comm, dim, 1, simplex, NULL, -1, &fe[3])); 396 PetscCall(PetscObjectSetName((PetscObject)fe[3], "potential")); 397 PetscCall(PetscFECopyQuadrature(fe[0], fe[3])); 398 PetscCall(PetscFECreateDefault(comm, dim, 1, simplex, NULL, -1, &fe[4])); 399 PetscCall(PetscObjectSetName((PetscObject)fe[4], "current")); 400 PetscCall(PetscFECopyQuadrature(fe[0], fe[4])); 401 402 PetscCall(PetscFECreateDefault(comm, dim, 1, simplex, NULL, -1, &feAux[0])); 403 PetscCall(PetscObjectSetName((PetscObject)feAux[0], "n_0")); 404 PetscCall(PetscFECopyQuadrature(fe[0], feAux[0])); 405 PetscCall(PetscFECreateDefault(comm, dim, 1, simplex, NULL, -1, &feAux[1])); 406 PetscCall(PetscObjectSetName((PetscObject)feAux[1], "vorticity_0")); 407 PetscCall(PetscFECopyQuadrature(fe[0], feAux[1])); 408 PetscCall(PetscFECreateDefault(comm, dim, 1, simplex, NULL, -1, &feAux[2])); 409 PetscCall(PetscObjectSetName((PetscObject)feAux[2], "flux_0")); 410 PetscCall(PetscFECopyQuadrature(fe[0], feAux[2])); 411 /* Set discretization and boundary conditions for each mesh */ 412 for (f = 0; f < Nf; ++f) PetscCall(DMSetField(dm, f, NULL, (PetscObject)fe[f])); 413 PetscCall(DMCreateDS(dm)); 414 PetscCall(SetupProblem(dm, ctx)); 415 while (cdm) { 416 PetscCall(SetupAuxDM(dm, NfAux, feAux, ctx)); 417 PetscCall(DMCopyDisc(dm, cdm)); 418 PetscCall(DMGetCoarseDM(cdm, &cdm)); 419 } 420 for (f = 0; f < Nf; ++f) PetscCall(PetscFEDestroy(&fe[f])); 421 for (f = 0; f < NfAux; ++f) PetscCall(PetscFEDestroy(&feAux[f])); 422 PetscFunctionReturn(0); 423 } 424 425 int main(int argc, char **argv) { 426 DM dm; 427 TS ts; 428 Vec u, r; 429 AppCtx ctx; 430 PetscReal t = 0.0; 431 PetscReal L2error = 0.0; 432 AppCtx *ctxarr[5]; 433 434 ctxarr[0] = ctxarr[1] = ctxarr[2] = ctxarr[3] = ctxarr[4] = &ctx; /* each variable could have a different context */ 435 s_ctx = &ctx; 436 PetscFunctionBeginUser; 437 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 438 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &ctx)); 439 /* create mesh and problem */ 440 PetscCall(CreateMesh(PETSC_COMM_WORLD, &ctx, &dm)); 441 PetscCall(DMSetApplicationContext(dm, &ctx)); 442 PetscCall(PetscMalloc1(5, &ctx.initialFuncs)); 443 PetscCall(SetupDiscretization(dm, &ctx)); 444 PetscCall(DMCreateGlobalVector(dm, &u)); 445 PetscCall(PetscObjectSetName((PetscObject)u, "u")); 446 PetscCall(VecDuplicate(u, &r)); 447 PetscCall(PetscObjectSetName((PetscObject)r, "r")); 448 /* create TS */ 449 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 450 PetscCall(TSSetDM(ts, dm)); 451 PetscCall(TSSetApplicationContext(ts, &ctx)); 452 PetscCall(DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &ctx)); 453 if (ctx.implicit) { 454 PetscCall(DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &ctx)); 455 PetscCall(DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &ctx)); 456 } else { 457 PetscCall(DMTSSetRHSFunctionLocal(dm, DMPlexTSComputeRHSFunctionFVM, &ctx)); 458 } 459 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); 460 PetscCall(TSSetFromOptions(ts)); 461 PetscCall(TSSetPostStep(ts, PostStep)); 462 /* make solution & solve */ 463 PetscCall(DMProjectFunction(dm, t, ctx.initialFuncs, (void **)ctxarr, INSERT_ALL_VALUES, u)); 464 PetscCall(TSSetSolution(ts, u)); 465 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view")); 466 PetscCall(PostStep(ts)); /* print the initial state */ 467 PetscCall(TSSolve(ts, u)); 468 PetscCall(TSGetTime(ts, &t)); 469 PetscCall(DMComputeL2Diff(dm, t, ctx.initialFuncs, (void **)ctxarr, u, &L2error)); 470 if (L2error < 1.0e-11) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: < 1.0e-11\n")); 471 else PetscCall(PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", (double)L2error)); 472 PetscCall(VecDestroy(&u)); 473 PetscCall(VecDestroy(&r)); 474 PetscCall(TSDestroy(&ts)); 475 PetscCall(DMDestroy(&dm)); 476 PetscCall(PetscFree(ctx.initialFuncs)); 477 PetscCall(PetscFinalize()); 478 return 0; 479 } 480 481 /*TEST 482 483 test: 484 suffix: 0 485 args: -debug 1 -dm_refine 1 -dm_plex_simplex 0 -dm_plex_box_faces 3,3 -dm_plex_box_bd periodic,none -dm_plex_box_upper 2.0,6.283185307179586 \ 486 -ts_max_steps 1 -ts_max_time 10. -ts_dt 1.0 487 test: 488 # Remapping with periodicity is broken 489 suffix: 1 490 args: -debug 1 -dm_plex_shape cylinder -dm_plex_dim 3 -dm_refine 1 -dm_refine_remap 0 -dm_plex_cylinder_bd periodic -dm_plex_boundary_label marker \ 491 -ts_max_steps 1 -ts_max_time 10. -ts_dt 1.0 492 493 TEST*/ 494