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 #include <assert.h>
21
22 typedef struct {
23 PetscInt debug; /* The debugging level */
24 PetscBool plotRef; /* Plot the reference fields */
25 /* Domain and mesh definition */
26 PetscInt dim; /* The topological mesh dimension */
27 char filename[2048]; /* The optional ExodusII file */
28 PetscBool cell_simplex; /* Simplicial mesh */
29 DMBoundaryType boundary_types[3];
30 PetscInt cells[3];
31 PetscInt refine;
32 /* geometry */
33 PetscReal domain_lo[3], domain_hi[3];
34 DMBoundaryType periodicity[3]; /* The domain periodicity */
35 PetscReal b0[3]; /* not used */
36 /* Problem definition */
37 PetscErrorCode (**initialFuncs)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
38 PetscReal mu, eta, beta;
39 PetscReal a,b,Jo,Jop,m,ke,kx,ky,DeltaPrime,eps;
40 /* solver */
41 PetscBool implicit;
42 } AppCtx;
43
44 static AppCtx *s_ctx;
45
poissonBracket(PetscInt dim,const PetscScalar df[],const PetscScalar dg[])46 static PetscScalar poissonBracket(PetscInt dim, const PetscScalar df[], const PetscScalar dg[])
47 {
48 PetscScalar ret = df[0]*dg[1] - df[1]*dg[0];
49 return ret;
50 }
51
52 enum field_idx {DENSITY,OMEGA,PSI,PHI,JZ};
53
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[])54 static void f0_n(PetscInt dim, PetscInt Nf, PetscInt NfAux,
55 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
56 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
57 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
58 {
59 const PetscScalar *pnDer = &u_x[uOff_x[DENSITY]];
60 const PetscScalar *ppsiDer = &u_x[uOff_x[PSI]];
61 const PetscScalar *pphiDer = &u_x[uOff_x[PHI]];
62 const PetscScalar *jzDer = &u_x[uOff_x[JZ]];
63 const PetscScalar *logRefDenDer = &a_x[aOff_x[DENSITY]];
64 f0[0] += - poissonBracket(dim,pnDer, pphiDer) - s_ctx->beta*poissonBracket(dim,jzDer, ppsiDer) - poissonBracket(dim,logRefDenDer, pphiDer);
65 if (u_t) f0[0] += u_t[DENSITY];
66 }
67
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[])68 static void f1_n(PetscInt dim, PetscInt Nf, PetscInt NfAux,
69 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
70 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
71 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
72 {
73 const PetscScalar *pnDer = &u_x[uOff_x[DENSITY]];
74 PetscInt d;
75
76 for (d = 0; d < dim-1; ++d) f1[d] = -s_ctx->mu*pnDer[d];
77 }
78
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[])79 static void f0_Omega(PetscInt dim, PetscInt Nf, PetscInt NfAux,
80 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
81 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
82 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
83 {
84 const PetscScalar *pOmegaDer = &u_x[uOff_x[OMEGA]];
85 const PetscScalar *ppsiDer = &u_x[uOff_x[PSI]];
86 const PetscScalar *pphiDer = &u_x[uOff_x[PHI]];
87 const PetscScalar *jzDer = &u_x[uOff_x[JZ]];
88
89 f0[0] += - poissonBracket(dim,pOmegaDer, pphiDer) - s_ctx->beta*poissonBracket(dim,jzDer, ppsiDer);
90 if (u_t) f0[0] += u_t[OMEGA];
91 }
92
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[])93 static void f1_Omega(PetscInt dim, PetscInt Nf, PetscInt NfAux,
94 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
95 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
96 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
97 {
98 const PetscScalar *pOmegaDer = &u_x[uOff_x[OMEGA]];
99 PetscInt d;
100
101 for (d = 0; d < dim-1; ++d) f1[d] = -s_ctx->mu*pOmegaDer[d];
102 }
103
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[])104 static void f0_psi(PetscInt dim, PetscInt Nf, PetscInt NfAux,
105 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
106 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
107 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
108 {
109 const PetscScalar *pnDer = &u_x[uOff_x[DENSITY]];
110 const PetscScalar *ppsiDer = &u_x[uOff_x[PSI]];
111 const PetscScalar *pphiDer = &u_x[uOff_x[PHI]];
112 const PetscScalar *refPsiDer = &a_x[aOff_x[PSI]];
113 const PetscScalar *logRefDenDer= &a_x[aOff_x[DENSITY]];
114 PetscScalar psiDer[3];
115 PetscScalar phi_n_Der[3];
116 PetscInt d;
117 if (dim < 2) {MPI_Abort(MPI_COMM_WORLD,1); return;} /* this is needed so that the clang static analyzer does not generate a warning about variables used by not set */
118 for (d = 0; d < dim; ++d) {
119 psiDer[d] = refPsiDer[d] + ppsiDer[d];
120 phi_n_Der[d] = pphiDer[d] - pnDer[d];
121 }
122 f0[0] = - poissonBracket(dim,psiDer, phi_n_Der) + poissonBracket(dim,logRefDenDer, ppsiDer);
123 if (u_t) f0[0] += u_t[PSI];
124 }
125
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[])126 static void f1_psi(PetscInt dim, PetscInt Nf, PetscInt NfAux,
127 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
128 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
129 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
130 {
131 const PetscScalar *ppsi = &u_x[uOff_x[PSI]];
132 PetscInt d;
133
134 for (d = 0; d < dim-1; ++d) f1[d] = -(s_ctx->eta/s_ctx->beta)*ppsi[d];
135 }
136
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[])137 static void f0_phi(PetscInt dim, PetscInt Nf, PetscInt NfAux,
138 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
139 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
140 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
141 {
142 f0[0] = -u[uOff[OMEGA]];
143 }
144
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[])145 static void f1_phi(PetscInt dim, PetscInt Nf, PetscInt NfAux,
146 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
147 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
148 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
149 {
150 const PetscScalar *pphi = &u_x[uOff_x[PHI]];
151 PetscInt d;
152
153 for (d = 0; d < dim-1; ++d) f1[d] = pphi[d];
154 }
155
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[])156 static void f0_jz(PetscInt dim, PetscInt Nf, PetscInt NfAux,
157 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
158 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
159 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
160 {
161 f0[0] = u[uOff[JZ]];
162 }
163
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[])164 static void f1_jz(PetscInt dim, PetscInt Nf, PetscInt NfAux,
165 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
166 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
167 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
168 {
169 const PetscScalar *ppsi = &u_x[uOff_x[PSI]];
170 const PetscScalar *refPsiDer = &a_x[aOff_x[PSI]]; /* aOff_x[PSI] == 2*PSI */
171 PetscInt d;
172
173 for (d = 0; d < dim-1; ++d) f1[d] = ppsi[d] + refPsiDer[d];
174 }
175
ProcessOptions(MPI_Comm comm,AppCtx * options)176 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
177 {
178 PetscBool flg;
179 PetscErrorCode ierr;
180 PetscInt ii, bd;
181 PetscFunctionBeginUser;
182 options->debug = 1;
183 options->plotRef = PETSC_FALSE;
184 options->dim = 2;
185 options->filename[0] = '\0';
186 options->cell_simplex = PETSC_FALSE;
187 options->implicit = PETSC_FALSE;
188 options->refine = 2;
189 options->domain_lo[0] = 0.0;
190 options->domain_lo[1] = 0.0;
191 options->domain_lo[2] = 0.0;
192 options->domain_hi[0] = 2.0;
193 options->domain_hi[1] = 2.0*PETSC_PI;
194 options->domain_hi[2] = 2.0;
195 options->periodicity[0] = DM_BOUNDARY_NONE;
196 options->periodicity[1] = DM_BOUNDARY_NONE;
197 options->periodicity[2] = DM_BOUNDARY_NONE;
198 options->mu = 0;
199 options->eta = 0;
200 options->beta = 1;
201 options->a = 1;
202 options->b = PETSC_PI;
203 options->Jop = 0;
204 options->m = 1;
205 options->eps = 1.e-6;
206
207 for (ii = 0; ii < options->dim; ++ii) options->cells[ii] = 4;
208 ierr = PetscOptionsBegin(comm, "", "Poisson Problem Options", "DMPLEX");CHKERRQ(ierr);
209 ierr = PetscOptionsInt("-debug", "The debugging level", "ex48.c", options->debug, &options->debug, NULL);CHKERRQ(ierr);
210 ierr = PetscOptionsBool("-plot_ref", "Plot the reference fields", "ex48.c", options->plotRef, &options->plotRef, NULL);CHKERRQ(ierr);
211 ierr = PetscOptionsInt("-dim", "The topological mesh dimension", "ex48.c", options->dim, &options->dim, NULL);CHKERRQ(ierr);
212 if (options->dim < 2 || options->dim > 3) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Dim %D must be 2 or 3",options->dim);CHKERRQ(ierr);
213 ierr = PetscOptionsInt("-dm_refine", "Hack to get refinement level for cylinder", "ex48.c", options->refine, &options->refine, NULL);CHKERRQ(ierr);
214 ierr = PetscOptionsReal("-mu", "mu", "ex48.c", options->mu, &options->mu, NULL);CHKERRQ(ierr);
215 ierr = PetscOptionsReal("-eta", "eta", "ex48.c", options->eta, &options->eta, NULL);CHKERRQ(ierr);
216 ierr = PetscOptionsReal("-beta", "beta", "ex48.c", options->beta, &options->beta, NULL);CHKERRQ(ierr);
217 ierr = PetscOptionsReal("-Jop", "Jop", "ex48.c", options->Jop, &options->Jop, NULL);CHKERRQ(ierr);
218 ierr = PetscOptionsReal("-m", "m", "ex48.c", options->m, &options->m, NULL);CHKERRQ(ierr);
219 ierr = PetscOptionsReal("-eps", "eps", "ex48.c", options->eps, &options->eps, NULL);CHKERRQ(ierr);
220 ierr = PetscOptionsString("-f", "Exodus.II filename to read", "ex48.c", options->filename, options->filename, sizeof(options->filename), &flg);CHKERRQ(ierr);
221 ierr = PetscOptionsBool("-cell_simplex", "Simplicial (true) or tensor (false) mesh", "ex48.c", options->cell_simplex, &options->cell_simplex, NULL);CHKERRQ(ierr);
222 ierr = PetscOptionsBool("-implicit", "Use implicit time integrator", "ex48.c", options->implicit, &options->implicit, NULL);CHKERRQ(ierr);
223 ii = options->dim;
224 ierr = PetscOptionsRealArray("-domain_hi", "Domain size", "ex48.c", options->domain_hi, &ii, NULL);CHKERRQ(ierr);
225 ii = options->dim;
226 ierr = PetscOptionsRealArray("-domain_lo", "Domain size", "ex48.c", options->domain_lo, &ii, NULL);CHKERRQ(ierr);
227 ii = options->dim;
228 bd = options->periodicity[0];
229 ierr = PetscOptionsEList("-x_periodicity", "The x-boundary periodicity", "ex48.c", DMBoundaryTypes, 5, DMBoundaryTypes[options->periodicity[0]], &bd, NULL);CHKERRQ(ierr);
230 options->periodicity[0] = (DMBoundaryType) bd;
231 bd = options->periodicity[1];
232 ierr = PetscOptionsEList("-y_periodicity", "The y-boundary periodicity", "ex48.c", DMBoundaryTypes, 5, DMBoundaryTypes[options->periodicity[1]], &bd, NULL);CHKERRQ(ierr);
233 options->periodicity[1] = (DMBoundaryType) bd;
234 bd = options->periodicity[2];
235 ierr = PetscOptionsEList("-z_periodicity", "The z-boundary periodicity", "ex48.c", DMBoundaryTypes, 5, DMBoundaryTypes[options->periodicity[2]], &bd, NULL);CHKERRQ(ierr);
236 options->periodicity[2] = (DMBoundaryType) bd;
237 ii = options->dim;
238 ierr = PetscOptionsIntArray("-cells", "Number of cells in each dimension", "ex48.c", options->cells, &ii, NULL);CHKERRQ(ierr);
239 ierr = PetscOptionsEnd();CHKERRQ(ierr);
240 options->a = (options->domain_hi[0]-options->domain_lo[0])/2.0;
241 options->b = (options->domain_hi[1]-options->domain_lo[1])/2.0;
242 for (ii = 0; ii < options->dim; ++ii) {
243 if (options->domain_hi[ii] <= options->domain_lo[ii]) SETERRQ3(comm,PETSC_ERR_ARG_WRONG,"Domain %D lo=%g hi=%g",ii,options->domain_lo[ii],options->domain_hi[ii]);
244 }
245 options->ke = PetscSqrtScalar(options->Jop);
246 if (options->Jop==0.0) {
247 options->Jo = 1.0/PetscPowScalar(options->a,2);
248 } else {
249 options->Jo = options->Jop*PetscCosReal(options->ke*options->a)/(1.0-PetscCosReal(options->ke*options->a));
250 }
251 options->ky = PETSC_PI*options->m/options->b;
252 if (PetscPowReal(options->ky, 2) < options->Jop) {
253 options->kx = PetscSqrtScalar(options->Jop-PetscPowScalar(options->ky,2));
254 options->DeltaPrime = -2.0*options->kx*options->a*PetscCosReal(options->kx*options->a)/PetscSinReal(options->kx*options->a);
255 } else if (PetscPowReal(options->ky, 2) > options->Jop) {
256 options->kx = PetscSqrtScalar(PetscPowScalar(options->ky,2)-options->Jop);
257 options->DeltaPrime = -2.0*options->kx*options->a*PetscCoshReal(options->kx*options->a)/PetscSinhReal(options->kx*options->a);
258 } else { /*they're equal (or there's a NaN), lim(x*cot(x))_x->0=1*/
259 options->kx = 0;
260 options->DeltaPrime = -2.0;
261 }
262 ierr = PetscPrintf(comm, "DeltaPrime=%g\n",options->DeltaPrime);CHKERRQ(ierr);
263
264 PetscFunctionReturn(0);
265 }
266
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)267 static void f_n(PetscInt dim, PetscInt Nf, PetscInt NfAux,
268 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
269 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
270 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar *f0)
271 {
272 const PetscScalar *pn = &u[uOff[DENSITY]];
273 *f0 = *pn;
274 }
275
PostStep(TS ts)276 static PetscErrorCode PostStep(TS ts)
277 {
278 PetscErrorCode ierr;
279 DM dm;
280 AppCtx *ctx;
281 PetscInt stepi,num;
282 Vec X;
283 PetscFunctionBegin;
284 ierr = TSGetApplicationContext(ts, &ctx);CHKERRQ(ierr); assert(ctx);
285 if (ctx->debug<1) PetscFunctionReturn(0);
286 ierr = TSGetSolution(ts, &X);CHKERRQ(ierr);
287 ierr = VecGetDM(X, &dm);CHKERRQ(ierr);
288 ierr = TSGetStepNumber(ts, &stepi);CHKERRQ(ierr);
289 ierr = DMGetOutputSequenceNumber(dm, &num, NULL);CHKERRQ(ierr);
290 if (num < 0) {ierr = DMSetOutputSequenceNumber(dm, 0, 0.0);CHKERRQ(ierr);}
291 ierr = PetscObjectSetName((PetscObject) X, "u");CHKERRQ(ierr);
292 ierr = VecViewFromOptions(X, NULL, "-vec_view");CHKERRQ(ierr);
293 /* print integrals */
294 {
295 PetscDS prob;
296 DM plex;
297 PetscScalar den, tt[5];
298 ierr = DMConvert(dm, DMPLEX, &plex);CHKERRQ(ierr);
299 ierr = DMGetDS(plex, &prob);CHKERRQ(ierr);
300 ierr = PetscDSSetObjective(prob, 0, &f_n);CHKERRQ(ierr);
301 ierr = DMPlexComputeIntegralFEM(plex,X,tt,ctx);CHKERRQ(ierr);
302 den = tt[0];
303 ierr = DMDestroy(&plex);CHKERRQ(ierr);
304 PetscPrintf(PetscObjectComm((PetscObject)dm), "%D) total perturbed mass = %g\n", stepi, (double) PetscRealPart(den));CHKERRQ(ierr);
305 }
306 PetscFunctionReturn(0);
307 }
308
CreateBCLabel(DM dm,const char name[])309 static PetscErrorCode CreateBCLabel(DM dm, const char name[])
310 {
311 DM plex;
312 DMLabel label;
313 PetscErrorCode ierr;
314
315 PetscFunctionBeginUser;
316 ierr = DMCreateLabel(dm, name);CHKERRQ(ierr);
317 ierr = DMGetLabel(dm, name, &label);CHKERRQ(ierr);
318 ierr = DMConvert(dm, DMPLEX, &plex);CHKERRQ(ierr);
319 ierr = DMPlexMarkBoundaryFaces(dm, 1, label);CHKERRQ(ierr);
320 ierr = DMDestroy(&plex);CHKERRQ(ierr);
321 PetscFunctionReturn(0);
322 }
323
CreateMesh(MPI_Comm comm,AppCtx * ctx,DM * dm)324 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *ctx, DM *dm)
325 {
326 PetscInt dim = ctx->dim;
327 const char *filename = ctx->filename;
328 size_t len;
329 PetscMPIInt numProcs;
330 PetscErrorCode ierr;
331
332 PetscFunctionBeginUser;
333 ierr = MPI_Comm_size(comm, &numProcs);CHKERRQ(ierr);
334 ierr = PetscStrlen(filename, &len);CHKERRQ(ierr);
335 if (len) {
336 ierr = DMPlexCreateFromFile(comm, filename, PETSC_TRUE, dm);CHKERRQ(ierr);
337 } else {
338 PetscInt d;
339
340 /* create DM */
341 if (ctx->cell_simplex && dim == 3) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Cannot mesh a cylinder with simplices");
342 if (dim==2) {
343 PetscInt refineRatio, totCells = 1;
344 if (ctx->cell_simplex) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Cannot mesh 2D with simplices");
345 refineRatio = PetscMax((PetscInt) (PetscPowReal(numProcs, 1.0/dim) + 0.1) - 1, 1);
346 for (d = 0; d < dim; ++d) {
347 if (ctx->cells[d] < refineRatio) ctx->cells[d] = refineRatio;
348 if (ctx->periodicity[d]==DM_BOUNDARY_PERIODIC && ctx->cells[d]*refineRatio <= 2) refineRatio = 2;
349 }
350 for (d = 0; d < dim; ++d) {
351 ctx->cells[d] *= refineRatio;
352 totCells *= ctx->cells[d];
353 }
354 if (totCells % numProcs) SETERRQ2(comm,PETSC_ERR_ARG_WRONG,"Total cells %D not divisible by processes %D", totCells, numProcs);
355 ierr = DMPlexCreateBoxMesh(comm, dim, PETSC_FALSE, ctx->cells, ctx->domain_lo, ctx->domain_hi, ctx->periodicity, PETSC_TRUE, dm);CHKERRQ(ierr);
356 } else {
357 if (ctx->periodicity[0]==DM_BOUNDARY_PERIODIC || ctx->periodicity[1]==DM_BOUNDARY_PERIODIC) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Cannot do periodic in x or y in a cylinder");
358 /* we stole dm_refine so clear it */
359 ierr = PetscOptionsClearValue(NULL,"-dm_refine");CHKERRQ(ierr);
360 ierr = DMPlexCreateHexCylinderMesh(comm, ctx->refine, ctx->periodicity[2], dm);CHKERRQ(ierr);
361 }
362 }
363 {
364 DM distributedMesh = NULL;
365 /* Distribute mesh over processes */
366 ierr = DMPlexDistribute(*dm, 0, NULL, &distributedMesh);CHKERRQ(ierr);
367 if (distributedMesh) {
368 ierr = DMDestroy(dm);CHKERRQ(ierr);
369 *dm = distributedMesh;
370 }
371 }
372 {
373 PetscBool hasLabel;
374 ierr = DMHasLabel(*dm, "marker", &hasLabel);CHKERRQ(ierr);
375 if (!hasLabel) {ierr = CreateBCLabel(*dm, "marker");CHKERRQ(ierr);}
376 }
377 {
378 char convType[256];
379 PetscBool flg;
380 ierr = PetscOptionsBegin(comm, "", "Mesh conversion options", "DMPLEX");CHKERRQ(ierr);
381 ierr = PetscOptionsFList("-dm_plex_convert_type","Convert DMPlex to another format","ex48",DMList,DMPLEX,convType,256,&flg);CHKERRQ(ierr);
382 ierr = PetscOptionsEnd();CHKERRQ(ierr);
383 if (flg) {
384 DM dmConv;
385 ierr = DMConvert(*dm,convType,&dmConv);CHKERRQ(ierr);
386 if (dmConv) {
387 ierr = DMDestroy(dm);CHKERRQ(ierr);
388 *dm = dmConv;
389 }
390 }
391 }
392 ierr = PetscObjectSetName((PetscObject) *dm, "Mesh");CHKERRQ(ierr);
393 ierr = DMSetFromOptions(*dm);CHKERRQ(ierr);
394 ierr = DMLocalizeCoordinates(*dm);CHKERRQ(ierr); /* needed for periodic */
395 PetscFunctionReturn(0);
396 }
397
log_n_0(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nf,PetscScalar * u,void * ctx)398 static PetscErrorCode log_n_0(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
399 {
400 AppCtx *lctx = (AppCtx*)ctx;
401 assert(ctx);
402 u[0] = (lctx->domain_hi-lctx->domain_lo)+x[0];
403 return 0;
404 }
405
Omega_0(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nf,PetscScalar * u,void * ctx)406 static PetscErrorCode Omega_0(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
407 {
408 u[0] = 0.0;
409 return 0;
410 }
411
psi_0(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nf,PetscScalar * u,void * ctx)412 static PetscErrorCode psi_0(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
413 {
414 AppCtx *lctx = (AppCtx*)ctx;
415 assert(ctx);
416 /* This sets up a symmetrix By flux aroound the mid point in x, which represents a current density flux along z. The stability
417 is analytically known and reported in ProcessOptions. */
418 if (lctx->ke!=0.0) {
419 u[0] = (PetscCosReal(lctx->ke*(x[0]-lctx->a))-PetscCosReal(lctx->ke*lctx->a))/(1.0-PetscCosReal(lctx->ke*lctx->a));
420 } else {
421 u[0] = 1.0-PetscPowScalar((x[0]-lctx->a)/lctx->a,2);
422 }
423 return 0;
424 }
425
initialSolution_n(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nf,PetscScalar * u,void * ctx)426 static PetscErrorCode initialSolution_n(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
427 {
428 u[0] = 0.0;
429 return 0;
430 }
431
initialSolution_Omega(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nf,PetscScalar * u,void * ctx)432 static PetscErrorCode initialSolution_Omega(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
433 {
434 u[0] = 0.0;
435 return 0;
436 }
437
initialSolution_psi(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nf,PetscScalar * u,void * a_ctx)438 static PetscErrorCode initialSolution_psi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *a_ctx)
439 {
440 AppCtx *ctx = (AppCtx*)a_ctx;
441 PetscScalar r = ctx->eps*(PetscScalar) (rand()) / (PetscScalar) (RAND_MAX);
442 assert(ctx);
443 if (x[0] == ctx->domain_lo[0] || x[0] == ctx->domain_hi[0]) r = 0;
444 u[0] = r;
445 /* PetscPrintf(PETSC_COMM_WORLD, "rand psi %lf\n",u[0]); */
446 return 0;
447 }
448
initialSolution_phi(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nf,PetscScalar * u,void * ctx)449 static PetscErrorCode initialSolution_phi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
450 {
451 u[0] = 0.0;
452 return 0;
453 }
454
initialSolution_jz(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nf,PetscScalar * u,void * ctx)455 static PetscErrorCode initialSolution_jz(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
456 {
457 u[0] = 0.0;
458 return 0;
459 }
460
SetupProblem(DM dm,AppCtx * ctx)461 static PetscErrorCode SetupProblem(DM dm, AppCtx *ctx)
462 {
463 PetscDS prob;
464 const PetscInt id = 1;
465 PetscErrorCode ierr, f;
466
467 PetscFunctionBeginUser;
468 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr);
469 ierr = PetscDSSetResidual(prob, 0, f0_n, f1_n);CHKERRQ(ierr);
470 ierr = PetscDSSetResidual(prob, 1, f0_Omega, f1_Omega);CHKERRQ(ierr);
471 ierr = PetscDSSetResidual(prob, 2, f0_psi, f1_psi);CHKERRQ(ierr);
472 ierr = PetscDSSetResidual(prob, 3, f0_phi, f1_phi);CHKERRQ(ierr);
473 ierr = PetscDSSetResidual(prob, 4, f0_jz, f1_jz);CHKERRQ(ierr);
474 ctx->initialFuncs[0] = initialSolution_n;
475 ctx->initialFuncs[1] = initialSolution_Omega;
476 ctx->initialFuncs[2] = initialSolution_psi;
477 ctx->initialFuncs[3] = initialSolution_phi;
478 ctx->initialFuncs[4] = initialSolution_jz;
479 for (f = 0; f < 5; ++f) {
480 ierr = PetscDSSetImplicit(prob, f, ctx->implicit);CHKERRQ(ierr);
481 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", "marker", f, 0, NULL, (void (*)(void)) ctx->initialFuncs[f], NULL, 1, &id, ctx);CHKERRQ(ierr);
482 }
483 ierr = PetscDSSetContext(prob, 0, ctx);CHKERRQ(ierr);
484 PetscFunctionReturn(0);
485 }
486
SetupEquilibriumFields(DM dm,DM dmAux,AppCtx * ctx)487 static PetscErrorCode SetupEquilibriumFields(DM dm, DM dmAux, AppCtx *ctx)
488 {
489 PetscErrorCode (*eqFuncs[3])(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar [], void *) = {log_n_0, Omega_0, psi_0};
490 Vec eq;
491 PetscErrorCode ierr;
492 AppCtx *ctxarr[3];
493
494 ctxarr[0] = ctxarr[1] = ctxarr[2] = ctx; /* each variable could have a different context */
495 PetscFunctionBegin;
496 ierr = DMCreateLocalVector(dmAux, &eq);CHKERRQ(ierr);
497 ierr = DMProjectFunctionLocal(dmAux, 0.0, eqFuncs, (void **)ctxarr, INSERT_ALL_VALUES, eq);CHKERRQ(ierr);
498 ierr = PetscObjectCompose((PetscObject) dm, "A", (PetscObject) eq);CHKERRQ(ierr);
499 if (ctx->plotRef) { /* plot reference functions */
500 PetscViewer viewer = NULL;
501 PetscBool isHDF5,isVTK;
502 char buf[256];
503 Vec global;
504 ierr = DMCreateGlobalVector(dmAux,&global);CHKERRQ(ierr);
505 ierr = VecSet(global,.0);CHKERRQ(ierr); /* BCs! */
506 ierr = DMLocalToGlobalBegin(dmAux,eq,INSERT_VALUES,global);CHKERRQ(ierr);
507 ierr = DMLocalToGlobalEnd(dmAux,eq,INSERT_VALUES,global);CHKERRQ(ierr);
508 ierr = PetscViewerCreate(PetscObjectComm((PetscObject)dmAux),&viewer);CHKERRQ(ierr);
509 #ifdef PETSC_HAVE_HDF5
510 ierr = PetscViewerSetType(viewer,PETSCVIEWERHDF5);CHKERRQ(ierr);
511 #else
512 ierr = PetscViewerSetType(viewer,PETSCVIEWERVTK);CHKERRQ(ierr);
513 #endif
514 ierr = PetscViewerSetFromOptions(viewer);CHKERRQ(ierr);
515 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&isHDF5);CHKERRQ(ierr);
516 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERVTK,&isVTK);CHKERRQ(ierr);
517 if (isHDF5) {
518 ierr = PetscSNPrintf(buf, 256, "uEquilibrium-%dD.h5", ctx->dim);CHKERRQ(ierr);
519 } else if (isVTK) {
520 ierr = PetscSNPrintf(buf, 256, "uEquilibrium-%dD.vtu", ctx->dim);CHKERRQ(ierr);
521 ierr = PetscViewerPushFormat(viewer,PETSC_VIEWER_VTK_VTU);CHKERRQ(ierr);
522 }
523 ierr = PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);CHKERRQ(ierr);
524 ierr = PetscViewerFileSetName(viewer,buf);CHKERRQ(ierr);
525 if (isHDF5) {ierr = DMView(dmAux,viewer);CHKERRQ(ierr);}
526 /* view equilibrium fields, this will overwrite fine grids with coarse grids! */
527 ierr = PetscObjectSetName((PetscObject) global, "u0");CHKERRQ(ierr);
528 ierr = VecView(global,viewer);CHKERRQ(ierr);
529 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
530 ierr = VecDestroy(&global);CHKERRQ(ierr);
531 }
532 ierr = VecDestroy(&eq);CHKERRQ(ierr);
533 PetscFunctionReturn(0);
534 }
535
SetupAuxDM(DM dm,PetscInt NfAux,PetscFE feAux[],AppCtx * user)536 static PetscErrorCode SetupAuxDM(DM dm, PetscInt NfAux, PetscFE feAux[], AppCtx *user)
537 {
538 DM dmAux, coordDM;
539 PetscInt f;
540 PetscErrorCode ierr;
541
542 PetscFunctionBegin;
543 /* MUST call DMGetCoordinateDM() in order to get p4est setup if present */
544 ierr = DMGetCoordinateDM(dm, &coordDM);CHKERRQ(ierr);
545 if (!feAux) PetscFunctionReturn(0);
546 ierr = DMClone(dm, &dmAux);CHKERRQ(ierr);
547 ierr = PetscObjectCompose((PetscObject) dm, "dmAux", (PetscObject) dmAux);CHKERRQ(ierr);
548 ierr = DMSetCoordinateDM(dmAux, coordDM);CHKERRQ(ierr);
549 for (f = 0; f < NfAux; ++f) {ierr = DMSetField(dmAux, f, NULL, (PetscObject) feAux[f]);CHKERRQ(ierr);}
550 ierr = DMCreateDS(dmAux);CHKERRQ(ierr);
551 ierr = SetupEquilibriumFields(dm, dmAux, user);CHKERRQ(ierr);
552 ierr = DMDestroy(&dmAux);CHKERRQ(ierr);
553 PetscFunctionReturn(0);
554 }
555
SetupDiscretization(DM dm,AppCtx * ctx)556 static PetscErrorCode SetupDiscretization(DM dm, AppCtx *ctx)
557 {
558 DM cdm = dm;
559 const PetscInt dim = ctx->dim;
560 PetscFE fe[5], feAux[3];
561 PetscInt Nf = 5, NfAux = 3, f;
562 PetscBool cell_simplex = ctx->cell_simplex;
563 MPI_Comm comm;
564 PetscErrorCode ierr;
565
566 PetscFunctionBeginUser;
567 /* Create finite element */
568 ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr);
569 ierr = PetscFECreateDefault(comm, dim, 1, cell_simplex, NULL, -1, &fe[0]);CHKERRQ(ierr);
570 ierr = PetscObjectSetName((PetscObject) fe[0], "density");CHKERRQ(ierr);
571 ierr = PetscFECreateDefault(comm, dim, 1, cell_simplex, NULL, -1, &fe[1]);CHKERRQ(ierr);
572 ierr = PetscObjectSetName((PetscObject) fe[1], "vorticity");CHKERRQ(ierr);
573 ierr = PetscFECopyQuadrature(fe[0], fe[1]);CHKERRQ(ierr);
574 ierr = PetscFECreateDefault(comm, dim, 1, cell_simplex, NULL, -1, &fe[2]);CHKERRQ(ierr);
575 ierr = PetscObjectSetName((PetscObject) fe[2], "flux");CHKERRQ(ierr);
576 ierr = PetscFECopyQuadrature(fe[0], fe[2]);CHKERRQ(ierr);
577 ierr = PetscFECreateDefault(comm, dim, 1, cell_simplex, NULL, -1, &fe[3]);CHKERRQ(ierr);
578 ierr = PetscObjectSetName((PetscObject) fe[3], "potential");CHKERRQ(ierr);
579 ierr = PetscFECopyQuadrature(fe[0], fe[3]);CHKERRQ(ierr);
580 ierr = PetscFECreateDefault(comm, dim, 1, cell_simplex, NULL, -1, &fe[4]);CHKERRQ(ierr);
581 ierr = PetscObjectSetName((PetscObject) fe[4], "current");CHKERRQ(ierr);
582 ierr = PetscFECopyQuadrature(fe[0], fe[4]);CHKERRQ(ierr);
583
584 ierr = PetscFECreateDefault(comm, dim, 1, cell_simplex, NULL, -1, &feAux[0]);CHKERRQ(ierr);
585 ierr = PetscObjectSetName((PetscObject) feAux[0], "n_0");CHKERRQ(ierr);
586 ierr = PetscFECopyQuadrature(fe[0], feAux[0]);CHKERRQ(ierr);
587 ierr = PetscFECreateDefault(comm, dim, 1, cell_simplex, NULL, -1, &feAux[1]);CHKERRQ(ierr);
588 ierr = PetscObjectSetName((PetscObject) feAux[1], "vorticity_0");CHKERRQ(ierr);
589 ierr = PetscFECopyQuadrature(fe[0], feAux[1]);CHKERRQ(ierr);
590 ierr = PetscFECreateDefault(comm, dim, 1, cell_simplex, NULL, -1, &feAux[2]);CHKERRQ(ierr);
591 ierr = PetscObjectSetName((PetscObject) feAux[2], "flux_0");CHKERRQ(ierr);
592 ierr = PetscFECopyQuadrature(fe[0], feAux[2]);CHKERRQ(ierr);
593 /* Set discretization and boundary conditions for each mesh */
594 for (f = 0; f < Nf; ++f) {ierr = DMSetField(dm, f, NULL, (PetscObject) fe[f]);CHKERRQ(ierr);}
595 ierr = DMCreateDS(dm);CHKERRQ(ierr);
596 ierr = SetupProblem(dm, ctx);CHKERRQ(ierr);
597 while (cdm) {
598 ierr = SetupAuxDM(dm, NfAux, feAux, ctx);CHKERRQ(ierr);
599 {
600 PetscBool hasLabel;
601
602 ierr = DMHasLabel(cdm, "marker", &hasLabel);CHKERRQ(ierr);
603 if (!hasLabel) {ierr = CreateBCLabel(cdm, "marker");CHKERRQ(ierr);}
604 }
605 ierr = DMCopyDisc(dm, cdm);CHKERRQ(ierr);
606 ierr = DMGetCoarseDM(cdm, &cdm);CHKERRQ(ierr);
607 }
608 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&fe[f]);CHKERRQ(ierr);}
609 for (f = 0; f < NfAux; ++f) {ierr = PetscFEDestroy(&feAux[f]);CHKERRQ(ierr);}
610 PetscFunctionReturn(0);
611 }
612
main(int argc,char ** argv)613 int main(int argc, char **argv)
614 {
615 DM dm;
616 TS ts;
617 Vec u, r;
618 AppCtx ctx;
619 PetscReal t = 0.0;
620 PetscReal L2error = 0.0;
621 PetscErrorCode ierr;
622 AppCtx *ctxarr[5];
623
624 ctxarr[0] = ctxarr[1] = ctxarr[2] = ctxarr[3] = ctxarr[4] = &ctx; /* each variable could have a different context */
625 s_ctx = &ctx;
626 ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr;
627 ierr = ProcessOptions(PETSC_COMM_WORLD, &ctx);CHKERRQ(ierr);
628 /* create mesh and problem */
629 ierr = CreateMesh(PETSC_COMM_WORLD, &ctx, &dm);CHKERRQ(ierr);
630 ierr = DMSetApplicationContext(dm, &ctx);CHKERRQ(ierr);
631 ierr = PetscMalloc1(5, &ctx.initialFuncs);CHKERRQ(ierr);
632 ierr = SetupDiscretization(dm, &ctx);CHKERRQ(ierr);
633 ierr = DMCreateGlobalVector(dm, &u);CHKERRQ(ierr);
634 ierr = PetscObjectSetName((PetscObject) u, "u");CHKERRQ(ierr);
635 ierr = VecDuplicate(u, &r);CHKERRQ(ierr);
636 ierr = PetscObjectSetName((PetscObject) r, "r");CHKERRQ(ierr);
637 /* create TS */
638 ierr = TSCreate(PETSC_COMM_WORLD, &ts);CHKERRQ(ierr);
639 ierr = TSSetDM(ts, dm);CHKERRQ(ierr);
640 ierr = TSSetApplicationContext(ts, &ctx);CHKERRQ(ierr);
641 ierr = DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &ctx);CHKERRQ(ierr);
642 if (ctx.implicit) {
643 ierr = DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &ctx);CHKERRQ(ierr);
644 ierr = DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &ctx);CHKERRQ(ierr);
645 } else {
646 ierr = DMTSSetRHSFunctionLocal(dm, DMPlexTSComputeRHSFunctionFVM, &ctx);CHKERRQ(ierr);
647 }
648 ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr);
649 ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
650 ierr = TSSetPostStep(ts, PostStep);CHKERRQ(ierr);
651 /* make solution & solve */
652 ierr = DMProjectFunction(dm, t, ctx.initialFuncs, (void **)ctxarr, INSERT_ALL_VALUES, u);CHKERRQ(ierr);
653 ierr = TSSetSolution(ts,u);CHKERRQ(ierr);
654 ierr = DMViewFromOptions(dm, NULL, "-dm_view");CHKERRQ(ierr);
655 ierr = PostStep(ts);CHKERRQ(ierr); /* print the initial state */
656 ierr = TSSolve(ts, u);CHKERRQ(ierr);
657 ierr = TSGetTime(ts, &t);CHKERRQ(ierr);
658 ierr = DMComputeL2Diff(dm, t, ctx.initialFuncs, (void **)ctxarr, u, &L2error);CHKERRQ(ierr);
659 if (L2error < 1.0e-11) {ierr = PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: < 1.0e-11\n");CHKERRQ(ierr);}
660 else {ierr = PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", L2error);CHKERRQ(ierr);}
661 #if 0
662 {
663 PetscReal res = 0.0;
664 /* Check discretization error */
665 ierr = VecViewFromOptions(u, NULL, "-initial_guess_view");CHKERRQ(ierr);
666 ierr = DMComputeL2Diff(dm, 0.0, ctx.exactFuncs, NULL, u, &error);CHKERRQ(ierr);
667 if (error < 1.0e-11) {ierr = PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: < 1.0e-11\n");CHKERRQ(ierr);}
668 else {ierr = PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", error);CHKERRQ(ierr);}
669 /* Check residual */
670 ierr = SNESComputeFunction(snes, u, r);CHKERRQ(ierr);
671 ierr = VecChop(r, 1.0e-10);CHKERRQ(ierr);
672 ierr = VecViewFromOptions(r, NULL, "-initial_residual_view");CHKERRQ(ierr);
673 ierr = VecNorm(r, NORM_2, &res);CHKERRQ(ierr);
674 ierr = PetscPrintf(PETSC_COMM_WORLD, "L_2 Residual: %g\n", res);CHKERRQ(ierr);
675 /* Check Jacobian */
676 {
677 Mat A;
678 Vec b;
679
680 ierr = SNESGetJacobian(snes, &A, NULL, NULL, NULL);CHKERRQ(ierr);
681 ierr = SNESComputeJacobian(snes, u, A, A);CHKERRQ(ierr);
682 ierr = VecDuplicate(u, &b);CHKERRQ(ierr);
683 ierr = VecSet(r, 0.0);CHKERRQ(ierr);
684 ierr = SNESComputeFunction(snes, r, b);CHKERRQ(ierr);
685 ierr = MatMult(A, u, r);CHKERRQ(ierr);
686 ierr = VecAXPY(r, 1.0, b);CHKERRQ(ierr);
687 ierr = VecDestroy(&b);CHKERRQ(ierr);
688 ierr = PetscPrintf(PETSC_COMM_WORLD, "Au - b = Au + F(0)\n");CHKERRQ(ierr);
689 ierr = VecChop(r, 1.0e-10);CHKERRQ(ierr);
690 ierr = VecViewFromOptions(r, NULL, "-linear_residual_view");CHKERRQ(ierr);
691 ierr = VecNorm(r, NORM_2, &res);CHKERRQ(ierr);
692 ierr = PetscPrintf(PETSC_COMM_WORLD, "Linear L_2 Residual: %g\n", res);CHKERRQ(ierr);
693 }
694 }
695 #endif
696 ierr = VecDestroy(&u);CHKERRQ(ierr);
697 ierr = VecDestroy(&r);CHKERRQ(ierr);
698 ierr = TSDestroy(&ts);CHKERRQ(ierr);
699 ierr = DMDestroy(&dm);CHKERRQ(ierr);
700 ierr = PetscFree(ctx.initialFuncs);CHKERRQ(ierr);
701 ierr = PetscFinalize();
702 return ierr;
703 }
704
705 /*TEST
706
707 test:
708 suffix: 0
709 args: -debug 1 -dim 2 -dm_refine 1 -x_periodicity PERIODIC -ts_max_steps 1 -ts_max_time 10. -ts_dt 1.0
710 test:
711 suffix: 1
712 args: -debug 1 -dim 3 -dm_refine 1 -z_periodicity PERIODIC -ts_max_steps 1 -ts_max_time 10. -ts_dt 1.0 -domain_lo -2,-1,-1 -domain_hi 2,1,1
713
714 TEST*/
715