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