1 static char help[] = "Landau collision operator driver\n\n"; 2 3 #include <petscts.h> 4 #include <petsclandau.h> 5 6 typedef struct { 7 PetscInt grid_target, batch_target, field_target; 8 PetscBool init; 9 } ex1Ctx; 10 11 /* 12 call back method for DMPlexLandauAddToFunction: 13 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 20 21 Input/Output Parameters: 22 + x - Vector to data to 23 24 */ 25 PetscErrorCode landau_field_add_to_callback(DM dm, Vec x, PetscInt local_field, PetscInt grid, PetscInt b_id, void *vctx) { 26 ex1Ctx *user = (ex1Ctx *)vctx; 27 28 PetscFunctionBegin; 29 if (grid == user->grid_target && b_id == user->batch_target && local_field == user->field_target) { 30 PetscScalar one = 1.0; 31 PetscCall(VecSet(x, one)); 32 if (!user->init) { 33 PetscCall(PetscObjectSetName((PetscObject)dm, "single")); 34 PetscCall(DMViewFromOptions(dm, NULL, "-ex1_dm_view")); 35 user->init = PETSC_TRUE; 36 } 37 PetscCall(PetscObjectSetName((PetscObject)x, "u")); 38 PetscCall(VecViewFromOptions(x, NULL, "-ex1_vec_view")); // this causes diffs with Kokkos, etc 39 PetscCall(PetscInfo(dm, "DMPlexLandauAddToFunction user 'add' method to grid %" PetscInt_FMT ", batch %" PetscInt_FMT " and local field %" PetscInt_FMT "\n", grid, b_id, local_field)); 40 } 41 PetscFunctionReturn(0); 42 } 43 44 int main(int argc, char **argv) { 45 DM pack; 46 Vec X, X_0; 47 PetscInt dim = 2; 48 TS ts; 49 Mat J; 50 SNES snes; 51 KSP ksp; 52 PC pc; 53 SNESLineSearch linesearch; 54 PetscReal time; 55 56 PetscFunctionBeginUser; 57 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 58 PetscCall(PetscOptionsGetInt(NULL, NULL, "-dim", &dim, NULL)); 59 /* Create a mesh */ 60 PetscCall(DMPlexLandauCreateVelocitySpace(PETSC_COMM_SELF, dim, "", &X, &J, &pack)); 61 PetscCall(DMSetUp(pack)); 62 PetscCall(VecDuplicate(X, &X_0)); 63 PetscCall(VecCopy(X, X_0)); 64 PetscCall(DMPlexLandauPrintNorms(X, 0)); 65 PetscCall(DMSetOutputSequenceNumber(pack, 0, 0.0)); 66 /* Create timestepping solver context */ 67 PetscCall(TSCreate(PETSC_COMM_SELF, &ts)); 68 PetscCall(TSSetDM(ts, pack)); 69 PetscCall(TSGetSNES(ts, &snes)); 70 PetscCall(SNESGetLineSearch(snes, &linesearch)); 71 PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHBASIC)); 72 PetscCall(TSSetIFunction(ts, NULL, DMPlexLandauIFunction, NULL)); 73 PetscCall(TSSetIJacobian(ts, J, J, DMPlexLandauIJacobian, NULL)); 74 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); 75 PetscCall(SNESGetKSP(snes, &ksp)); 76 PetscCall(KSPGetPC(ksp, &pc)); 77 PetscCall(TSSetFromOptions(ts)); 78 PetscCall(TSSetSolution(ts, X)); 79 PetscCall(TSSolve(ts, X)); 80 PetscCall(DMPlexLandauPrintNorms(X, 1)); 81 PetscCall(TSGetTime(ts, &time)); 82 PetscCall(DMSetOutputSequenceNumber(pack, 1, time)); 83 PetscCall(VecAXPY(X, -1, X_0)); 84 { /* test add field method */ 85 ex1Ctx *user; 86 PetscCall(PetscNew(&user)); 87 user->grid_target = 1; // 2nd ion species 88 user->field_target = 1; 89 PetscCall(DMPlexLandauAddToFunction(pack, X, landau_field_add_to_callback, user)); 90 PetscCall(PetscFree(user)); 91 } 92 /* clean up */ 93 PetscCall(DMPlexLandauDestroyVelocitySpace(&pack)); 94 PetscCall(TSDestroy(&ts)); 95 PetscCall(VecDestroy(&X)); 96 PetscCall(VecDestroy(&X_0)); 97 PetscCall(PetscFinalize()); 98 return 0; 99 } 100 101 /*TEST 102 testset: 103 requires: p4est !complex double defined(PETSC_USE_DMLANDAU_2D) 104 output_file: output/ex1_0.out 105 filter: grep -v "% type: seq" 106 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 107 test: 108 suffix: cpu 109 args: -dm_landau_device_type cpu 110 test: 111 suffix: kokkos 112 requires: kokkos_kernels 113 args: -dm_landau_device_type kokkos -dm_mat_type aijkokkos -dm_vec_type kokkos 114 test: 115 suffix: cuda 116 requires: cuda 117 args: -dm_landau_device_type cuda -dm_mat_type aijcusparse -dm_vec_type cuda -mat_cusparse_use_cpu_solve 118 119 TEST*/ 120