Actual source code: ex1.c

  1: static char help[] = "Landau collision operator driver\n\n";

  3: #include <petscts.h>
  4: #include <petsclandau.h>

  6: typedef struct {
  7:   PetscInt  grid_target, batch_target, field_target;
  8:   PetscBool init;
  9: } ex1Ctx;

 11: /*
 12:  call back method for DMPlexLandauAccess:

 14: Input Parameters:
 15:  .   dm - a DM for this field
 16:  -   local_field - the local index in the grid for this field
 17:  .   grid - the grid index
 18:  +   b_id - the batch index
 19:  -   vctx - a user context

 21:  Input/Output Parameters:
 22:  +   x - Vector to data to

 24:  */
 25: PetscErrorCode landau_field_print_access_callback(DM dm, Vec x, PetscInt local_field, PetscInt grid, PetscInt b_id, void *vctx)
 26: {
 27:   ex1Ctx    *user = (ex1Ctx *)vctx;
 28:   LandauCtx *ctx;

 30:   DMGetApplicationContext(dm, &ctx);
 31:   if (grid == user->grid_target && b_id == user->batch_target && local_field == user->field_target) {
 32:     PetscScalar one = 1.0e10;

 34:     VecSet(x, one);
 35:     if (!user->init) {
 36:       PetscObjectSetName((PetscObject)dm, "single");
 37:       DMViewFromOptions(dm, NULL, "-ex1_dm_view"); // DMCreateSubDM does seem to give the DM's fild the name from the original DM
 38:       user->init = PETSC_TRUE;
 39:     }
 40:     PetscObjectSetName((PetscObject)x, "u");      // this gives the vector a nicer name, DMCreateSubDM could do this for us and get the correct name
 41:     VecViewFromOptions(x, NULL, "-ex1_vec_view"); // this causes diffs with Kokkos, etc
 42:     PetscInfo(dm, "DMPlexLandauAccess user 'add' method to grid %" PetscInt_FMT ", batch %" PetscInt_FMT " and species %" PetscInt_FMT "\n", grid, b_id, ctx->species_offset[grid] + local_field);
 43:   }
 44:   return 0;
 45: }

 47: int main(int argc, char **argv)
 48: {
 49:   DM             pack;
 50:   Vec            X, X_0;
 51:   PetscInt       dim = 2;
 52:   TS             ts;
 53:   Mat            J;
 54:   SNES           snes;
 55:   KSP            ksp;
 56:   PC             pc;
 57:   SNESLineSearch linesearch;
 58:   PetscReal      time;

 61:   PetscInitialize(&argc, &argv, NULL, help);
 62:   PetscOptionsGetInt(NULL, NULL, "-dim", &dim, NULL);
 63:   /* Create a mesh */
 64:   DMPlexLandauCreateVelocitySpace(PETSC_COMM_SELF, dim, "", &X, &J, &pack);
 65:   DMSetUp(pack);
 66:   VecDuplicate(X, &X_0);
 67:   VecCopy(X, X_0);
 68:   DMPlexLandauPrintNorms(X, 0);
 69:   DMSetOutputSequenceNumber(pack, 0, 0.0);
 70:   /* Create timestepping solver context */
 71:   TSCreate(PETSC_COMM_SELF, &ts);
 72:   TSSetDM(ts, pack);
 73:   TSGetSNES(ts, &snes);
 74:   SNESGetLineSearch(snes, &linesearch);
 75:   SNESLineSearchSetType(linesearch, SNESLINESEARCHBASIC);
 76:   TSSetIFunction(ts, NULL, DMPlexLandauIFunction, NULL);
 77:   TSSetIJacobian(ts, J, J, DMPlexLandauIJacobian, NULL);
 78:   TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
 79:   SNESGetKSP(snes, &ksp);
 80:   KSPGetPC(ksp, &pc);
 81:   TSSetFromOptions(ts);
 82:   TSSetSolution(ts, X);
 83:   TSSolve(ts, X);
 84:   DMPlexLandauPrintNorms(X, 1);
 85:   TSGetTime(ts, &time);
 86:   DMSetOutputSequenceNumber(pack, 1, time);
 87:   VecAXPY(X, -1, X_0);
 88:   { /* test add field method */
 89:     ex1Ctx    *user;
 90:     LandauCtx *ctx;
 91:     DMGetApplicationContext(pack, &ctx);
 92:     PetscNew(&user);
 93:     user->grid_target  = 1; // 2nd ion species
 94:     user->field_target = 1;
 95:     DMPlexLandauAccess(pack, X, landau_field_print_access_callback, user);
 96:     PetscFree(user);
 97:   }
 98:   /* clean up */
 99:   DMPlexLandauDestroyVelocitySpace(&pack);
100:   TSDestroy(&ts);
101:   VecDestroy(&X);
102:   VecDestroy(&X_0);
103:   PetscFinalize();
104:   return 0;
105: }

107: /*TEST
108:   testset:
109:     requires: p4est !complex double defined(PETSC_USE_DMLANDAU_2D)
110:     output_file: output/ex1_0.out
111:     filter: grep -v "%  type: seq"
112:     args: -dm_landau_num_species_grid 1,2 -petscspace_degree 3 -petscspace_poly_tensor 1 -dm_landau_type p4est -dm_landau_ion_masses 2,4 -dm_landau_ion_charges 1,18 -dm_landau_thermal_temps 5,5,.5 -dm_landau_n 1.00018,1,1e-5 -dm_landau_n_0 1e20 -ts_monitor -snes_rtol 1.e-14 -snes_stol 1.e-14 -snes_monitor -snes_converged_reason -ts_type arkimex -ts_arkimex_type 1bee -ts_max_snes_failures -1 -ts_rtol 1e-1 -ts_dt 1.e-1 -ts_max_time 1 -ts_adapt_clip .5,1.25 -ts_adapt_scale_solve_failed 0.75 -ts_adapt_time_step_increase_delay 5 -ts_max_steps 1 -pc_type lu -ksp_type preonly -dm_landau_amr_levels_max 2,1 -ex1_dm_view  -ex1_vec_view ::ascii_matlab
113:     test:
114:       suffix: cpu
115:       args: -dm_landau_device_type cpu
116:     test:
117:       suffix: kokkos
118:       requires: kokkos_kernels
119:       args: -dm_landau_device_type kokkos -dm_mat_type aijkokkos -dm_vec_type kokkos
120:     test:
121:       suffix: cuda
122:       requires: cuda
123:       args: -dm_landau_device_type cuda -dm_mat_type aijcusparse -dm_vec_type cuda -mat_cusparse_use_cpu_solve

125: TEST*/