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 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 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 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 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 */ 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 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> */ 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 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 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 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 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 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 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 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 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 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 520 static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec G, void *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 */ 567 static PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat J, Mat P, void *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 599 static PetscErrorCode RHSFunctionX(TS ts, PetscReal t, Vec V, Vec Xres, void *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 621 static PetscErrorCode RHSFunctionV(TS ts, PetscReal t, Vec X, Vec Vres, void *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) */ 661 PetscErrorCode RHSJacobianS(TS ts, PetscReal t, Vec U, Mat S, void *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 684 PetscErrorCode RHSObjectiveF(TS ts, PetscReal t, Vec U, PetscScalar *F, void *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 */ 741 PetscErrorCode RHSFunctionG(TS ts, PetscReal t, Vec U, Vec G, void *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 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 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, (void **)&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 861 PetscErrorCode circleSingleX(PetscInt dim, PetscReal time, const PetscReal dummy[], PetscInt p, PetscScalar x[], void *ctx) 862 { 863 x[0] = p + 1.; 864 x[1] = 0.; 865 return PETSC_SUCCESS; 866 } 867 868 PetscErrorCode circleSingleV(PetscInt dim, PetscReal time, const PetscReal dummy[], PetscInt p, PetscScalar v[], void *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 */ 876 PetscErrorCode circleMultipleX(PetscInt dim, PetscReal time, const PetscReal dummy[], PetscInt p, PetscScalar x[], void *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 */ 888 PetscErrorCode circleMultipleV(PetscInt dim, PetscReal time, const PetscReal dummy[], PetscInt p, PetscScalar v[], void *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 */ 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 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 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 1014 static PetscErrorCode EnergyMonitor(TS ts, PetscInt step, PetscReal t, Vec U, void *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 1065 static PetscErrorCode SetUpMigrateParticles(TS ts, PetscInt n, PetscReal t, Vec x, PetscBool *flg, void *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 1090 static PetscErrorCode MigrateParticles(TS ts, PetscInt nv, Vec vecsin[], Vec vecsout[], void *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 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