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 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 108 test: 109 suffix: cpu 110 args: -dm_landau_device_type cpu 111 test: 112 suffix: kokkos 113 requires: kokkos_kernels 114 args: -dm_landau_device_type kokkos -dm_mat_type aijkokkos -dm_vec_type kokkos 115 test: 116 suffix: cuda 117 requires: cuda 118 args: -dm_landau_device_type cuda -dm_mat_type aijcusparse -dm_vec_type cuda -mat_cusparse_use_cpu_solve 119 120 TEST*/ 121