1 static char help[] = "Vlasov-Poisson example of central orbits\n";
2
3 /*
4 To visualize the orbit, we can used
5
6 -ts_monitor_sp_swarm -ts_monitor_sp_swarm_retain -1 -ts_monitor_sp_swarm_phase 0 -draw_size 500,500
7
8 and we probably want it to run fast and not check convergence
9
10 -convest_num_refine 0 -ts_time_step 0.01 -ts_max_steps 100 -ts_max_time 100 -output_step 10
11 */
12
13 #include <petscts.h>
14 #include <petscdmplex.h>
15 #include <petscdmswarm.h>
16 #include <petsc/private/dmpleximpl.h> /* For norm and dot */
17 #include <petscfe.h>
18 #include <petscds.h>
19 #include <petsc/private/petscfeimpl.h> /* For interpolation */
20 #include <petscksp.h>
21 #include <petscsnes.h>
22
23 PETSC_EXTERN PetscErrorCode circleSingleX(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar[], void *);
24 PETSC_EXTERN PetscErrorCode circleSingleV(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar[], void *);
25 PETSC_EXTERN PetscErrorCode circleMultipleX(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar[], void *);
26 PETSC_EXTERN PetscErrorCode circleMultipleV(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar[], void *);
27
28 const char *EMTypes[] = {"primal", "mixed", "coulomb", "none", "EMType", "EM_", NULL};
29 typedef enum {
30 EM_PRIMAL,
31 EM_MIXED,
32 EM_COULOMB,
33 EM_NONE
34 } EMType;
35
36 typedef struct {
37 PetscBool error; /* Flag for printing the error */
38 PetscInt ostep; /* print the energy at each ostep time steps */
39 PetscReal timeScale; /* Nondimensionalizing time scale */
40 PetscReal sigma; /* Linear charge per box length */
41 EMType em; /* Type of electrostatic model */
42 SNES snes;
43 } AppCtx;
44
ProcessOptions(MPI_Comm comm,AppCtx * options)45 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
46 {
47 PetscFunctionBeginUser;
48 options->error = PETSC_FALSE;
49 options->ostep = 100;
50 options->timeScale = 1.0e-6;
51 options->sigma = 1.;
52 options->em = EM_COULOMB;
53
54 PetscOptionsBegin(comm, "", "Central Orbit Options", "DMSWARM");
55 PetscCall(PetscOptionsBool("-error", "Flag to print the error", "ex6.c", options->error, &options->error, NULL));
56 PetscCall(PetscOptionsInt("-output_step", "Number of time steps between output", "ex6.c", options->ostep, &options->ostep, NULL));
57 PetscCall(PetscOptionsReal("-sigma", "Linear charge per box length", "ex6.c", options->sigma, &options->sigma, NULL));
58 PetscCall(PetscOptionsReal("-timeScale", "Nondimensionalizing time scale", "ex6.c", options->timeScale, &options->timeScale, NULL));
59 PetscCall(PetscOptionsEnum("-em_type", "Type of electrostatic solver", "ex6.c", EMTypes, (PetscEnum)options->em, (PetscEnum *)&options->em, NULL));
60 PetscOptionsEnd();
61 PetscFunctionReturn(PETSC_SUCCESS);
62 }
63
CreateMesh(MPI_Comm comm,AppCtx * user,DM * dm)64 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
65 {
66 PetscFunctionBeginUser;
67 PetscCall(DMCreate(comm, dm));
68 PetscCall(DMSetType(*dm, DMPLEX));
69 PetscCall(DMSetFromOptions(*dm));
70 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
71 PetscFunctionReturn(PETSC_SUCCESS);
72 }
73
laplacian_f1(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[])74 static void laplacian_f1(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[])
75 {
76 PetscInt d;
77 for (d = 0; d < dim; ++d) f1[d] = u_x[d];
78 }
79
laplacian_g3(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[])80 static void laplacian_g3(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[])
81 {
82 PetscInt d;
83 for (d = 0; d < dim; ++d) g3[d * dim + d] = 1.0;
84 }
85
86 /*
87 / I grad\ /q\ = /0\
88 \-div 0 / \u/ \f/
89 */
f0_q(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[])90 static void f0_q(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[])
91 {
92 for (PetscInt c = 0; c < dim; ++c) f0[c] += u[uOff[0] + c];
93 }
94
f1_q(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[])95 static void f1_q(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[])
96 {
97 for (PetscInt d = 0; d < dim; ++d) f1[d * dim + d] += u[uOff[1]];
98 }
99
100 /* <t, q> */
g0_qq(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[])101 static void g0_qq(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[])
102 {
103 PetscInt c;
104 for (c = 0; c < dim; ++c) g0[c * dim + c] += 1.0;
105 }
106
g2_qu(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 g2[])107 static void g2_qu(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 g2[])
108 {
109 for (PetscInt d = 0; d < dim; ++d) g2[d * dim + d] += 1.0;
110 }
111
g1_uq(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[])112 static void g1_uq(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[])
113 {
114 for (PetscInt d = 0; d < dim; ++d) g1[d * dim + d] += 1.0;
115 }
116
CreateFEM(DM dm,AppCtx * user)117 static PetscErrorCode CreateFEM(DM dm, AppCtx *user)
118 {
119 PetscFE feu, feq;
120 PetscDS ds;
121 PetscBool simplex;
122 PetscInt dim;
123
124 PetscFunctionBeginUser;
125 PetscCall(DMGetDimension(dm, &dim));
126 PetscCall(DMPlexIsSimplex(dm, &simplex));
127 if (user->em == EM_MIXED) {
128 DMLabel label;
129
130 PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, dim, simplex, "field_", PETSC_DETERMINE, &feq));
131 PetscCall(PetscObjectSetName((PetscObject)feq, "field"));
132 PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, "potential_", PETSC_DETERMINE, &feu));
133 PetscCall(PetscObjectSetName((PetscObject)feu, "potential"));
134 PetscCall(PetscFECopyQuadrature(feq, feu));
135 PetscCall(DMSetField(dm, 0, NULL, (PetscObject)feq));
136 PetscCall(DMSetField(dm, 1, NULL, (PetscObject)feu));
137 PetscCall(DMCreateDS(dm));
138 PetscCall(PetscFEDestroy(&feu));
139 PetscCall(PetscFEDestroy(&feq));
140
141 PetscCall(DMGetLabel(dm, "marker", &label));
142 PetscCall(DMGetDS(dm, &ds));
143 PetscCall(PetscDSSetResidual(ds, 0, f0_q, f1_q));
144 PetscCall(PetscDSSetJacobian(ds, 0, 0, g0_qq, NULL, NULL, NULL));
145 PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_qu, NULL));
146 PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, g1_uq, NULL, NULL));
147 } else if (user->em == EM_PRIMAL) {
148 PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, NULL, PETSC_DETERMINE, &feu));
149 PetscCall(PetscObjectSetName((PetscObject)feu, "potential"));
150 PetscCall(DMSetField(dm, 0, NULL, (PetscObject)feu));
151 PetscCall(DMCreateDS(dm));
152 PetscCall(PetscFEDestroy(&feu));
153 PetscCall(DMGetDS(dm, &ds));
154 PetscCall(PetscDSSetResidual(ds, 0, NULL, laplacian_f1));
155 PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, laplacian_g3));
156 }
157 PetscFunctionReturn(PETSC_SUCCESS);
158 }
159
CreatePoisson(DM dm,AppCtx * user)160 static PetscErrorCode CreatePoisson(DM dm, AppCtx *user)
161 {
162 SNES snes;
163 Mat J;
164 MatNullSpace nullSpace;
165
166 PetscFunctionBeginUser;
167 PetscCall(CreateFEM(dm, user));
168 PetscCall(SNESCreate(PetscObjectComm((PetscObject)dm), &snes));
169 PetscCall(SNESSetOptionsPrefix(snes, "em_"));
170 PetscCall(SNESSetDM(snes, dm));
171 PetscCall(DMPlexSetSNESLocalFEM(dm, PETSC_FALSE, user));
172 PetscCall(SNESSetFromOptions(snes));
173
174 PetscCall(DMCreateMatrix(dm, &J));
175 PetscCall(MatNullSpaceCreate(PetscObjectComm((PetscObject)dm), PETSC_TRUE, 0, NULL, &nullSpace));
176 PetscCall(MatSetNullSpace(J, nullSpace));
177 PetscCall(MatNullSpaceDestroy(&nullSpace));
178 PetscCall(SNESSetJacobian(snes, J, J, NULL, NULL));
179 PetscCall(MatDestroy(&J));
180 user->snes = snes;
181 PetscFunctionReturn(PETSC_SUCCESS);
182 }
183
CreateSwarm(DM dm,AppCtx * user,DM * sw)184 static PetscErrorCode CreateSwarm(DM dm, AppCtx *user, DM *sw)
185 {
186 PetscReal v0[1] = {1.};
187 PetscInt dim;
188
189 PetscFunctionBeginUser;
190 PetscCall(DMGetDimension(dm, &dim));
191 PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), sw));
192 PetscCall(DMSetType(*sw, DMSWARM));
193 PetscCall(DMSetDimension(*sw, dim));
194 PetscCall(DMSwarmSetType(*sw, DMSWARM_PIC));
195 PetscCall(DMSwarmSetCellDM(*sw, dm));
196 PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "w_q", 1, PETSC_SCALAR));
197 PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "velocity", dim, PETSC_REAL));
198 PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "species", 1, PETSC_INT));
199 PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "initCoordinates", dim, PETSC_REAL));
200 PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "initVelocity", dim, PETSC_REAL));
201 PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "E_field", dim, PETSC_REAL));
202 PetscCall(DMSwarmFinalizeFieldRegister(*sw));
203 PetscCall(DMSwarmComputeLocalSizeFromOptions(*sw));
204 PetscCall(DMSwarmInitializeCoordinates(*sw));
205 PetscCall(DMSwarmInitializeVelocitiesFromOptions(*sw, v0));
206 PetscCall(DMSetFromOptions(*sw));
207 PetscCall(DMSetApplicationContext(*sw, user));
208 PetscCall(PetscObjectSetName((PetscObject)*sw, "Particles"));
209 PetscCall(DMViewFromOptions(*sw, NULL, "-sw_view"));
210 {
211 Vec gc, gc0, gv, gv0;
212
213 PetscCall(DMSwarmCreateGlobalVectorFromField(*sw, DMSwarmPICField_coor, &gc));
214 PetscCall(DMSwarmCreateGlobalVectorFromField(*sw, "initCoordinates", &gc0));
215 PetscCall(VecCopy(gc, gc0));
216 PetscCall(DMSwarmDestroyGlobalVectorFromField(*sw, DMSwarmPICField_coor, &gc));
217 PetscCall(DMSwarmDestroyGlobalVectorFromField(*sw, "initCoordinates", &gc0));
218 PetscCall(DMSwarmCreateGlobalVectorFromField(*sw, "velocity", &gv));
219 PetscCall(DMSwarmCreateGlobalVectorFromField(*sw, "initVelocity", &gv0));
220 PetscCall(VecCopy(gv, gv0));
221 PetscCall(DMSwarmDestroyGlobalVectorFromField(*sw, "velocity", &gv));
222 PetscCall(DMSwarmDestroyGlobalVectorFromField(*sw, "initVelocity", &gv0));
223 }
224 {
225 const char *fieldnames[2] = {DMSwarmPICField_coor, "velocity"};
226
227 PetscCall(DMSwarmVectorDefineFields(*sw, 2, fieldnames));
228 }
229 PetscFunctionReturn(PETSC_SUCCESS);
230 }
231
ComputeFieldAtParticles_Coulomb(SNES snes,DM sw,PetscReal E[])232 static PetscErrorCode ComputeFieldAtParticles_Coulomb(SNES snes, DM sw, PetscReal E[])
233 {
234 PetscReal *coords;
235 PetscInt dim, d, Np, p, q;
236 PetscMPIInt size;
237
238 PetscFunctionBegin;
239 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)snes), &size));
240 PetscCheck(size == 1, PetscObjectComm((PetscObject)snes), PETSC_ERR_SUP, "Coulomb code only works in serial");
241 PetscCall(DMGetDimension(sw, &dim));
242 PetscCall(DMSwarmGetLocalSize(sw, &Np));
243
244 PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
245 for (p = 0; p < Np; ++p) {
246 PetscReal *pcoord = &coords[p * dim];
247 PetscReal *pE = &E[p * dim];
248 /* Calculate field at particle p due to particle q */
249 for (q = 0; q < Np; ++q) {
250 PetscReal *qcoord = &coords[q * dim];
251 PetscReal rpq[3], r;
252
253 if (p == q) continue;
254 for (d = 0; d < dim; ++d) rpq[d] = pcoord[d] - qcoord[d];
255 r = DMPlex_NormD_Internal(dim, rpq);
256 for (d = 0; d < dim; ++d) pE[d] += rpq[d] / PetscPowRealInt(r, 3);
257 }
258 }
259 PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
260 PetscFunctionReturn(PETSC_SUCCESS);
261 }
262
ComputeFieldAtParticles_Primal(SNES snes,DM sw,PetscReal E[])263 static PetscErrorCode ComputeFieldAtParticles_Primal(SNES snes, DM sw, PetscReal E[])
264 {
265 DM dm;
266 PetscDS ds;
267 PetscFE fe;
268 Mat M_p;
269 Vec phi, locPhi, rho, f;
270 PetscReal *coords, chargeTol = 1e-13;
271 PetscInt dim, d, cStart, cEnd, c, Np;
272 const char **oldFields;
273 PetscInt Nf;
274 const char **tmp;
275
276 PetscFunctionBegin;
277 PetscCall(DMGetDimension(sw, &dim));
278 PetscCall(DMSwarmGetLocalSize(sw, &Np));
279 PetscCall(SNESGetDM(snes, &dm));
280
281 PetscCall(DMSwarmVectorGetField(sw, &Nf, &tmp));
282 PetscCall(PetscMalloc1(Nf, &oldFields));
283 for (PetscInt f = 0; f < Nf; ++f) PetscCall(PetscStrallocpy(tmp[f], (char **)&oldFields[f]));
284 PetscCall(DMSwarmVectorDefineField(sw, "w_q"));
285 PetscCall(DMCreateMassMatrix(sw, dm, &M_p));
286 PetscCall(DMSwarmVectorDefineFields(sw, Nf, oldFields));
287 for (PetscInt f = 0; f < Nf; ++f) PetscCall(PetscFree(oldFields[f]));
288 PetscCall(PetscFree(oldFields));
289
290 /* Create the charges rho */
291 PetscCall(DMGetGlobalVector(dm, &rho));
292 PetscCall(PetscObjectSetName((PetscObject)rho, "rho"));
293 PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "w_q", &f));
294 PetscCall(PetscObjectSetName((PetscObject)f, "particle weight"));
295 PetscCall(MatMultTranspose(M_p, f, rho));
296 PetscCall(MatViewFromOptions(M_p, NULL, "-mp_view"));
297 PetscCall(VecViewFromOptions(f, NULL, "-weights_view"));
298 PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f));
299 PetscCall(MatDestroy(&M_p));
300 {
301 PetscScalar sum;
302 PetscInt n;
303 PetscReal phi_0 = 1.; /* (sigma*sigma*sigma)*(timeScale*timeScale)/(m_e*q_e*epsi_0)*/
304
305 /* Remove constant from rho */
306 PetscCall(VecGetSize(rho, &n));
307 PetscCall(VecSum(rho, &sum));
308 PetscCall(VecShift(rho, -sum / n));
309 PetscCall(VecSum(rho, &sum));
310 PetscCheck(PetscAbsScalar(sum) < chargeTol, PetscObjectComm((PetscObject)sw), PETSC_ERR_PLIB, "Charge should have no DC component: %g", (double)PetscRealPart(sum));
311 /* Nondimensionalize rho */
312 PetscCall(VecScale(rho, phi_0));
313 }
314 PetscCall(VecViewFromOptions(rho, NULL, "-poisson_rho_view"));
315
316 PetscCall(DMGetGlobalVector(dm, &phi));
317 PetscCall(PetscObjectSetName((PetscObject)phi, "potential"));
318 PetscCall(VecSet(phi, 0.0));
319 PetscCall(SNESSolve(snes, rho, phi));
320 PetscCall(DMRestoreGlobalVector(dm, &rho));
321 PetscCall(VecViewFromOptions(phi, NULL, "-phi_view"));
322
323 PetscCall(DMGetLocalVector(dm, &locPhi));
324 PetscCall(DMGlobalToLocalBegin(dm, phi, INSERT_VALUES, locPhi));
325 PetscCall(DMGlobalToLocalEnd(dm, phi, INSERT_VALUES, locPhi));
326 PetscCall(DMRestoreGlobalVector(dm, &phi));
327
328 PetscCall(DMGetDS(dm, &ds));
329 PetscCall(PetscDSGetDiscretization(ds, 0, (PetscObject *)&fe));
330 PetscCall(DMSwarmSortGetAccess(sw));
331 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
332 PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
333 for (c = cStart; c < cEnd; ++c) {
334 PetscTabulation tab;
335 PetscScalar *clPhi = NULL;
336 PetscReal *pcoord, *refcoord;
337 PetscReal v[3], J[9], invJ[9], detJ;
338 PetscInt *points;
339 PetscInt Ncp, cp;
340
341 PetscCall(DMSwarmSortGetPointsPerCell(sw, c, &Ncp, &points));
342 PetscCall(DMGetWorkArray(dm, Ncp * dim, MPIU_REAL, &pcoord));
343 PetscCall(DMGetWorkArray(dm, Ncp * dim, MPIU_REAL, &refcoord));
344 for (cp = 0; cp < Ncp; ++cp)
345 for (d = 0; d < dim; ++d) pcoord[cp * dim + d] = coords[points[cp] * dim + d];
346 PetscCall(DMPlexCoordinatesToReference(dm, c, Ncp, pcoord, refcoord));
347 PetscCall(PetscFECreateTabulation(fe, 1, Ncp, refcoord, 1, &tab));
348 PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v, J, invJ, &detJ));
349 PetscCall(DMPlexVecGetClosure(dm, NULL, locPhi, c, NULL, &clPhi));
350 for (cp = 0; cp < Ncp; ++cp) {
351 const PetscReal *basisDer = tab->T[1];
352 const PetscInt p = points[cp];
353
354 for (d = 0; d < dim; ++d) E[p * dim + d] = 0.;
355 PetscCall(PetscFEFreeInterpolateGradient_Static(fe, basisDer, clPhi, dim, invJ, NULL, cp, &E[p * dim]));
356 for (d = 0; d < dim; ++d) E[p * dim + d] *= -1.0;
357 }
358 PetscCall(DMPlexVecRestoreClosure(dm, NULL, locPhi, c, NULL, &clPhi));
359 PetscCall(DMRestoreWorkArray(dm, Ncp * dim, MPIU_REAL, &pcoord));
360 PetscCall(DMRestoreWorkArray(dm, Ncp * dim, MPIU_REAL, &refcoord));
361 PetscCall(PetscTabulationDestroy(&tab));
362 PetscCall(DMSwarmSortRestorePointsPerCell(sw, c, &Ncp, &points));
363 }
364 PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
365 PetscCall(DMSwarmSortRestoreAccess(sw));
366 PetscCall(DMRestoreLocalVector(dm, &locPhi));
367 PetscFunctionReturn(PETSC_SUCCESS);
368 }
369
ComputeFieldAtParticles_Mixed(SNES snes,DM sw,PetscReal E[])370 static PetscErrorCode ComputeFieldAtParticles_Mixed(SNES snes, DM sw, PetscReal E[])
371 {
372 DM dm, potential_dm;
373 IS potential_IS;
374 PetscDS ds;
375 PetscFE fe;
376 PetscFEGeom feGeometry;
377 Mat M_p;
378 Vec phi, locPhi, rho, f, temp_rho;
379 PetscQuadrature q;
380 PetscReal *coords, chargeTol = 1e-13;
381 PetscInt dim, d, cStart, cEnd, c, Np, pot_field = 1;
382 const char **oldFields;
383 PetscInt Nf;
384 const char **tmp;
385
386 PetscFunctionBegin;
387 PetscCall(DMGetDimension(sw, &dim));
388 PetscCall(DMSwarmGetLocalSize(sw, &Np));
389
390 /* Create the charges rho */
391 PetscCall(SNESGetDM(snes, &dm));
392 PetscCall(DMGetGlobalVector(dm, &rho));
393 PetscCall(PetscObjectSetName((PetscObject)rho, "rho"));
394 PetscCall(DMCreateSubDM(dm, 1, &pot_field, &potential_IS, &potential_dm));
395
396 PetscCall(DMSwarmVectorGetField(sw, &Nf, &tmp));
397 PetscCall(PetscMalloc1(Nf, &oldFields));
398 for (PetscInt f = 0; f < Nf; ++f) PetscCall(PetscStrallocpy(tmp[f], (char **)&oldFields[f]));
399 PetscCall(DMSwarmVectorDefineField(sw, "w_q"));
400 PetscCall(DMCreateMassMatrix(sw, potential_dm, &M_p));
401 PetscCall(DMSwarmVectorDefineFields(sw, Nf, oldFields));
402 for (PetscInt f = 0; f < Nf; ++f) PetscCall(PetscFree(oldFields[f]));
403 PetscCall(PetscFree(oldFields));
404
405 PetscCall(MatViewFromOptions(M_p, NULL, "-mp_view"));
406 PetscCall(DMGetGlobalVector(potential_dm, &temp_rho));
407 PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "w_q", &f));
408 PetscCall(PetscObjectSetName((PetscObject)f, "particle weight"));
409 PetscCall(VecViewFromOptions(f, NULL, "-weights_view"));
410 PetscCall(MatMultTranspose(M_p, f, temp_rho));
411 PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f));
412 PetscCall(MatDestroy(&M_p));
413 PetscCall(PetscObjectSetName((PetscObject)rho, "rho"));
414 PetscCall(VecViewFromOptions(rho, NULL, "-poisson_rho_view"));
415 PetscCall(VecISCopy(rho, potential_IS, SCATTER_FORWARD, temp_rho));
416 PetscCall(DMRestoreGlobalVector(potential_dm, &temp_rho));
417 PetscCall(DMDestroy(&potential_dm));
418 PetscCall(ISDestroy(&potential_IS));
419 {
420 PetscScalar sum;
421 PetscInt n;
422 PetscReal phi_0 = 1.; /*(sigma*sigma*sigma)*(timeScale*timeScale)/(m_e*q_e*epsi_0);*/
423
424 /*Remove constant from rho*/
425 PetscCall(VecGetSize(rho, &n));
426 PetscCall(VecSum(rho, &sum));
427 PetscCall(VecShift(rho, -sum / n));
428 PetscCall(VecSum(rho, &sum));
429 PetscCheck(PetscAbsScalar(sum) < chargeTol, PetscObjectComm((PetscObject)sw), PETSC_ERR_PLIB, "Charge should have no DC component: %g", (double)PetscRealPart(sum));
430 /* Nondimensionalize rho */
431 PetscCall(VecScale(rho, phi_0));
432 }
433 PetscCall(DMGetGlobalVector(dm, &phi));
434 PetscCall(PetscObjectSetName((PetscObject)phi, "potential"));
435 PetscCall(VecSet(phi, 0.0));
436 PetscCall(SNESSolve(snes, NULL, phi));
437 PetscCall(DMRestoreGlobalVector(dm, &rho));
438 PetscCall(VecViewFromOptions(phi, NULL, "-phi_view"));
439
440 PetscCall(DMGetLocalVector(dm, &locPhi));
441 PetscCall(DMGlobalToLocalBegin(dm, phi, INSERT_VALUES, locPhi));
442 PetscCall(DMGlobalToLocalEnd(dm, phi, INSERT_VALUES, locPhi));
443 PetscCall(DMRestoreGlobalVector(dm, &phi));
444
445 PetscCall(DMGetDS(dm, &ds));
446 PetscCall(PetscDSGetDiscretization(ds, 0, (PetscObject *)&fe));
447 PetscCall(DMSwarmSortGetAccess(sw));
448 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
449 PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
450 for (c = cStart; c < cEnd; ++c) {
451 PetscTabulation tab;
452 PetscScalar *clPhi = NULL;
453 PetscReal *pcoord, *refcoord;
454 PetscReal v[3], J[9], invJ[9], detJ;
455 PetscInt *points;
456 PetscInt Ncp, cp;
457
458 PetscCall(DMSwarmSortGetPointsPerCell(sw, c, &Ncp, &points));
459 PetscCall(DMGetWorkArray(dm, Ncp * dim, MPIU_REAL, &pcoord));
460 PetscCall(DMGetWorkArray(dm, Ncp * dim, MPIU_REAL, &refcoord));
461 for (cp = 0; cp < Ncp; ++cp)
462 for (d = 0; d < dim; ++d) pcoord[cp * dim + d] = coords[points[cp] * dim + d];
463 PetscCall(DMPlexCoordinatesToReference(dm, c, Ncp, pcoord, refcoord));
464 PetscCall(PetscFECreateTabulation(fe, 1, Ncp, refcoord, 1, &tab));
465 PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v, J, invJ, &detJ));
466 PetscCall(DMPlexVecGetClosure(dm, NULL, locPhi, c, NULL, &clPhi));
467 for (cp = 0; cp < Ncp; ++cp) {
468 const PetscInt p = points[cp];
469
470 for (d = 0; d < dim; ++d) E[p * dim + d] = 0.;
471 PetscCall(PetscFEGetQuadrature(fe, &q));
472 PetscCall(PetscFECreateCellGeometry(fe, q, &feGeometry));
473 PetscCall(PetscFEInterpolateAtPoints_Static(fe, tab, clPhi, &feGeometry, cp, &E[p * dim]));
474 PetscCall(PetscFEDestroyCellGeometry(fe, &feGeometry));
475 }
476 PetscCall(DMPlexVecRestoreClosure(dm, NULL, locPhi, c, NULL, &clPhi));
477 PetscCall(DMRestoreWorkArray(dm, Ncp * dim, MPIU_REAL, &pcoord));
478 PetscCall(DMRestoreWorkArray(dm, Ncp * dim, MPIU_REAL, &refcoord));
479 PetscCall(PetscTabulationDestroy(&tab));
480 PetscCall(DMSwarmSortRestorePointsPerCell(sw, c, &Ncp, &points));
481 }
482 PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
483 PetscCall(DMSwarmSortRestoreAccess(sw));
484 PetscCall(DMRestoreLocalVector(dm, &locPhi));
485 PetscFunctionReturn(PETSC_SUCCESS);
486 }
487
ComputeFieldAtParticles(SNES snes,DM sw,PetscReal E[])488 static PetscErrorCode ComputeFieldAtParticles(SNES snes, DM sw, PetscReal E[])
489 {
490 AppCtx *ctx;
491 PetscInt dim, Np;
492
493 PetscFunctionBegin;
494 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
495 PetscValidHeaderSpecific(sw, DM_CLASSID, 2);
496 PetscAssertPointer(E, 3);
497 PetscCall(DMGetDimension(sw, &dim));
498 PetscCall(DMSwarmGetLocalSize(sw, &Np));
499 PetscCall(DMGetApplicationContext(sw, &ctx));
500 PetscCall(PetscArrayzero(E, Np * dim));
501
502 switch (ctx->em) {
503 case EM_PRIMAL:
504 PetscCall(ComputeFieldAtParticles_Primal(snes, sw, E));
505 break;
506 case EM_COULOMB:
507 PetscCall(ComputeFieldAtParticles_Coulomb(snes, sw, E));
508 break;
509 case EM_MIXED:
510 PetscCall(ComputeFieldAtParticles_Mixed(snes, sw, E));
511 break;
512 case EM_NONE:
513 break;
514 default:
515 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No solver for electrostatic model %s", EMTypes[ctx->em]);
516 }
517 PetscFunctionReturn(PETSC_SUCCESS);
518 }
519
RHSFunction(TS ts,PetscReal t,Vec U,Vec G,PetscCtx ctx)520 static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec G, PetscCtx ctx)
521 {
522 DM sw;
523 SNES snes = ((AppCtx *)ctx)->snes;
524 const PetscReal *coords, *vel;
525 const PetscScalar *u;
526 PetscScalar *g;
527 PetscReal *E;
528 PetscInt dim, d, Np, p;
529
530 PetscFunctionBeginUser;
531 PetscCall(TSGetDM(ts, &sw));
532 PetscCall(DMGetDimension(sw, &dim));
533 PetscCall(DMSwarmGetField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
534 PetscCall(DMSwarmGetField(sw, "initVelocity", NULL, NULL, (void **)&vel));
535 PetscCall(DMSwarmGetField(sw, "E_field", NULL, NULL, (void **)&E));
536 PetscCall(VecGetLocalSize(U, &Np));
537 PetscCall(VecGetArrayRead(U, &u));
538 PetscCall(VecGetArray(G, &g));
539
540 PetscCall(ComputeFieldAtParticles(snes, sw, E));
541
542 Np /= 2 * dim;
543 for (p = 0; p < Np; ++p) {
544 const PetscReal x0 = coords[p * dim + 0];
545 const PetscReal vy0 = vel[p * dim + 1];
546 const PetscReal omega = vy0 / x0;
547
548 for (d = 0; d < dim; ++d) {
549 g[(p * 2 + 0) * dim + d] = u[(p * 2 + 1) * dim + d];
550 g[(p * 2 + 1) * dim + d] = E[p * dim + d] - PetscSqr(omega) * u[(p * 2 + 0) * dim + d];
551 }
552 }
553 PetscCall(DMSwarmRestoreField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
554 PetscCall(DMSwarmRestoreField(sw, "initVelocity", NULL, NULL, (void **)&vel));
555 PetscCall(DMSwarmRestoreField(sw, "E_field", NULL, NULL, (void **)&E));
556 PetscCall(VecRestoreArrayRead(U, &u));
557 PetscCall(VecRestoreArray(G, &g));
558 PetscFunctionReturn(PETSC_SUCCESS);
559 }
560
561 /* J_{ij} = dF_i/dx_j
562 J_p = ( 0 1)
563 (-w^2 0)
564 TODO Now there is another term with w^2 from the electric field. I think we will need to invert the operator.
565 Perhaps we can approximate the Jacobian using only the cellwise P-P gradient from Coulomb
566 */
RHSJacobian(TS ts,PetscReal t,Vec U,Mat J,Mat P,PetscCtx ctx)567 static PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat J, Mat P, PetscCtx ctx)
568 {
569 DM sw;
570 const PetscReal *coords, *vel;
571 PetscInt dim, d, Np, p, rStart;
572
573 PetscFunctionBeginUser;
574 PetscCall(TSGetDM(ts, &sw));
575 PetscCall(DMGetDimension(sw, &dim));
576 PetscCall(VecGetLocalSize(U, &Np));
577 PetscCall(MatGetOwnershipRange(J, &rStart, NULL));
578 PetscCall(DMSwarmGetField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
579 PetscCall(DMSwarmGetField(sw, "initVelocity", NULL, NULL, (void **)&vel));
580 Np /= 2 * dim;
581 for (p = 0; p < Np; ++p) {
582 const PetscReal x0 = coords[p * dim + 0];
583 const PetscReal vy0 = vel[p * dim + 1];
584 const PetscReal omega = vy0 / x0;
585 PetscScalar vals[4] = {0., 1., -PetscSqr(omega), 0.};
586
587 for (d = 0; d < dim; ++d) {
588 const PetscInt rows[2] = {(p * 2 + 0) * dim + d + rStart, (p * 2 + 1) * dim + d + rStart};
589 PetscCall(MatSetValues(J, 2, rows, 2, rows, vals, INSERT_VALUES));
590 }
591 }
592 PetscCall(DMSwarmRestoreField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
593 PetscCall(DMSwarmRestoreField(sw, "initVelocity", NULL, NULL, (void **)&vel));
594 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
595 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
596 PetscFunctionReturn(PETSC_SUCCESS);
597 }
598
RHSFunctionX(TS ts,PetscReal t,Vec V,Vec Xres,PetscCtx ctx)599 static PetscErrorCode RHSFunctionX(TS ts, PetscReal t, Vec V, Vec Xres, PetscCtx ctx)
600 {
601 DM sw;
602 const PetscScalar *v;
603 PetscScalar *xres;
604 PetscInt Np, p, dim, d;
605
606 PetscFunctionBeginUser;
607 PetscCall(TSGetDM(ts, &sw));
608 PetscCall(DMGetDimension(sw, &dim));
609 PetscCall(VecGetLocalSize(Xres, &Np));
610 Np /= dim;
611 PetscCall(VecGetArrayRead(V, &v));
612 PetscCall(VecGetArray(Xres, &xres));
613 for (p = 0; p < Np; ++p) {
614 for (d = 0; d < dim; ++d) xres[p * dim + d] = v[p * dim + d];
615 }
616 PetscCall(VecRestoreArrayRead(V, &v));
617 PetscCall(VecRestoreArray(Xres, &xres));
618 PetscFunctionReturn(PETSC_SUCCESS);
619 }
620
RHSFunctionV(TS ts,PetscReal t,Vec X,Vec Vres,PetscCtx ctx)621 static PetscErrorCode RHSFunctionV(TS ts, PetscReal t, Vec X, Vec Vres, PetscCtx ctx)
622 {
623 DM sw;
624 SNES snes = ((AppCtx *)ctx)->snes;
625 const PetscScalar *x;
626 const PetscReal *coords, *vel;
627 PetscScalar *vres;
628 PetscReal *E;
629 PetscInt Np, p, dim, d;
630
631 PetscFunctionBeginUser;
632 PetscCall(TSGetDM(ts, &sw));
633 PetscCall(DMGetDimension(sw, &dim));
634 PetscCall(DMSwarmGetField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
635 PetscCall(DMSwarmGetField(sw, "initVelocity", NULL, NULL, (void **)&vel));
636 PetscCall(DMSwarmGetField(sw, "E_field", NULL, NULL, (void **)&E));
637 PetscCall(VecGetLocalSize(Vres, &Np));
638 PetscCall(VecGetArrayRead(X, &x));
639 PetscCall(VecGetArray(Vres, &vres));
640 PetscCheck(dim == 2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Dimension must be 2");
641
642 PetscCall(ComputeFieldAtParticles(snes, sw, E));
643
644 Np /= dim;
645 for (p = 0; p < Np; ++p) {
646 const PetscReal x0 = coords[p * dim + 0];
647 const PetscReal vy0 = vel[p * dim + 1];
648 const PetscReal omega = vy0 / x0;
649
650 for (d = 0; d < dim; ++d) vres[p * dim + d] = E[p * dim + d] - PetscSqr(omega) * x[p * dim + d];
651 }
652 PetscCall(VecRestoreArrayRead(X, &x));
653 PetscCall(VecRestoreArray(Vres, &vres));
654 PetscCall(DMSwarmRestoreField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
655 PetscCall(DMSwarmRestoreField(sw, "initVelocity", NULL, NULL, (void **)&vel));
656 PetscCall(DMSwarmRestoreField(sw, "E_field", NULL, NULL, (void **)&E));
657 PetscFunctionReturn(PETSC_SUCCESS);
658 }
659
660 /* Discrete Gradients Formulation: S, F, gradF (G) */
RHSJacobianS(TS ts,PetscReal t,Vec U,Mat S,PetscCtx ctx)661 PetscErrorCode RHSJacobianS(TS ts, PetscReal t, Vec U, Mat S, PetscCtx ctx)
662 {
663 PetscScalar vals[4] = {0., 1., -1., 0.};
664 DM sw;
665 PetscInt dim, d, Np, p, rStart;
666
667 PetscFunctionBeginUser;
668 PetscCall(TSGetDM(ts, &sw));
669 PetscCall(DMGetDimension(sw, &dim));
670 PetscCall(VecGetLocalSize(U, &Np));
671 PetscCall(MatGetOwnershipRange(S, &rStart, NULL));
672 Np /= 2 * dim;
673 for (p = 0; p < Np; ++p) {
674 for (d = 0; d < dim; ++d) {
675 const PetscInt rows[2] = {(p * 2 + 0) * dim + d + rStart, (p * 2 + 1) * dim + d + rStart};
676 PetscCall(MatSetValues(S, 2, rows, 2, rows, vals, INSERT_VALUES));
677 }
678 }
679 PetscCall(MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY));
680 PetscCall(MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY));
681 PetscFunctionReturn(PETSC_SUCCESS);
682 }
683
RHSObjectiveF(TS ts,PetscReal t,Vec U,PetscScalar * F,PetscCtx ctx)684 PetscErrorCode RHSObjectiveF(TS ts, PetscReal t, Vec U, PetscScalar *F, PetscCtx ctx)
685 {
686 SNES snes = ((AppCtx *)ctx)->snes;
687 DM dm, sw;
688 const PetscScalar *u, *phi_vals;
689 PetscInt dim, Np, cStart, cEnd;
690 PetscReal *vel, *coords, m_p = 1., q_p = -1.;
691 Vec phi;
692
693 PetscFunctionBeginUser;
694 PetscCall(TSGetDM(ts, &sw));
695 PetscCall(DMGetDimension(sw, &dim));
696 PetscCall(SNESGetDM(snes, &dm));
697 PetscCall(VecGetArrayRead(U, &u));
698 PetscCall(VecGetLocalSize(U, &Np));
699 PetscCall(DMGetGlobalVector(dm, &phi));
700 PetscCall(VecViewFromOptions(phi, NULL, "-phi_view_dg"));
701 PetscCall(PetscObjectSetName((PetscObject)phi, "potential"));
702 PetscInt phi_size;
703 PetscCall(VecGetSize(phi, &phi_size));
704 PetscCall(VecGetArrayRead(phi, &phi_vals));
705 PetscCall(DMSwarmGetField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
706 PetscCall(DMSwarmGetField(sw, "initVelocity", NULL, NULL, (void **)&vel));
707
708 PetscCall(DMSwarmSortGetAccess(sw));
709 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
710 Np /= 2 * dim;
711 for (PetscInt c = cStart; c < cEnd; ++c) {
712 PetscInt *points;
713 PetscInt Ncp;
714 PetscReal E = phi_vals[c];
715
716 PetscCall(DMSwarmSortGetPointsPerCell(sw, c, &Ncp, &points));
717 for (PetscInt cp = 0; cp < Ncp; ++cp) {
718 const PetscInt p = points[cp];
719 const PetscReal x0 = coords[p * dim + 0];
720 const PetscReal vy0 = vel[p * dim + 1];
721 const PetscReal omega = vy0 / x0;
722 const PetscReal v2 = DMPlex_DotRealD_Internal(dim, &u[(p * 2 + 1) * dim], &u[(p * 2 + 1) * dim]);
723 const PetscReal x2 = DMPlex_DotRealD_Internal(dim, &u[(p * 2 + 0) * dim], &u[(p * 2 + 0) * dim]);
724 E += 0.5 * q_p * m_p * (v2) + 0.5 * PetscSqr(omega) * (x2);
725
726 *F += E;
727 }
728 PetscCall(DMSwarmSortRestorePointsPerCell(sw, c, &Ncp, &points));
729 }
730 PetscCall(DMSwarmSortRestoreAccess(sw));
731 PetscCall(DMSwarmRestoreField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
732 PetscCall(DMSwarmRestoreField(sw, "initVelocity", NULL, NULL, (void **)&vel));
733 PetscCall(VecRestoreArrayRead(phi, &phi_vals));
734 PetscCall(DMRestoreGlobalVector(dm, &phi));
735 // PetscCall(DMSwarmRestoreField(sw, "potential", NULL, NULL, (void **)&pot));
736 PetscCall(VecRestoreArrayRead(U, &u));
737 PetscFunctionReturn(PETSC_SUCCESS);
738 }
739
740 /* dF/dx = q E dF/dv = v */
RHSFunctionG(TS ts,PetscReal t,Vec U,Vec G,PetscCtx ctx)741 PetscErrorCode RHSFunctionG(TS ts, PetscReal t, Vec U, Vec G, PetscCtx ctx)
742 {
743 DM sw;
744 SNES snes = ((AppCtx *)ctx)->snes;
745 const PetscReal *coords, *vel;
746 const PetscScalar *u;
747 PetscScalar *g;
748 PetscReal *E, m_p = 1., q_p = -1.;
749 PetscInt dim, d, Np, p;
750
751 PetscFunctionBeginUser;
752 PetscCall(TSGetDM(ts, &sw));
753 PetscCall(DMGetDimension(sw, &dim));
754 PetscCall(DMSwarmGetField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
755 PetscCall(DMSwarmGetField(sw, "initVelocity", NULL, NULL, (void **)&vel));
756 PetscCall(DMSwarmGetField(sw, "E_field", NULL, NULL, (void **)&E));
757 PetscCall(DMSwarmGetLocalSize(sw, &Np));
758 PetscCall(VecGetArrayRead(U, &u));
759 PetscCall(VecGetArray(G, &g));
760
761 int COMPUTEFIELD;
762 PetscCall(PetscLogEventRegister("COMPFIELDATPART", TS_CLASSID, &COMPUTEFIELD));
763 PetscCall(PetscLogEventBegin(COMPUTEFIELD, 0, 0, 0, 0));
764 PetscCall(ComputeFieldAtParticles(snes, sw, E));
765 PetscCall(PetscLogEventEnd(COMPUTEFIELD, 0, 0, 0, 0));
766 for (p = 0; p < Np; ++p) {
767 const PetscReal x0 = coords[p * dim + 0];
768 const PetscReal vy0 = vel[p * dim + 1];
769 const PetscReal omega = vy0 / x0;
770 for (d = 0; d < dim; ++d) {
771 g[(p * 2 + 0) * dim + d] = -(q_p / m_p) * E[p * dim + d] + PetscSqr(omega) * u[(p * 2 + 0) * dim + d];
772 g[(p * 2 + 1) * dim + d] = u[(p * 2 + 1) * dim + d];
773 }
774 }
775 PetscCall(DMSwarmRestoreField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
776 PetscCall(DMSwarmRestoreField(sw, "initVelocity", NULL, NULL, (void **)&vel));
777 PetscCall(DMSwarmRestoreField(sw, "E_field", NULL, NULL, (void **)&E));
778 PetscCall(VecRestoreArrayRead(U, &u));
779 PetscCall(VecRestoreArray(G, &g));
780 PetscFunctionReturn(PETSC_SUCCESS);
781 }
782
CreateSolution(TS ts)783 static PetscErrorCode CreateSolution(TS ts)
784 {
785 DM sw;
786 Vec u;
787 PetscInt dim, Np;
788
789 PetscFunctionBegin;
790 PetscCall(TSGetDM(ts, &sw));
791 PetscCall(DMGetDimension(sw, &dim));
792 PetscCall(DMSwarmGetLocalSize(sw, &Np));
793 PetscCall(VecCreate(PETSC_COMM_WORLD, &u));
794 PetscCall(VecSetBlockSize(u, dim));
795 PetscCall(VecSetSizes(u, 2 * Np * dim, PETSC_DECIDE));
796 PetscCall(VecSetUp(u));
797 PetscCall(TSSetSolution(ts, u));
798 PetscCall(VecDestroy(&u));
799 PetscFunctionReturn(PETSC_SUCCESS);
800 }
801
SetProblem(TS ts)802 static PetscErrorCode SetProblem(TS ts)
803 {
804 AppCtx *user;
805 DM sw;
806
807 PetscFunctionBegin;
808 PetscCall(TSGetDM(ts, &sw));
809 PetscCall(DMGetApplicationContext(sw, &user));
810 // Define unified system for (X, V)
811 {
812 Mat J;
813 PetscInt dim, Np;
814
815 PetscCall(DMGetDimension(sw, &dim));
816 PetscCall(DMSwarmGetLocalSize(sw, &Np));
817 PetscCall(MatCreate(PETSC_COMM_WORLD, &J));
818 PetscCall(MatSetSizes(J, 2 * Np * dim, 2 * Np * dim, PETSC_DECIDE, PETSC_DECIDE));
819 PetscCall(MatSetBlockSize(J, 2 * dim));
820 PetscCall(MatSetFromOptions(J));
821 PetscCall(MatSetUp(J));
822 PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, user));
823 PetscCall(TSSetRHSJacobian(ts, J, J, RHSJacobian, user));
824 PetscCall(MatDestroy(&J));
825 }
826 /* Define split system for X and V */
827 {
828 Vec u;
829 IS isx, isv, istmp;
830 const PetscInt *idx;
831 PetscInt dim, Np, rstart;
832
833 PetscCall(TSGetSolution(ts, &u));
834 PetscCall(DMGetDimension(sw, &dim));
835 PetscCall(DMSwarmGetLocalSize(sw, &Np));
836 PetscCall(VecGetOwnershipRange(u, &rstart, NULL));
837 PetscCall(ISCreateStride(PETSC_COMM_WORLD, Np, (rstart / dim) + 0, 2, &istmp));
838 PetscCall(ISGetIndices(istmp, &idx));
839 PetscCall(ISCreateBlock(PETSC_COMM_WORLD, dim, Np, idx, PETSC_COPY_VALUES, &isx));
840 PetscCall(ISRestoreIndices(istmp, &idx));
841 PetscCall(ISDestroy(&istmp));
842 PetscCall(ISCreateStride(PETSC_COMM_WORLD, Np, (rstart / dim) + 1, 2, &istmp));
843 PetscCall(ISGetIndices(istmp, &idx));
844 PetscCall(ISCreateBlock(PETSC_COMM_WORLD, dim, Np, idx, PETSC_COPY_VALUES, &isv));
845 PetscCall(ISRestoreIndices(istmp, &idx));
846 PetscCall(ISDestroy(&istmp));
847 PetscCall(TSRHSSplitSetIS(ts, "position", isx));
848 PetscCall(TSRHSSplitSetIS(ts, "momentum", isv));
849 PetscCall(ISDestroy(&isx));
850 PetscCall(ISDestroy(&isv));
851 PetscCall(TSRHSSplitSetRHSFunction(ts, "position", NULL, RHSFunctionX, user));
852 PetscCall(TSRHSSplitSetRHSFunction(ts, "momentum", NULL, RHSFunctionV, user));
853 }
854 // Define symplectic formulation U_t = S . G, where G = grad F
855 {
856 PetscCall(TSDiscGradSetFormulation(ts, RHSJacobianS, RHSObjectiveF, RHSFunctionG, user));
857 }
858 PetscFunctionReturn(PETSC_SUCCESS);
859 }
860
circleSingleX(PetscInt dim,PetscReal time,const PetscReal unused[],PetscInt p,PetscScalar x[],PetscCtx ctx)861 PetscErrorCode circleSingleX(PetscInt dim, PetscReal time, const PetscReal unused[], PetscInt p, PetscScalar x[], PetscCtx ctx)
862 {
863 x[0] = p + 1.;
864 x[1] = 0.;
865 return PETSC_SUCCESS;
866 }
867
circleSingleV(PetscInt dim,PetscReal time,const PetscReal unused[],PetscInt p,PetscScalar v[],PetscCtx ctx)868 PetscErrorCode circleSingleV(PetscInt dim, PetscReal time, const PetscReal unused[], PetscInt p, PetscScalar v[], PetscCtx ctx)
869 {
870 v[0] = 0.;
871 v[1] = PetscSqrtReal(1000. / (p + 1.));
872 return PETSC_SUCCESS;
873 }
874
875 /* Put 5 particles into each circle */
circleMultipleX(PetscInt dim,PetscReal time,const PetscReal unused[],PetscInt p,PetscScalar x[],PetscCtx ctx)876 PetscErrorCode circleMultipleX(PetscInt dim, PetscReal time, const PetscReal unused[], PetscInt p, PetscScalar x[], PetscCtx ctx)
877 {
878 const PetscInt n = 5;
879 const PetscReal r0 = (p / n) + 1.;
880 const PetscReal th0 = (2. * PETSC_PI * (p % n)) / n;
881
882 x[0] = r0 * PetscCosReal(th0);
883 x[1] = r0 * PetscSinReal(th0);
884 return PETSC_SUCCESS;
885 }
886
887 /* Put 5 particles into each circle */
circleMultipleV(PetscInt dim,PetscReal time,const PetscReal unused[],PetscInt p,PetscScalar v[],PetscCtx ctx)888 PetscErrorCode circleMultipleV(PetscInt dim, PetscReal time, const PetscReal unused[], PetscInt p, PetscScalar v[], PetscCtx ctx)
889 {
890 const PetscInt n = 5;
891 const PetscReal r0 = (p / n) + 1.;
892 const PetscReal th0 = (2. * PETSC_PI * (p % n)) / n;
893 const PetscReal omega = PetscSqrtReal(1000. / r0) / r0;
894
895 v[0] = -r0 * omega * PetscSinReal(th0);
896 v[1] = r0 * omega * PetscCosReal(th0);
897 return PETSC_SUCCESS;
898 }
899
900 /*
901 InitializeSolveAndSwarm - Set the solution values to the swarm coordinates and velocities, and also possibly set the initial values.
902
903 Input Parameters:
904 + ts - The TS
905 - useInitial - Flag to also set the initial conditions to the current coordinates and velocities and setup the problem
906
907 Output Parameter:
908 . u - The initialized solution vector
909
910 Level: advanced
911
912 .seealso: InitializeSolve()
913 */
InitializeSolveAndSwarm(TS ts,PetscBool useInitial)914 static PetscErrorCode InitializeSolveAndSwarm(TS ts, PetscBool useInitial)
915 {
916 DM sw;
917 Vec u, gc, gv, gc0, gv0;
918 IS isx, isv;
919 AppCtx *user;
920
921 PetscFunctionBeginUser;
922 PetscCall(TSGetDM(ts, &sw));
923 PetscCall(DMGetApplicationContext(sw, &user));
924 if (useInitial) {
925 PetscReal v0[1] = {1.};
926
927 PetscCall(DMSwarmInitializeCoordinates(sw));
928 PetscCall(DMSwarmInitializeVelocitiesFromOptions(sw, v0));
929 PetscCall(DMSwarmMigrate(sw, PETSC_TRUE));
930 PetscCall(TSReset(ts));
931 PetscCall(CreateSolution(ts));
932 PetscCall(SetProblem(ts));
933 }
934 PetscCall(TSGetSolution(ts, &u));
935 PetscCall(TSRHSSplitGetIS(ts, "position", &isx));
936 PetscCall(TSRHSSplitGetIS(ts, "momentum", &isv));
937 PetscCall(DMSwarmCreateGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc));
938 PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "initCoordinates", &gc0));
939 if (useInitial) PetscCall(VecCopy(gc, gc0));
940 PetscCall(VecISCopy(u, isx, SCATTER_FORWARD, gc));
941 PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc));
942 PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "initCoordinates", &gc0));
943 PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "velocity", &gv));
944 PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "initVelocity", &gv0));
945 if (useInitial) PetscCall(VecCopy(gv, gv0));
946 PetscCall(VecISCopy(u, isv, SCATTER_FORWARD, gv));
947 PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "velocity", &gv));
948 PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "initVelocity", &gv0));
949 PetscFunctionReturn(PETSC_SUCCESS);
950 }
951
InitializeSolve(TS ts,Vec u)952 static PetscErrorCode InitializeSolve(TS ts, Vec u)
953 {
954 PetscFunctionBegin;
955 PetscCall(TSSetSolution(ts, u));
956 PetscCall(InitializeSolveAndSwarm(ts, PETSC_TRUE));
957 PetscFunctionReturn(PETSC_SUCCESS);
958 }
959
ComputeError(TS ts,Vec U,Vec E)960 static PetscErrorCode ComputeError(TS ts, Vec U, Vec E)
961 {
962 MPI_Comm comm;
963 DM sw;
964 AppCtx *user;
965 const PetscScalar *u;
966 const PetscReal *coords, *vel;
967 PetscScalar *e;
968 PetscReal t;
969 PetscInt dim, Np, p;
970
971 PetscFunctionBeginUser;
972 PetscCall(PetscObjectGetComm((PetscObject)ts, &comm));
973 PetscCall(TSGetDM(ts, &sw));
974 PetscCall(DMGetApplicationContext(sw, &user));
975 PetscCall(DMGetDimension(sw, &dim));
976 PetscCall(TSGetSolveTime(ts, &t));
977 PetscCall(VecGetArray(E, &e));
978 PetscCall(VecGetArrayRead(U, &u));
979 PetscCall(VecGetLocalSize(U, &Np));
980 PetscCall(DMSwarmGetField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
981 PetscCall(DMSwarmGetField(sw, "initVelocity", NULL, NULL, (void **)&vel));
982 Np /= 2 * dim;
983 for (p = 0; p < Np; ++p) {
984 /* TODO generalize initial conditions and project into plane instead of assuming x-y */
985 const PetscReal r0 = DMPlex_NormD_Internal(dim, &coords[p * dim]);
986 const PetscReal th0 = PetscAtan2Real(coords[p * dim + 1], coords[p * dim + 0]);
987 const PetscReal v0 = DMPlex_NormD_Internal(dim, &vel[p * dim]);
988 const PetscReal omega = v0 / r0;
989 const PetscReal ct = PetscCosReal(omega * t + th0);
990 const PetscReal st = PetscSinReal(omega * t + th0);
991 const PetscScalar *x = &u[(p * 2 + 0) * dim];
992 const PetscScalar *v = &u[(p * 2 + 1) * dim];
993 const PetscReal xe[3] = {r0 * ct, r0 * st, 0.0};
994 const PetscReal ve[3] = {-v0 * st, v0 * ct, 0.0};
995 PetscInt d;
996
997 for (d = 0; d < dim; ++d) {
998 e[(p * 2 + 0) * dim + d] = x[d] - xe[d];
999 e[(p * 2 + 1) * dim + d] = v[d] - ve[d];
1000 }
1001 if (user->error) {
1002 const PetscReal en = 0.5 * DMPlex_DotRealD_Internal(dim, v, v);
1003 const PetscReal exen = 0.5 * PetscSqr(v0);
1004 PetscCall(PetscPrintf(comm, "t %.4g: p%" PetscInt_FMT " error [%.2f %.2f] sol [(%.6lf %.6lf) (%.6lf %.6lf)] exact [(%.6lf %.6lf) (%.6lf %.6lf)] energy/exact energy %g / %g (%.10lf%%)\n", (double)t, p, (double)DMPlex_NormD_Internal(dim, &e[(p * 2 + 0) * dim]), (double)DMPlex_NormD_Internal(dim, &e[(p * 2 + 1) * dim]), (double)x[0], (double)x[1], (double)v[0], (double)v[1], (double)xe[0], (double)xe[1], (double)ve[0], (double)ve[1], (double)en, (double)exen, (double)(PetscAbsReal(exen - en) * 100. / exen)));
1005 }
1006 }
1007 PetscCall(DMSwarmRestoreField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
1008 PetscCall(DMSwarmRestoreField(sw, "initVelocity", NULL, NULL, (void **)&vel));
1009 PetscCall(VecRestoreArrayRead(U, &u));
1010 PetscCall(VecRestoreArray(E, &e));
1011 PetscFunctionReturn(PETSC_SUCCESS);
1012 }
1013
EnergyMonitor(TS ts,PetscInt step,PetscReal t,Vec U,PetscCtx ctx)1014 static PetscErrorCode EnergyMonitor(TS ts, PetscInt step, PetscReal t, Vec U, PetscCtx ctx)
1015 {
1016 const PetscInt ostep = ((AppCtx *)ctx)->ostep;
1017 const EMType em = ((AppCtx *)ctx)->em;
1018 DM sw;
1019 const PetscScalar *u;
1020 PetscReal *coords, *E;
1021 PetscReal enKin = 0., enEM = 0.;
1022 PetscInt dim, d, Np, p, q;
1023
1024 PetscFunctionBeginUser;
1025 if (step % ostep == 0) {
1026 PetscCall(TSGetDM(ts, &sw));
1027 PetscCall(DMGetDimension(sw, &dim));
1028 PetscCall(VecGetArrayRead(U, &u));
1029 PetscCall(VecGetLocalSize(U, &Np));
1030 Np /= 2 * dim;
1031 PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
1032 PetscCall(DMSwarmGetField(sw, "E_field", NULL, NULL, (void **)&E));
1033 if (!step) PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "Time Step Part Energy\n"));
1034 for (p = 0; p < Np; ++p) {
1035 const PetscReal v2 = DMPlex_DotRealD_Internal(dim, &u[(p * 2 + 1) * dim], &u[(p * 2 + 1) * dim]);
1036 PetscReal *pcoord = &coords[p * dim];
1037
1038 PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)ts), "%.6lf %4" PetscInt_FMT " %5" PetscInt_FMT " %10.4lf\n", (double)t, step, p, (double)(0.5 * v2)));
1039 enKin += 0.5 * v2;
1040 if (em == EM_NONE) {
1041 continue;
1042 } else if (em == EM_COULOMB) {
1043 for (q = p + 1; q < Np; ++q) {
1044 PetscReal *qcoord = &coords[q * dim];
1045 PetscReal rpq[3], r;
1046 for (d = 0; d < dim; ++d) rpq[d] = pcoord[d] - qcoord[d];
1047 r = DMPlex_NormD_Internal(dim, rpq);
1048 enEM += 1. / r;
1049 }
1050 } else if (em == EM_PRIMAL || em == EM_MIXED) {
1051 for (d = 0; d < dim; ++d) enEM += E[p * dim + d];
1052 }
1053 }
1054 PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)ts), "%.6lf %4" PetscInt_FMT " KE\t %10.4lf\n", (double)t, step, (double)enKin));
1055 PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)ts), "%.6lf %4" PetscInt_FMT " PE\t %1.10g\n", (double)t, step, (double)enEM));
1056 PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)ts), "%.6lf %4" PetscInt_FMT " E\t %10.4lf\n", (double)t, step, (double)(enKin + enEM)));
1057 PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
1058 PetscCall(DMSwarmRestoreField(sw, "E_field", NULL, NULL, (void **)&E));
1059 PetscCall(PetscSynchronizedFlush(PetscObjectComm((PetscObject)ts), NULL));
1060 PetscCall(VecRestoreArrayRead(U, &u));
1061 }
1062 PetscFunctionReturn(PETSC_SUCCESS);
1063 }
1064
SetUpMigrateParticles(TS ts,PetscInt n,PetscReal t,Vec x,PetscBool * flg,PetscCtx ctx)1065 static PetscErrorCode SetUpMigrateParticles(TS ts, PetscInt n, PetscReal t, Vec x, PetscBool *flg, PetscCtx ctx)
1066 {
1067 DM sw;
1068
1069 PetscFunctionBeginUser;
1070 *flg = PETSC_TRUE;
1071 PetscCall(TSGetDM(ts, &sw));
1072 PetscCall(DMViewFromOptions(sw, NULL, "-migrate_view_pre"));
1073 {
1074 Vec u, gc, gv;
1075 IS isx, isv;
1076
1077 PetscCall(TSGetSolution(ts, &u));
1078 PetscCall(TSRHSSplitGetIS(ts, "position", &isx));
1079 PetscCall(TSRHSSplitGetIS(ts, "momentum", &isv));
1080 PetscCall(DMSwarmCreateGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc));
1081 PetscCall(VecISCopy(u, isx, SCATTER_REVERSE, gc));
1082 PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc));
1083 PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "velocity", &gv));
1084 PetscCall(VecISCopy(u, isv, SCATTER_REVERSE, gv));
1085 PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "velocity", &gv));
1086 }
1087 PetscFunctionReturn(PETSC_SUCCESS);
1088 }
1089
MigrateParticles(TS ts,PetscInt nv,Vec vecsin[],Vec vecsout[],PetscCtx ctx)1090 static PetscErrorCode MigrateParticles(TS ts, PetscInt nv, Vec vecsin[], Vec vecsout[], PetscCtx ctx)
1091 {
1092 DM sw;
1093
1094 PetscFunctionBeginUser;
1095 PetscCall(TSGetDM(ts, &sw));
1096 PetscCall(DMSwarmMigrate(sw, PETSC_TRUE));
1097 PetscCall(CreateSolution(ts));
1098 PetscCall(SetProblem(ts));
1099 PetscCall(InitializeSolveAndSwarm(ts, PETSC_FALSE));
1100 PetscFunctionReturn(PETSC_SUCCESS);
1101 }
1102
main(int argc,char ** argv)1103 int main(int argc, char **argv)
1104 {
1105 DM dm, sw;
1106 TS ts;
1107 Vec u;
1108 AppCtx user;
1109
1110 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
1111 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
1112 PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
1113 PetscCall(CreatePoisson(dm, &user));
1114 PetscCall(CreateSwarm(dm, &user, &sw));
1115 PetscCall(DMSetApplicationContext(sw, &user));
1116
1117 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
1118 PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
1119 PetscCall(TSSetDM(ts, sw));
1120 PetscCall(TSSetMaxTime(ts, 0.1));
1121 PetscCall(TSSetTimeStep(ts, 0.00001));
1122 PetscCall(TSSetMaxSteps(ts, 100));
1123 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
1124 PetscCall(TSMonitorSet(ts, EnergyMonitor, &user, NULL));
1125 PetscCall(TSSetFromOptions(ts));
1126 PetscCall(TSSetComputeInitialCondition(ts, InitializeSolve));
1127 PetscCall(TSSetComputeExactError(ts, ComputeError));
1128 PetscCall(TSSetResize(ts, PETSC_FALSE, SetUpMigrateParticles, MigrateParticles, NULL));
1129
1130 PetscCall(CreateSolution(ts));
1131 PetscCall(TSGetSolution(ts, &u));
1132 PetscCall(TSComputeInitialCondition(ts, u));
1133 PetscCall(TSSolve(ts, NULL));
1134
1135 PetscCall(SNESDestroy(&user.snes));
1136 PetscCall(TSDestroy(&ts));
1137 PetscCall(DMDestroy(&sw));
1138 PetscCall(DMDestroy(&dm));
1139 PetscCall(PetscFinalize());
1140 return 0;
1141 }
1142
1143 /*TEST
1144
1145 build:
1146 requires: double !complex
1147
1148 testset:
1149 requires: defined(PETSC_HAVE_EXECUTABLE_EXPORT)
1150 args: -dm_plex_dim 2 -dm_plex_simplex 0 -dm_plex_box_faces 1,1 -dm_plex_box_lower -5,-5 -dm_plex_box_upper 5,5 \
1151 -dm_swarm_num_particles 2 -dm_swarm_coordinate_function circleSingleX -dm_swarm_velocity_function circleSingleV \
1152 -ts_type basicsymplectic\
1153 -dm_view -output_step 50 -ts_time_step 0.01 -ts_max_time 10.0 -ts_max_steps 10
1154 test:
1155 suffix: none_bsi_2d_1
1156 args: -ts_basicsymplectic_type 1 -em_type none -error
1157 test:
1158 suffix: none_bsi_2d_2
1159 args: -ts_basicsymplectic_type 2 -em_type none -error
1160 test:
1161 suffix: none_bsi_2d_3
1162 args: -ts_basicsymplectic_type 3 -em_type none -error
1163 test:
1164 suffix: none_bsi_2d_4
1165 args: -ts_basicsymplectic_type 4 -em_type none -error
1166 test:
1167 suffix: coulomb_bsi_2d_1
1168 args: -ts_basicsymplectic_type 1
1169 test:
1170 suffix: coulomb_bsi_2d_2
1171 args: -ts_basicsymplectic_type 2
1172 test:
1173 suffix: coulomb_bsi_2d_3
1174 args: -ts_basicsymplectic_type 3
1175 test:
1176 suffix: coulomb_bsi_2d_4
1177 args: -ts_basicsymplectic_type 4
1178
1179 testset:
1180 requires: defined(PETSC_HAVE_EXECUTABLE_EXPORT)
1181 args: -dm_plex_dim 2 -dm_plex_simplex 0 -dm_plex_box_faces 1,1 -dm_plex_box_lower -5,-5 -dm_plex_box_upper 5,5 \
1182 -dm_swarm_num_particles 2 -dm_swarm_coordinate_function circleSingleX -dm_swarm_velocity_function circleSingleV \
1183 -ts_type basicsymplectic\
1184 -em_type primal -em_pc_type svd\
1185 -dm_view -output_step 50 -error -ts_time_step 0.01 -ts_max_time 10.0 -ts_max_steps 10\
1186 -petscspace_degree 2 -petscfe_default_quadrature_order 3 -sigma 1.0e-8 -timeScale 2.0e-14
1187 test:
1188 suffix: poisson_bsi_2d_1
1189 args: -ts_basicsymplectic_type 1
1190 test:
1191 suffix: poisson_bsi_2d_2
1192 args: -ts_basicsymplectic_type 2
1193 test:
1194 suffix: poisson_bsi_2d_3
1195 args: -ts_basicsymplectic_type 3
1196 test:
1197 suffix: poisson_bsi_2d_4
1198 args: -ts_basicsymplectic_type 4
1199
1200 testset:
1201 requires: defined(PETSC_HAVE_EXECUTABLE_EXPORT)
1202 args: -dm_plex_dim 2 -dm_plex_simplex 0 -dm_plex_box_faces 1,1 -dm_plex_box_lower -5,-5 -dm_plex_box_upper 5,5 \
1203 -dm_swarm_num_particles 2 -dm_swarm_coordinate_function circleSingleX -dm_swarm_velocity_function circleSingleV \
1204 -ts_convergence_estimate -convest_num_refine 2 -em_type primal \
1205 -mat_type baij -em_ksp_error_if_not_converged -em_pc_type svd \
1206 -dm_view -output_step 50 -error -ts_time_step 0.01 -ts_max_time 10.0 -ts_max_steps 10 \
1207 -sigma 1.0e-8 -timeScale 2.0e-14
1208 test:
1209 suffix: im_2d_0
1210 args: -ts_type theta -ts_theta_theta 0.5
1211 test:
1212 suffix: dg_2d_none
1213 args: -ts_type discgrad -ts_discgrad_type none -snes_type qn
1214 test:
1215 suffix: dg_2d_average
1216 args: -ts_type discgrad -ts_discgrad_type average -snes_type qn
1217 test:
1218 suffix: dg_2d_gonzalez
1219 args: -ts_type discgrad -ts_discgrad_type gonzalez -snes_fd -snes_type newtonls -snes_fd -pc_type lu
1220
1221 testset:
1222 requires: defined(PETSC_HAVE_EXECUTABLE_EXPORT)
1223 args: -dm_plex_dim 2 -dm_plex_simplex 0 -dm_plex_box_faces 10,10 -dm_plex_box_lower -5,-5 -dm_plex_box_upper 5,5 -petscpartitioner_type simple \
1224 -dm_swarm_num_particles 2 -dm_swarm_coordinate_function circleSingleX -dm_swarm_velocity_function circleSingleV -dm_swarm_num_species 1\
1225 -ts_type basicsymplectic -ts_convergence_estimate -convest_num_refine 2 \
1226 -em_snes_type ksponly -em_pc_type svd -em_type primal -petscspace_degree 1\
1227 -dm_view -output_step 50\
1228 -pc_type svd -sigma 1.0e-8 -timeScale 2.0e-14 -ts_time_step 0.01 -ts_max_time 1.0 -ts_max_steps 10
1229 test:
1230 suffix: bsi_2d_mesh_1
1231 args: -ts_basicsymplectic_type 4
1232 test:
1233 suffix: bsi_2d_mesh_1_par_2
1234 nsize: 2
1235 args: -ts_basicsymplectic_type 4
1236 test:
1237 suffix: bsi_2d_mesh_1_par_3
1238 nsize: 3
1239 args: -ts_basicsymplectic_type 4
1240 test:
1241 suffix: bsi_2d_mesh_1_par_4
1242 nsize: 4
1243 args: -ts_basicsymplectic_type 4 -dm_swarm_num_particles 0,0,2,0
1244
1245 testset:
1246 requires: defined(PETSC_HAVE_EXECUTABLE_EXPORT)
1247 args: -dm_plex_dim 2 -dm_plex_simplex 0 -dm_plex_box_faces 10,10 -dm_plex_box_lower -5,-5 -dm_plex_box_upper 5,5 \
1248 -dm_swarm_num_particles 10 -dm_swarm_coordinate_function circleMultipleX -dm_swarm_velocity_function circleMultipleV \
1249 -ts_convergence_estimate -convest_num_refine 2 \
1250 -em_pc_type lu\
1251 -dm_view -output_step 50 -error\
1252 -sigma 1.0e-8 -timeScale 2.0e-14 -ts_time_step 0.01 -ts_max_time 10.0 -ts_max_steps 10
1253 test:
1254 suffix: bsi_2d_multiple_1
1255 args: -ts_type basicsymplectic -ts_basicsymplectic_type 1
1256 test:
1257 suffix: bsi_2d_multiple_2
1258 args: -ts_type basicsymplectic -ts_basicsymplectic_type 2
1259 test:
1260 suffix: bsi_2d_multiple_3
1261 args: -ts_type basicsymplectic -ts_basicsymplectic_type 3 -ts_time_step 0.001
1262 test:
1263 suffix: im_2d_multiple_0
1264 args: -ts_type theta -ts_theta_theta 0.5 \
1265 -mat_type baij -em_ksp_error_if_not_converged -em_pc_type lu
1266
1267 testset:
1268 requires: defined(PETSC_HAVE_EXECUTABLE_EXPORT)
1269 args: -dm_plex_dim 2 -dm_plex_simplex 0 -dm_plex_box_faces 2,2 -dm_plex_box_lower -5,-5 -dm_plex_box_upper 5,5 \
1270 -dm_swarm_num_particles 2 -dm_swarm_coordinate_function circleSingleX -dm_swarm_velocity_function circleSingleV \
1271 -em_pc_type fieldsplit -ksp_rtol 1e-10 -em_ksp_type preonly -em_type mixed -em_ksp_error_if_not_converged\
1272 -dm_view -output_step 50 -error -dm_refine 0\
1273 -pc_type svd -sigma 1.0e-8 -timeScale 2.0e-14 -ts_time_step 0.01 -ts_max_time 10.0 -ts_max_steps 10
1274 test:
1275 suffix: bsi_4_rt_1
1276 args: -ts_type basicsymplectic -ts_basicsymplectic_type 4\
1277 -pc_fieldsplit_detect_saddle_point\
1278 -pc_type fieldsplit\
1279 -pc_fieldsplit_type schur\
1280 -pc_fieldsplit_schur_precondition full \
1281 -field_petscspace_degree 2\
1282 -field_petscfe_default_quadrature_order 1\
1283 -field_petscspace_type sum \
1284 -field_petscspace_variables 2 \
1285 -field_petscspace_components 2 \
1286 -field_petscspace_sum_spaces 2 \
1287 -field_petscspace_sum_concatenate true \
1288 -field_sumcomp_0_petscspace_variables 2 \
1289 -field_sumcomp_0_petscspace_type tensor \
1290 -field_sumcomp_0_petscspace_tensor_spaces 2 \
1291 -field_sumcomp_0_petscspace_tensor_uniform false \
1292 -field_sumcomp_0_tensorcomp_0_petscspace_degree 1 \
1293 -field_sumcomp_0_tensorcomp_1_petscspace_degree 0 \
1294 -field_sumcomp_1_petscspace_variables 2 \
1295 -field_sumcomp_1_petscspace_type tensor \
1296 -field_sumcomp_1_petscspace_tensor_spaces 2 \
1297 -field_sumcomp_1_petscspace_tensor_uniform false \
1298 -field_sumcomp_1_tensorcomp_0_petscspace_degree 0 \
1299 -field_sumcomp_1_tensorcomp_1_petscspace_degree 1 \
1300 -field_petscdualspace_form_degree -1 \
1301 -field_petscdualspace_order 1 \
1302 -field_petscdualspace_lagrange_trimmed true\
1303 -ksp_gmres_restart 500
1304
1305 TEST*/
1306