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