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_dt 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, user, user, 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 PetscFunctionReturn(PETSC_SUCCESS); 225 } 226 227 static PetscErrorCode ComputeFieldAtParticles_Coulomb(SNES snes, DM sw, PetscReal E[]) 228 { 229 PetscReal *coords; 230 PetscInt dim, d, Np, p, q; 231 PetscMPIInt size; 232 233 PetscFunctionBegin; 234 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)snes), &size)); 235 PetscCheck(size == 1, PetscObjectComm((PetscObject)snes), PETSC_ERR_SUP, "Coulomb code only works in serial"); 236 PetscCall(DMGetDimension(sw, &dim)); 237 PetscCall(DMSwarmGetLocalSize(sw, &Np)); 238 239 PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords)); 240 for (p = 0; p < Np; ++p) { 241 PetscReal *pcoord = &coords[p * dim]; 242 PetscReal *pE = &E[p * dim]; 243 /* Calculate field at particle p due to particle q */ 244 for (q = 0; q < Np; ++q) { 245 PetscReal *qcoord = &coords[q * dim]; 246 PetscReal rpq[3], r; 247 248 if (p == q) continue; 249 for (d = 0; d < dim; ++d) rpq[d] = pcoord[d] - qcoord[d]; 250 r = DMPlex_NormD_Internal(dim, rpq); 251 for (d = 0; d < dim; ++d) pE[d] += rpq[d] / PetscPowRealInt(r, 3); 252 } 253 } 254 PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords)); 255 PetscFunctionReturn(PETSC_SUCCESS); 256 } 257 258 static PetscErrorCode ComputeFieldAtParticles_Primal(SNES snes, DM sw, PetscReal E[]) 259 { 260 DM dm; 261 PetscDS ds; 262 PetscFE fe; 263 Mat M_p; 264 Vec phi, locPhi, rho, f; 265 PetscReal *coords, chargeTol = 1e-13; 266 PetscInt dim, d, cStart, cEnd, c, Np; 267 268 PetscFunctionBegin; 269 PetscCall(DMGetDimension(sw, &dim)); 270 PetscCall(DMSwarmGetLocalSize(sw, &Np)); 271 272 /* Create the charges rho */ 273 PetscCall(SNESGetDM(snes, &dm)); 274 PetscCall(DMCreateMassMatrix(sw, dm, &M_p)); 275 PetscCall(DMGetGlobalVector(dm, &rho)); 276 PetscCall(PetscObjectSetName((PetscObject)rho, "rho")); 277 PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "w_q", &f)); 278 PetscCall(PetscObjectSetName((PetscObject)f, "particle weight")); 279 PetscCall(MatMultTranspose(M_p, f, rho)); 280 PetscCall(MatViewFromOptions(M_p, NULL, "-mp_view")); 281 PetscCall(VecViewFromOptions(f, NULL, "-weights_view")); 282 PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f)); 283 PetscCall(MatDestroy(&M_p)); 284 { 285 PetscScalar sum; 286 PetscInt n; 287 PetscReal phi_0 = 1.; /* (sigma*sigma*sigma)*(timeScale*timeScale)/(m_e*q_e*epsi_0)*/ 288 289 /* Remove constant from rho */ 290 PetscCall(VecGetSize(rho, &n)); 291 PetscCall(VecSum(rho, &sum)); 292 PetscCall(VecShift(rho, -sum / n)); 293 PetscCall(VecSum(rho, &sum)); 294 PetscCheck(PetscAbsScalar(sum) < chargeTol, PetscObjectComm((PetscObject)sw), PETSC_ERR_PLIB, "Charge should have no DC component: %g", (double)PetscRealPart(sum)); 295 /* Nondimensionalize rho */ 296 PetscCall(VecScale(rho, phi_0)); 297 } 298 PetscCall(VecViewFromOptions(rho, NULL, "-poisson_rho_view")); 299 300 PetscCall(DMGetGlobalVector(dm, &phi)); 301 PetscCall(PetscObjectSetName((PetscObject)phi, "potential")); 302 PetscCall(VecSet(phi, 0.0)); 303 PetscCall(SNESSolve(snes, rho, phi)); 304 PetscCall(DMRestoreGlobalVector(dm, &rho)); 305 PetscCall(VecViewFromOptions(phi, NULL, "-phi_view")); 306 307 PetscCall(DMGetLocalVector(dm, &locPhi)); 308 PetscCall(DMGlobalToLocalBegin(dm, phi, INSERT_VALUES, locPhi)); 309 PetscCall(DMGlobalToLocalEnd(dm, phi, INSERT_VALUES, locPhi)); 310 PetscCall(DMRestoreGlobalVector(dm, &phi)); 311 312 PetscCall(DMGetDS(dm, &ds)); 313 PetscCall(PetscDSGetDiscretization(ds, 0, (PetscObject *)&fe)); 314 PetscCall(DMSwarmSortGetAccess(sw)); 315 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 316 PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords)); 317 for (c = cStart; c < cEnd; ++c) { 318 PetscTabulation tab; 319 PetscScalar *clPhi = NULL; 320 PetscReal *pcoord, *refcoord; 321 PetscReal v[3], J[9], invJ[9], detJ; 322 PetscInt *points; 323 PetscInt Ncp, cp; 324 325 PetscCall(DMSwarmSortGetPointsPerCell(sw, c, &Ncp, &points)); 326 PetscCall(DMGetWorkArray(dm, Ncp * dim, MPIU_REAL, &pcoord)); 327 PetscCall(DMGetWorkArray(dm, Ncp * dim, MPIU_REAL, &refcoord)); 328 for (cp = 0; cp < Ncp; ++cp) 329 for (d = 0; d < dim; ++d) pcoord[cp * dim + d] = coords[points[cp] * dim + d]; 330 PetscCall(DMPlexCoordinatesToReference(dm, c, Ncp, pcoord, refcoord)); 331 PetscCall(PetscFECreateTabulation(fe, 1, Ncp, refcoord, 1, &tab)); 332 PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v, J, invJ, &detJ)); 333 PetscCall(DMPlexVecGetClosure(dm, NULL, locPhi, c, NULL, &clPhi)); 334 for (cp = 0; cp < Ncp; ++cp) { 335 const PetscReal *basisDer = tab->T[1]; 336 const PetscInt p = points[cp]; 337 338 for (d = 0; d < dim; ++d) E[p * dim + d] = 0.; 339 PetscCall(PetscFEFreeInterpolateGradient_Static(fe, basisDer, clPhi, dim, invJ, NULL, cp, &E[p * dim])); 340 for (d = 0; d < dim; ++d) E[p * dim + d] *= -1.0; 341 } 342 PetscCall(DMPlexVecRestoreClosure(dm, NULL, locPhi, c, NULL, &clPhi)); 343 PetscCall(DMRestoreWorkArray(dm, Ncp * dim, MPIU_REAL, &pcoord)); 344 PetscCall(DMRestoreWorkArray(dm, Ncp * dim, MPIU_REAL, &refcoord)); 345 PetscCall(PetscTabulationDestroy(&tab)); 346 PetscCall(PetscFree(points)); 347 } 348 PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords)); 349 PetscCall(DMSwarmSortRestoreAccess(sw)); 350 PetscCall(DMRestoreLocalVector(dm, &locPhi)); 351 PetscFunctionReturn(PETSC_SUCCESS); 352 } 353 354 static PetscErrorCode ComputeFieldAtParticles_Mixed(SNES snes, DM sw, PetscReal E[]) 355 { 356 DM dm, potential_dm; 357 IS potential_IS; 358 PetscDS ds; 359 PetscFE fe; 360 PetscFEGeom feGeometry; 361 Mat M_p; 362 Vec phi, locPhi, rho, f, temp_rho; 363 PetscQuadrature q; 364 PetscReal *coords, chargeTol = 1e-13; 365 PetscInt dim, d, cStart, cEnd, c, Np, fields = 1; 366 367 PetscFunctionBegin; 368 PetscCall(DMGetDimension(sw, &dim)); 369 PetscCall(DMSwarmGetLocalSize(sw, &Np)); 370 371 /* Create the charges rho */ 372 PetscCall(SNESGetDM(snes, &dm)); 373 PetscCall(DMGetGlobalVector(dm, &rho)); 374 PetscCall(PetscObjectSetName((PetscObject)rho, "rho")); 375 376 PetscCall(DMCreateSubDM(dm, 1, &fields, &potential_IS, &potential_dm)); 377 PetscCall(DMCreateMassMatrix(sw, potential_dm, &M_p)); 378 PetscCall(MatViewFromOptions(M_p, NULL, "-mp_view")); 379 PetscCall(DMGetGlobalVector(potential_dm, &temp_rho)); 380 PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "w_q", &f)); 381 PetscCall(PetscObjectSetName((PetscObject)f, "particle weight")); 382 PetscCall(VecViewFromOptions(f, NULL, "-weights_view")); 383 PetscCall(MatMultTranspose(M_p, f, temp_rho)); 384 PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f)); 385 PetscCall(MatDestroy(&M_p)); 386 PetscCall(PetscObjectSetName((PetscObject)rho, "rho")); 387 PetscCall(VecViewFromOptions(rho, NULL, "-poisson_rho_view")); 388 PetscCall(VecISCopy(rho, potential_IS, SCATTER_FORWARD, temp_rho)); 389 PetscCall(DMRestoreGlobalVector(potential_dm, &temp_rho)); 390 PetscCall(DMDestroy(&potential_dm)); 391 PetscCall(ISDestroy(&potential_IS)); 392 { 393 PetscScalar sum; 394 PetscInt n; 395 PetscReal phi_0 = 1.; /*(sigma*sigma*sigma)*(timeScale*timeScale)/(m_e*q_e*epsi_0);*/ 396 397 /*Remove constant from rho*/ 398 PetscCall(VecGetSize(rho, &n)); 399 PetscCall(VecSum(rho, &sum)); 400 PetscCall(VecShift(rho, -sum / n)); 401 PetscCall(VecSum(rho, &sum)); 402 PetscCheck(PetscAbsScalar(sum) < chargeTol, PetscObjectComm((PetscObject)sw), PETSC_ERR_PLIB, "Charge should have no DC component: %g", (double)PetscRealPart(sum)); 403 /* Nondimensionalize rho */ 404 PetscCall(VecScale(rho, phi_0)); 405 } 406 PetscCall(DMGetGlobalVector(dm, &phi)); 407 PetscCall(PetscObjectSetName((PetscObject)phi, "potential")); 408 PetscCall(VecSet(phi, 0.0)); 409 PetscCall(SNESSolve(snes, NULL, phi)); 410 PetscCall(DMRestoreGlobalVector(dm, &rho)); 411 PetscCall(VecViewFromOptions(phi, NULL, "-phi_view")); 412 413 PetscCall(DMGetLocalVector(dm, &locPhi)); 414 PetscCall(DMGlobalToLocalBegin(dm, phi, INSERT_VALUES, locPhi)); 415 PetscCall(DMGlobalToLocalEnd(dm, phi, INSERT_VALUES, locPhi)); 416 PetscCall(DMRestoreGlobalVector(dm, &phi)); 417 418 PetscCall(DMGetDS(dm, &ds)); 419 PetscCall(PetscDSGetDiscretization(ds, 0, (PetscObject *)&fe)); 420 PetscCall(DMSwarmSortGetAccess(sw)); 421 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 422 PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords)); 423 for (c = cStart; c < cEnd; ++c) { 424 PetscTabulation tab; 425 PetscScalar *clPhi = NULL; 426 PetscReal *pcoord, *refcoord; 427 PetscReal v[3], J[9], invJ[9], detJ; 428 PetscInt *points; 429 PetscInt Ncp, cp; 430 431 PetscCall(DMSwarmSortGetPointsPerCell(sw, c, &Ncp, &points)); 432 PetscCall(DMGetWorkArray(dm, Ncp * dim, MPIU_REAL, &pcoord)); 433 PetscCall(DMGetWorkArray(dm, Ncp * dim, MPIU_REAL, &refcoord)); 434 for (cp = 0; cp < Ncp; ++cp) 435 for (d = 0; d < dim; ++d) pcoord[cp * dim + d] = coords[points[cp] * dim + d]; 436 PetscCall(DMPlexCoordinatesToReference(dm, c, Ncp, pcoord, refcoord)); 437 PetscCall(PetscFECreateTabulation(fe, 1, Ncp, refcoord, 1, &tab)); 438 PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v, J, invJ, &detJ)); 439 PetscCall(DMPlexVecGetClosure(dm, NULL, locPhi, c, NULL, &clPhi)); 440 for (cp = 0; cp < Ncp; ++cp) { 441 const PetscInt p = points[cp]; 442 443 for (d = 0; d < dim; ++d) E[p * dim + d] = 0.; 444 PetscCall(PetscFEGetQuadrature(fe, &q)); 445 PetscCall(PetscFECreateCellGeometry(fe, q, &feGeometry)); 446 PetscCall(PetscFEInterpolateAtPoints_Static(fe, tab, clPhi, &feGeometry, cp, &E[p * dim])); 447 PetscCall(PetscFEDestroyCellGeometry(fe, &feGeometry)); 448 } 449 PetscCall(DMPlexVecRestoreClosure(dm, NULL, locPhi, c, NULL, &clPhi)); 450 PetscCall(DMRestoreWorkArray(dm, Ncp * dim, MPIU_REAL, &pcoord)); 451 PetscCall(DMRestoreWorkArray(dm, Ncp * dim, MPIU_REAL, &refcoord)); 452 PetscCall(PetscTabulationDestroy(&tab)); 453 PetscCall(PetscFree(points)); 454 } 455 PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords)); 456 PetscCall(DMSwarmSortRestoreAccess(sw)); 457 PetscCall(DMRestoreLocalVector(dm, &locPhi)); 458 PetscFunctionReturn(PETSC_SUCCESS); 459 } 460 461 static PetscErrorCode ComputeFieldAtParticles(SNES snes, DM sw, PetscReal E[]) 462 { 463 AppCtx *ctx; 464 PetscInt dim, Np; 465 466 PetscFunctionBegin; 467 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 468 PetscValidHeaderSpecific(sw, DM_CLASSID, 2); 469 PetscValidRealPointer(E, 3); 470 PetscCall(DMGetDimension(sw, &dim)); 471 PetscCall(DMSwarmGetLocalSize(sw, &Np)); 472 PetscCall(DMGetApplicationContext(sw, &ctx)); 473 PetscCall(PetscArrayzero(E, Np * dim)); 474 475 switch (ctx->em) { 476 case EM_PRIMAL: 477 PetscCall(ComputeFieldAtParticles_Primal(snes, sw, E)); 478 break; 479 case EM_COULOMB: 480 PetscCall(ComputeFieldAtParticles_Coulomb(snes, sw, E)); 481 break; 482 case EM_MIXED: 483 PetscCall(ComputeFieldAtParticles_Mixed(snes, sw, E)); 484 break; 485 case EM_NONE: 486 break; 487 default: 488 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No solver for electrostatic model %s", EMTypes[ctx->em]); 489 } 490 491 PetscFunctionReturn(PETSC_SUCCESS); 492 } 493 494 static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec G, void *ctx) 495 { 496 DM sw; 497 SNES snes = ((AppCtx *)ctx)->snes; 498 const PetscReal *coords, *vel; 499 const PetscScalar *u; 500 PetscScalar *g; 501 PetscReal *E; 502 PetscInt dim, d, Np, p; 503 504 PetscFunctionBeginUser; 505 PetscCall(TSGetDM(ts, &sw)); 506 PetscCall(DMGetDimension(sw, &dim)); 507 PetscCall(DMSwarmGetField(sw, "initCoordinates", NULL, NULL, (void **)&coords)); 508 PetscCall(DMSwarmGetField(sw, "initVelocity", NULL, NULL, (void **)&vel)); 509 PetscCall(DMSwarmGetField(sw, "E_field", NULL, NULL, (void **)&E)); 510 PetscCall(VecGetLocalSize(U, &Np)); 511 PetscCall(VecGetArrayRead(U, &u)); 512 PetscCall(VecGetArray(G, &g)); 513 514 PetscCall(ComputeFieldAtParticles(snes, sw, E)); 515 516 Np /= 2 * dim; 517 for (p = 0; p < Np; ++p) { 518 const PetscReal x0 = coords[p * dim + 0]; 519 const PetscReal vy0 = vel[p * dim + 1]; 520 const PetscReal omega = vy0 / x0; 521 522 for (d = 0; d < dim; ++d) { 523 g[(p * 2 + 0) * dim + d] = u[(p * 2 + 1) * dim + d]; 524 g[(p * 2 + 1) * dim + d] = E[p * dim + d] - PetscSqr(omega) * u[(p * 2 + 0) * dim + d]; 525 } 526 } 527 PetscCall(DMSwarmRestoreField(sw, "initCoordinates", NULL, NULL, (void **)&coords)); 528 PetscCall(DMSwarmRestoreField(sw, "initVelocity", NULL, NULL, (void **)&vel)); 529 PetscCall(DMSwarmRestoreField(sw, "E_field", NULL, NULL, (void **)&E)); 530 PetscCall(VecRestoreArrayRead(U, &u)); 531 PetscCall(VecRestoreArray(G, &g)); 532 PetscFunctionReturn(PETSC_SUCCESS); 533 } 534 535 /* J_{ij} = dF_i/dx_j 536 J_p = ( 0 1) 537 (-w^2 0) 538 TODO Now there is another term with w^2 from the electric field. I think we will need to invert the operator. 539 Perhaps we can approximate the Jacobian using only the cellwise P-P gradient from Coulomb 540 */ 541 static PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat J, Mat P, void *ctx) 542 { 543 DM sw; 544 const PetscReal *coords, *vel; 545 PetscInt dim, d, Np, p, rStart; 546 547 PetscFunctionBeginUser; 548 PetscCall(TSGetDM(ts, &sw)); 549 PetscCall(DMGetDimension(sw, &dim)); 550 PetscCall(VecGetLocalSize(U, &Np)); 551 PetscCall(MatGetOwnershipRange(J, &rStart, NULL)); 552 PetscCall(DMSwarmGetField(sw, "initCoordinates", NULL, NULL, (void **)&coords)); 553 PetscCall(DMSwarmGetField(sw, "initVelocity", NULL, NULL, (void **)&vel)); 554 Np /= 2 * dim; 555 for (p = 0; p < Np; ++p) { 556 const PetscReal x0 = coords[p * dim + 0]; 557 const PetscReal vy0 = vel[p * dim + 1]; 558 const PetscReal omega = vy0 / x0; 559 PetscScalar vals[4] = {0., 1., -PetscSqr(omega), 0.}; 560 561 for (d = 0; d < dim; ++d) { 562 const PetscInt rows[2] = {(p * 2 + 0) * dim + d + rStart, (p * 2 + 1) * dim + d + rStart}; 563 PetscCall(MatSetValues(J, 2, rows, 2, rows, vals, INSERT_VALUES)); 564 } 565 } 566 PetscCall(DMSwarmRestoreField(sw, "initCoordinates", NULL, NULL, (void **)&coords)); 567 PetscCall(DMSwarmRestoreField(sw, "initVelocity", NULL, NULL, (void **)&vel)); 568 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY)); 569 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY)); 570 PetscFunctionReturn(PETSC_SUCCESS); 571 } 572 573 static PetscErrorCode RHSFunctionX(TS ts, PetscReal t, Vec V, Vec Xres, void *ctx) 574 { 575 DM sw; 576 const PetscScalar *v; 577 PetscScalar *xres; 578 PetscInt Np, p, dim, d; 579 580 PetscFunctionBeginUser; 581 PetscCall(TSGetDM(ts, &sw)); 582 PetscCall(DMGetDimension(sw, &dim)); 583 PetscCall(VecGetLocalSize(Xres, &Np)); 584 Np /= dim; 585 PetscCall(VecGetArrayRead(V, &v)); 586 PetscCall(VecGetArray(Xres, &xres)); 587 for (p = 0; p < Np; ++p) { 588 for (d = 0; d < dim; ++d) xres[p * dim + d] = v[p * dim + d]; 589 } 590 PetscCall(VecRestoreArrayRead(V, &v)); 591 PetscCall(VecRestoreArray(Xres, &xres)); 592 PetscFunctionReturn(PETSC_SUCCESS); 593 } 594 595 static PetscErrorCode RHSFunctionV(TS ts, PetscReal t, Vec X, Vec Vres, void *ctx) 596 { 597 DM sw; 598 SNES snes = ((AppCtx *)ctx)->snes; 599 const PetscScalar *x; 600 const PetscReal *coords, *vel; 601 PetscScalar *vres; 602 PetscReal *E; 603 PetscInt Np, p, dim, d; 604 605 PetscFunctionBeginUser; 606 PetscCall(TSGetDM(ts, &sw)); 607 PetscCall(DMGetDimension(sw, &dim)); 608 PetscCall(DMSwarmGetField(sw, "initCoordinates", NULL, NULL, (void **)&coords)); 609 PetscCall(DMSwarmGetField(sw, "initVelocity", NULL, NULL, (void **)&vel)); 610 PetscCall(DMSwarmGetField(sw, "E_field", NULL, NULL, (void **)&E)); 611 PetscCall(VecGetLocalSize(Vres, &Np)); 612 PetscCall(VecGetArrayRead(X, &x)); 613 PetscCall(VecGetArray(Vres, &vres)); 614 PetscCheck(dim == 2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Dimension must be 2"); 615 616 PetscCall(ComputeFieldAtParticles(snes, sw, E)); 617 618 Np /= dim; 619 for (p = 0; p < Np; ++p) { 620 const PetscReal x0 = coords[p * dim + 0]; 621 const PetscReal vy0 = vel[p * dim + 1]; 622 const PetscReal omega = vy0 / x0; 623 624 for (d = 0; d < dim; ++d) vres[p * dim + d] = E[p * dim + d] - PetscSqr(omega) * x[p * dim + d]; 625 } 626 PetscCall(VecRestoreArrayRead(X, &x)); 627 PetscCall(VecRestoreArray(Vres, &vres)); 628 PetscCall(DMSwarmRestoreField(sw, "initCoordinates", NULL, NULL, (void **)&coords)); 629 PetscCall(DMSwarmRestoreField(sw, "initVelocity", NULL, NULL, (void **)&vel)); 630 PetscCall(DMSwarmRestoreField(sw, "E_field", NULL, NULL, (void **)&E)); 631 PetscFunctionReturn(PETSC_SUCCESS); 632 } 633 634 static PetscErrorCode CreateSolution(TS ts) 635 { 636 DM sw; 637 Vec u; 638 PetscInt dim, Np; 639 640 PetscFunctionBegin; 641 PetscCall(TSGetDM(ts, &sw)); 642 PetscCall(DMGetDimension(sw, &dim)); 643 PetscCall(DMSwarmGetLocalSize(sw, &Np)); 644 PetscCall(VecCreate(PETSC_COMM_WORLD, &u)); 645 PetscCall(VecSetBlockSize(u, dim)); 646 PetscCall(VecSetSizes(u, 2 * Np * dim, PETSC_DECIDE)); 647 PetscCall(VecSetUp(u)); 648 PetscCall(TSSetSolution(ts, u)); 649 PetscCall(VecDestroy(&u)); 650 PetscFunctionReturn(PETSC_SUCCESS); 651 } 652 653 static PetscErrorCode SetProblem(TS ts) 654 { 655 AppCtx *user; 656 DM sw; 657 658 PetscFunctionBegin; 659 PetscCall(TSGetDM(ts, &sw)); 660 PetscCall(DMGetApplicationContext(sw, (void **)&user)); 661 // Define unified system for (X, V) 662 { 663 Mat J; 664 PetscInt dim, Np; 665 666 PetscCall(DMGetDimension(sw, &dim)); 667 PetscCall(DMSwarmGetLocalSize(sw, &Np)); 668 PetscCall(MatCreate(PETSC_COMM_WORLD, &J)); 669 PetscCall(MatSetSizes(J, 2 * Np * dim, 2 * Np * dim, PETSC_DECIDE, PETSC_DECIDE)); 670 PetscCall(MatSetBlockSize(J, 2 * dim)); 671 PetscCall(MatSetFromOptions(J)); 672 PetscCall(MatSetUp(J)); 673 PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, user)); 674 PetscCall(TSSetRHSJacobian(ts, J, J, RHSJacobian, user)); 675 PetscCall(MatDestroy(&J)); 676 } 677 /* Define split system for X and V */ 678 { 679 Vec u; 680 IS isx, isv, istmp; 681 const PetscInt *idx; 682 PetscInt dim, Np, rstart; 683 684 PetscCall(TSGetSolution(ts, &u)); 685 PetscCall(DMGetDimension(sw, &dim)); 686 PetscCall(DMSwarmGetLocalSize(sw, &Np)); 687 PetscCall(VecGetOwnershipRange(u, &rstart, NULL)); 688 PetscCall(ISCreateStride(PETSC_COMM_WORLD, Np, (rstart / dim) + 0, 2, &istmp)); 689 PetscCall(ISGetIndices(istmp, &idx)); 690 PetscCall(ISCreateBlock(PETSC_COMM_WORLD, dim, Np, idx, PETSC_COPY_VALUES, &isx)); 691 PetscCall(ISRestoreIndices(istmp, &idx)); 692 PetscCall(ISDestroy(&istmp)); 693 PetscCall(ISCreateStride(PETSC_COMM_WORLD, Np, (rstart / dim) + 1, 2, &istmp)); 694 PetscCall(ISGetIndices(istmp, &idx)); 695 PetscCall(ISCreateBlock(PETSC_COMM_WORLD, dim, Np, idx, PETSC_COPY_VALUES, &isv)); 696 PetscCall(ISRestoreIndices(istmp, &idx)); 697 PetscCall(ISDestroy(&istmp)); 698 PetscCall(TSRHSSplitSetIS(ts, "position", isx)); 699 PetscCall(TSRHSSplitSetIS(ts, "momentum", isv)); 700 PetscCall(ISDestroy(&isx)); 701 PetscCall(ISDestroy(&isv)); 702 PetscCall(TSRHSSplitSetRHSFunction(ts, "position", NULL, RHSFunctionX, user)); 703 PetscCall(TSRHSSplitSetRHSFunction(ts, "momentum", NULL, RHSFunctionV, user)); 704 } 705 PetscFunctionReturn(PETSC_SUCCESS); 706 } 707 708 static PetscErrorCode DMSwarmTSRedistribute(TS ts) 709 { 710 DM sw; 711 Vec u; 712 PetscReal t, maxt, dt; 713 PetscInt n, maxn; 714 715 PetscFunctionBegin; 716 PetscCall(TSGetDM(ts, &sw)); 717 PetscCall(TSGetTime(ts, &t)); 718 PetscCall(TSGetMaxTime(ts, &maxt)); 719 PetscCall(TSGetTimeStep(ts, &dt)); 720 PetscCall(TSGetStepNumber(ts, &n)); 721 PetscCall(TSGetMaxSteps(ts, &maxn)); 722 723 PetscCall(TSReset(ts)); 724 PetscCall(TSSetDM(ts, sw)); 725 /* TODO Check whether TS was set from options */ 726 PetscCall(TSSetFromOptions(ts)); 727 PetscCall(TSSetTime(ts, t)); 728 PetscCall(TSSetMaxTime(ts, maxt)); 729 PetscCall(TSSetTimeStep(ts, dt)); 730 PetscCall(TSSetStepNumber(ts, n)); 731 PetscCall(TSSetMaxSteps(ts, maxn)); 732 733 PetscCall(CreateSolution(ts)); 734 PetscCall(SetProblem(ts)); 735 PetscCall(TSGetSolution(ts, &u)); 736 PetscFunctionReturn(PETSC_SUCCESS); 737 } 738 739 PetscErrorCode circleSingleX(PetscInt dim, PetscReal time, const PetscReal dummy[], PetscInt p, PetscScalar x[], void *ctx) 740 { 741 x[0] = p + 1.; 742 x[1] = 0.; 743 return PETSC_SUCCESS; 744 } 745 746 PetscErrorCode circleSingleV(PetscInt dim, PetscReal time, const PetscReal dummy[], PetscInt p, PetscScalar v[], void *ctx) 747 { 748 v[0] = 0.; 749 v[1] = PetscSqrtReal(1000. / (p + 1.)); 750 return PETSC_SUCCESS; 751 } 752 753 /* Put 5 particles into each circle */ 754 PetscErrorCode circleMultipleX(PetscInt dim, PetscReal time, const PetscReal dummy[], PetscInt p, PetscScalar x[], void *ctx) 755 { 756 const PetscInt n = 5; 757 const PetscReal r0 = (p / n) + 1.; 758 const PetscReal th0 = (2. * PETSC_PI * (p % n)) / n; 759 760 x[0] = r0 * PetscCosReal(th0); 761 x[1] = r0 * PetscSinReal(th0); 762 return PETSC_SUCCESS; 763 } 764 765 /* Put 5 particles into each circle */ 766 PetscErrorCode circleMultipleV(PetscInt dim, PetscReal time, const PetscReal dummy[], PetscInt p, PetscScalar v[], void *ctx) 767 { 768 const PetscInt n = 5; 769 const PetscReal r0 = (p / n) + 1.; 770 const PetscReal th0 = (2. * PETSC_PI * (p % n)) / n; 771 const PetscReal omega = PetscSqrtReal(1000. / r0) / r0; 772 773 v[0] = -r0 * omega * PetscSinReal(th0); 774 v[1] = r0 * omega * PetscCosReal(th0); 775 return PETSC_SUCCESS; 776 } 777 778 /* 779 InitializeSolveAndSwarm - Set the solution values to the swarm coordinates and velocities, and also possibly set the initial values. 780 781 Input Parameters: 782 + ts - The TS 783 - useInitial - Flag to also set the initial conditions to the current coordinates and velocities and setup the problem 784 785 Output Parameter: 786 . u - The initialized solution vector 787 788 Level: advanced 789 790 .seealso: InitializeSolve() 791 */ 792 static PetscErrorCode InitializeSolveAndSwarm(TS ts, PetscBool useInitial) 793 { 794 DM sw; 795 Vec u, gc, gv, gc0, gv0; 796 IS isx, isv; 797 AppCtx *user; 798 799 PetscFunctionBeginUser; 800 PetscCall(TSGetDM(ts, &sw)); 801 PetscCall(DMGetApplicationContext(sw, &user)); 802 if (useInitial) { 803 PetscReal v0[1] = {1.}; 804 805 PetscCall(DMSwarmInitializeCoordinates(sw)); 806 PetscCall(DMSwarmInitializeVelocitiesFromOptions(sw, v0)); 807 PetscCall(DMSwarmMigrate(sw, PETSC_TRUE)); 808 PetscCall(DMSwarmTSRedistribute(ts)); 809 } 810 PetscCall(TSGetSolution(ts, &u)); 811 PetscCall(TSRHSSplitGetIS(ts, "position", &isx)); 812 PetscCall(TSRHSSplitGetIS(ts, "momentum", &isv)); 813 PetscCall(DMSwarmCreateGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc)); 814 PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "initCoordinates", &gc0)); 815 if (useInitial) PetscCall(VecCopy(gc, gc0)); 816 PetscCall(VecISCopy(u, isx, SCATTER_FORWARD, gc)); 817 PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc)); 818 PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "initCoordinates", &gc0)); 819 PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "velocity", &gv)); 820 PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "initVelocity", &gv0)); 821 if (useInitial) PetscCall(VecCopy(gv, gv0)); 822 PetscCall(VecISCopy(u, isv, SCATTER_FORWARD, gv)); 823 PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "velocity", &gv)); 824 PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "initVelocity", &gv0)); 825 PetscFunctionReturn(PETSC_SUCCESS); 826 } 827 828 static PetscErrorCode InitializeSolve(TS ts, Vec u) 829 { 830 PetscFunctionBegin; 831 PetscCall(TSSetSolution(ts, u)); 832 PetscCall(InitializeSolveAndSwarm(ts, PETSC_TRUE)); 833 PetscFunctionReturn(PETSC_SUCCESS); 834 } 835 836 static PetscErrorCode ComputeError(TS ts, Vec U, Vec E) 837 { 838 MPI_Comm comm; 839 DM sw; 840 AppCtx *user; 841 const PetscScalar *u; 842 const PetscReal *coords, *vel; 843 PetscScalar *e; 844 PetscReal t; 845 PetscInt dim, Np, p; 846 847 PetscFunctionBeginUser; 848 PetscCall(PetscObjectGetComm((PetscObject)ts, &comm)); 849 PetscCall(TSGetDM(ts, &sw)); 850 PetscCall(DMGetApplicationContext(sw, &user)); 851 PetscCall(DMGetDimension(sw, &dim)); 852 PetscCall(TSGetSolveTime(ts, &t)); 853 PetscCall(VecGetArray(E, &e)); 854 PetscCall(VecGetArrayRead(U, &u)); 855 PetscCall(VecGetLocalSize(U, &Np)); 856 PetscCall(DMSwarmGetField(sw, "initCoordinates", NULL, NULL, (void **)&coords)); 857 PetscCall(DMSwarmGetField(sw, "initVelocity", NULL, NULL, (void **)&vel)); 858 Np /= 2 * dim; 859 for (p = 0; p < Np; ++p) { 860 /* TODO generalize initial conditions and project into plane instead of assuming x-y */ 861 const PetscReal r0 = DMPlex_NormD_Internal(dim, &coords[p * dim]); 862 const PetscReal th0 = PetscAtan2Real(coords[p * dim + 1], coords[p * dim + 0]); 863 const PetscReal v0 = DMPlex_NormD_Internal(dim, &vel[p * dim]); 864 const PetscReal omega = v0 / r0; 865 const PetscReal ct = PetscCosReal(omega * t + th0); 866 const PetscReal st = PetscSinReal(omega * t + th0); 867 const PetscScalar *x = &u[(p * 2 + 0) * dim]; 868 const PetscScalar *v = &u[(p * 2 + 1) * dim]; 869 const PetscReal xe[3] = {r0 * ct, r0 * st, 0.0}; 870 const PetscReal ve[3] = {-v0 * st, v0 * ct, 0.0}; 871 PetscInt d; 872 873 for (d = 0; d < dim; ++d) { 874 e[(p * 2 + 0) * dim + d] = x[d] - xe[d]; 875 e[(p * 2 + 1) * dim + d] = v[d] - ve[d]; 876 } 877 if (user->error) { 878 const PetscReal en = 0.5 * DMPlex_DotRealD_Internal(dim, v, v); 879 const PetscReal exen = 0.5 * PetscSqr(v0); 880 PetscCall(PetscPrintf(comm, "t %.4g: p%" PetscInt_FMT " error [%.2g %.2g] 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))); 881 } 882 } 883 PetscCall(DMSwarmRestoreField(sw, "initCoordinates", NULL, NULL, (void **)&coords)); 884 PetscCall(DMSwarmRestoreField(sw, "initVelocity", NULL, NULL, (void **)&vel)); 885 PetscCall(VecRestoreArrayRead(U, &u)); 886 PetscCall(VecRestoreArray(E, &e)); 887 PetscFunctionReturn(PETSC_SUCCESS); 888 } 889 890 static PetscErrorCode EnergyMonitor(TS ts, PetscInt step, PetscReal t, Vec U, void *ctx) 891 { 892 const PetscInt ostep = ((AppCtx *)ctx)->ostep; 893 const EMType em = ((AppCtx *)ctx)->em; 894 DM sw; 895 const PetscScalar *u; 896 PetscReal *coords, *E; 897 PetscReal enKin = 0., enEM = 0.; 898 PetscInt dim, d, Np, p, q; 899 900 PetscFunctionBeginUser; 901 if (step % ostep == 0) { 902 PetscCall(TSGetDM(ts, &sw)); 903 PetscCall(DMGetDimension(sw, &dim)); 904 PetscCall(VecGetArrayRead(U, &u)); 905 PetscCall(VecGetLocalSize(U, &Np)); 906 Np /= 2 * dim; 907 PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords)); 908 PetscCall(DMSwarmGetField(sw, "E_field", NULL, NULL, (void **)&E)); 909 if (!step) PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "Time Step Part Energy\n")); 910 for (p = 0; p < Np; ++p) { 911 const PetscReal v2 = DMPlex_DotRealD_Internal(dim, &u[(p * 2 + 1) * dim], &u[(p * 2 + 1) * dim]); 912 PetscReal *pcoord = &coords[p * dim]; 913 914 PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)ts), "%.6lf %4" PetscInt_FMT " %5" PetscInt_FMT " %10.4lf\n", (double)t, step, p, (double)(0.5 * v2))); 915 enKin += 0.5 * v2; 916 if (em == EM_NONE) { 917 continue; 918 } else if (em == EM_COULOMB) { 919 for (q = p + 1; q < Np; ++q) { 920 PetscReal *qcoord = &coords[q * dim]; 921 PetscReal rpq[3], r; 922 for (d = 0; d < dim; ++d) rpq[d] = pcoord[d] - qcoord[d]; 923 r = DMPlex_NormD_Internal(dim, rpq); 924 enEM += 1. / r; 925 } 926 } else if (em == EM_PRIMAL || em == EM_MIXED) { 927 for (d = 0; d < dim; ++d) enEM += E[p * dim + d]; 928 } 929 } 930 PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)ts), "%.6lf %4" PetscInt_FMT " KE\t %10.4lf\n", (double)t, step, (double)enKin)); 931 PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)ts), "%.6lf %4" PetscInt_FMT " PE\t %1.10g\n", (double)t, step, (double)enEM)); 932 PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)ts), "%.6lf %4" PetscInt_FMT " E\t %10.4lf\n", (double)t, step, (double)(enKin + enEM))); 933 PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords)); 934 PetscCall(DMSwarmRestoreField(sw, "E_field", NULL, NULL, (void **)&E)); 935 PetscCall(PetscSynchronizedFlush(PetscObjectComm((PetscObject)ts), NULL)); 936 PetscCall(VecRestoreArrayRead(U, &u)); 937 } 938 PetscFunctionReturn(PETSC_SUCCESS); 939 } 940 941 static PetscErrorCode MigrateParticles(TS ts) 942 { 943 DM sw; 944 945 PetscFunctionBeginUser; 946 PetscCall(TSGetDM(ts, &sw)); 947 PetscCall(DMViewFromOptions(sw, NULL, "-migrate_view_pre")); 948 { 949 Vec u, gc, gv; 950 IS isx, isv; 951 952 PetscCall(TSGetSolution(ts, &u)); 953 PetscCall(TSRHSSplitGetIS(ts, "position", &isx)); 954 PetscCall(TSRHSSplitGetIS(ts, "momentum", &isv)); 955 PetscCall(DMSwarmCreateGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc)); 956 PetscCall(VecISCopy(u, isx, SCATTER_REVERSE, gc)); 957 PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc)); 958 PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "velocity", &gv)); 959 PetscCall(VecISCopy(u, isv, SCATTER_REVERSE, gv)); 960 PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "velocity", &gv)); 961 } 962 PetscCall(DMSwarmMigrate(sw, PETSC_TRUE)); 963 PetscCall(DMSwarmTSRedistribute(ts)); 964 PetscCall(InitializeSolveAndSwarm(ts, PETSC_FALSE)); 965 PetscFunctionReturn(PETSC_SUCCESS); 966 } 967 968 int main(int argc, char **argv) 969 { 970 DM dm, sw; 971 TS ts; 972 Vec u; 973 AppCtx user; 974 975 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 976 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user)); 977 PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm)); 978 PetscCall(CreatePoisson(dm, &user)); 979 PetscCall(CreateSwarm(dm, &user, &sw)); 980 PetscCall(DMSetApplicationContext(sw, &user)); 981 982 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 983 PetscCall(TSSetProblemType(ts, TS_NONLINEAR)); 984 PetscCall(TSSetDM(ts, sw)); 985 PetscCall(TSSetMaxTime(ts, 0.1)); 986 PetscCall(TSSetTimeStep(ts, 0.00001)); 987 PetscCall(TSSetMaxSteps(ts, 100)); 988 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP)); 989 PetscCall(TSMonitorSet(ts, EnergyMonitor, &user, NULL)); 990 PetscCall(TSSetFromOptions(ts)); 991 PetscCall(DMSwarmVectorDefineField(sw, "velocity")); 992 PetscCall(TSSetComputeInitialCondition(ts, InitializeSolve)); 993 PetscCall(TSSetComputeExactError(ts, ComputeError)); 994 PetscCall(TSSetPostStep(ts, MigrateParticles)); 995 996 PetscCall(CreateSolution(ts)); 997 PetscCall(TSGetSolution(ts, &u)); 998 PetscCall(TSComputeInitialCondition(ts, u)); 999 PetscCall(TSSolve(ts, NULL)); 1000 1001 PetscCall(SNESDestroy(&user.snes)); 1002 PetscCall(TSDestroy(&ts)); 1003 PetscCall(DMDestroy(&sw)); 1004 PetscCall(DMDestroy(&dm)); 1005 PetscCall(PetscFinalize()); 1006 return 0; 1007 } 1008 1009 /*TEST 1010 1011 build: 1012 requires: double !complex 1013 1014 testset: 1015 requires: defined(PETSC_HAVE_EXECUTABLE_EXPORT) 1016 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 \ 1017 -dm_swarm_num_particles 2 -dm_swarm_coordinate_function circleSingleX -dm_swarm_velocity_function circleSingleV \ 1018 -ts_type basicsymplectic\ 1019 -dm_view -output_step 50 -ts_dt 0.01 -ts_max_time 10.0 -ts_max_steps 10 1020 test: 1021 suffix: none_bsi_2d_1 1022 args: -ts_basicsymplectic_type 1 -em_type none -error 1023 test: 1024 suffix: none_bsi_2d_2 1025 args: -ts_basicsymplectic_type 2 -em_type none -error 1026 test: 1027 suffix: none_bsi_2d_3 1028 args: -ts_basicsymplectic_type 3 -em_type none -error 1029 test: 1030 suffix: none_bsi_2d_4 1031 args: -ts_basicsymplectic_type 4 -em_type none -error 1032 test: 1033 suffix: coulomb_bsi_2d_1 1034 args: -ts_basicsymplectic_type 1 1035 test: 1036 suffix: coulomb_bsi_2d_2 1037 args: -ts_basicsymplectic_type 2 1038 test: 1039 suffix: coulomb_bsi_2d_3 1040 args: -ts_basicsymplectic_type 3 1041 test: 1042 suffix: coulomb_bsi_2d_4 1043 args: -ts_basicsymplectic_type 4 1044 1045 testset: 1046 requires: defined(PETSC_HAVE_EXECUTABLE_EXPORT) 1047 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 \ 1048 -dm_swarm_num_particles 2 -dm_swarm_coordinate_function circleSingleX -dm_swarm_velocity_function circleSingleV \ 1049 -ts_type basicsymplectic\ 1050 -em_type primal -em_pc_type svd\ 1051 -dm_view -output_step 50 -error -ts_dt 0.01 -ts_max_time 10.0 -ts_max_steps 10\ 1052 -petscspace_degree 2 -petscfe_default_quadrature_order 3 -sigma 1.0e-8 -timeScale 2.0e-14 1053 test: 1054 suffix: poisson_bsi_2d_1 1055 args: -ts_basicsymplectic_type 1 1056 test: 1057 suffix: poisson_bsi_2d_2 1058 args: -ts_basicsymplectic_type 2 1059 test: 1060 suffix: poisson_bsi_2d_3 1061 args: -ts_basicsymplectic_type 3 1062 test: 1063 suffix: poisson_bsi_2d_4 1064 args: -ts_basicsymplectic_type 4 1065 1066 testset: 1067 requires: defined(PETSC_HAVE_EXECUTABLE_EXPORT) 1068 args: -dm_swarm_num_particles 2 -dm_swarm_coordinate_function circleSingleX -dm_swarm_velocity_function circleSingleV \ 1069 -ts_type theta -ts_theta_theta 0.5 -ts_convergence_estimate -convest_num_refine 2 \ 1070 -mat_type baij -ksp_error_if_not_converged -em_pc_type svd\ 1071 -dm_view -output_step 50 -error -ts_dt 0.01 -ts_max_time 10.0 -ts_max_steps 10\ 1072 -pc_type svd -sigma 1.0e-8 -timeScale 2.0e-14 1073 test: 1074 suffix: im_2d_0 1075 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 1076 1077 testset: 1078 requires: defined(PETSC_HAVE_EXECUTABLE_EXPORT) 1079 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 \ 1080 -dm_swarm_num_particles 2 -dm_swarm_coordinate_function circleSingleX -dm_swarm_velocity_function circleSingleV -dm_swarm_num_species 1\ 1081 -ts_type basicsymplectic -ts_convergence_estimate -convest_num_refine 2 \ 1082 -em_snes_type ksponly -em_pc_type svd -em_type primal -petscspace_degree 1\ 1083 -dm_view -output_step 50\ 1084 -pc_type svd -sigma 1.0e-8 -timeScale 2.0e-14 -ts_dt 0.01 -ts_max_time 1.0 -ts_max_steps 10 1085 test: 1086 suffix: bsi_2d_mesh_1 1087 args: -ts_basicsymplectic_type 4 1088 test: 1089 suffix: bsi_2d_mesh_1_par_2 1090 nsize: 2 1091 args: -ts_basicsymplectic_type 4 1092 test: 1093 suffix: bsi_2d_mesh_1_par_3 1094 nsize: 3 1095 args: -ts_basicsymplectic_type 4 1096 test: 1097 suffix: bsi_2d_mesh_1_par_4 1098 nsize: 4 1099 args: -ts_basicsymplectic_type 4 -dm_swarm_num_particles 0,0,2,0 1100 1101 testset: 1102 requires: defined(PETSC_HAVE_EXECUTABLE_EXPORT) 1103 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 \ 1104 -dm_swarm_num_particles 10 -dm_swarm_coordinate_function circleMultipleX -dm_swarm_velocity_function circleMultipleV \ 1105 -ts_convergence_estimate -convest_num_refine 2 \ 1106 -em_pc_type lu\ 1107 -dm_view -output_step 50 -error\ 1108 -pc_type svd -sigma 1.0e-8 -timeScale 2.0e-14 -ts_dt 0.01 -ts_max_time 10.0 -ts_max_steps 10 1109 test: 1110 suffix: bsi_2d_multiple_1 1111 args: -ts_type basicsymplectic -ts_basicsymplectic_type 1 1112 test: 1113 suffix: bsi_2d_multiple_2 1114 args: -ts_type basicsymplectic -ts_basicsymplectic_type 2 1115 test: 1116 suffix: bsi_2d_multiple_3 1117 args: -ts_type basicsymplectic -ts_basicsymplectic_type 3 -ts_dt 0.001 1118 test: 1119 suffix: im_2d_multiple_0 1120 args: -ts_type theta -ts_theta_theta 0.5 \ 1121 -mat_type baij -ksp_error_if_not_converged -em_pc_type lu 1122 1123 testset: 1124 requires: defined(PETSC_HAVE_EXECUTABLE_EXPORT) 1125 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 \ 1126 -dm_swarm_num_particles 2 -dm_swarm_coordinate_function circleSingleX -dm_swarm_velocity_function circleSingleV \ 1127 -em_pc_type fieldsplit -ksp_rtol 1e-10 -em_ksp_type preonly -em_type mixed -em_ksp_error_if_not_converged\ 1128 -dm_view -output_step 50 -error -dm_refine 0\ 1129 -pc_type svd -sigma 1.0e-8 -timeScale 2.0e-14 -ts_dt 0.01 -ts_max_time 10.0 -ts_max_steps 10 1130 test: 1131 suffix: bsi_4_rt_1 1132 args: -ts_type basicsymplectic -ts_basicsymplectic_type 4\ 1133 -pc_fieldsplit_detect_saddle_point\ 1134 -pc_type fieldsplit\ 1135 -pc_fieldsplit_type schur\ 1136 -pc_fieldsplit_schur_precondition full \ 1137 -field_petscspace_degree 2\ 1138 -field_petscfe_default_quadrature_order 1\ 1139 -field_petscspace_type sum \ 1140 -field_petscspace_variables 2 \ 1141 -field_petscspace_components 2 \ 1142 -field_petscspace_sum_spaces 2 \ 1143 -field_petscspace_sum_concatenate true \ 1144 -field_sumcomp_0_petscspace_variables 2 \ 1145 -field_sumcomp_0_petscspace_type tensor \ 1146 -field_sumcomp_0_petscspace_tensor_spaces 2 \ 1147 -field_sumcomp_0_petscspace_tensor_uniform false \ 1148 -field_sumcomp_0_tensorcomp_0_petscspace_degree 1 \ 1149 -field_sumcomp_0_tensorcomp_1_petscspace_degree 0 \ 1150 -field_sumcomp_1_petscspace_variables 2 \ 1151 -field_sumcomp_1_petscspace_type tensor \ 1152 -field_sumcomp_1_petscspace_tensor_spaces 2 \ 1153 -field_sumcomp_1_petscspace_tensor_uniform false \ 1154 -field_sumcomp_1_tensorcomp_0_petscspace_degree 0 \ 1155 -field_sumcomp_1_tensorcomp_1_petscspace_degree 1 \ 1156 -field_petscdualspace_form_degree -1 \ 1157 -field_petscdualspace_order 1 \ 1158 -field_petscdualspace_lagrange_trimmed true\ 1159 -ksp_gmres_restart 500 1160 1161 TEST*/ 1162