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*/