1! test phase space (Maxwellian) mesh construction (serial) 2! 3! Contributed by Mark Adams 4program DMPlexTestLandauInterface 5#include <petsc/finclude/petscts.h> 6#include <petsc/finclude/petscdmplex.h> 7 use petscts 8 use petscdmplex 9 implicit none 10 11 external DMPlexLandauIFunction 12 external DMPlexLandauIJacobian 13 DM dm 14 PetscInt dim 15 PetscInt ii 16 PetscErrorCode ierr 17 TS ts 18 Vec X, X_0 19 Mat J 20 SNES snes 21 KSP ksp 22 PC pc 23 SNESLineSearch linesearch 24 PetscReal mone 25 PetscScalar scalar 26 27 PetscCallA(PetscInitialize(ierr)) 28 29 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 30 ! Create mesh (DM), read in parameters, create and add f_0 (X) 31 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 32 dim = 2 33 PetscCallA(DMPlexLandauCreateVelocitySpace(PETSC_COMM_WORLD, dim, '', X, J, dm, ierr)) 34 PetscCallA(DMSetUp(dm, ierr)) 35 PetscCallA(VecDuplicate(X, X_0, ierr)) 36 PetscCallA(VecCopy(X, X_0, ierr)) 37 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 38 ! View 39 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 40 ii = 0 41 PetscCallA(DMPlexLandauPrintNorms(X, ii, ierr)) 42 mone = 0 43 PetscCallA(DMSetOutputSequenceNumber(dm, ii, mone, ierr)) 44 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 45 ! Create timestepping solver context 46 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 47 PetscCallA(TSCreate(PETSC_COMM_SELF, ts, ierr)) 48 PetscCallA(TSSetOptionsPrefix(ts, 'ex1_', ierr)) ! should get this from the dm or give it to the dm 49 PetscCallA(TSSetDM(ts, dm, ierr)) 50 PetscCallA(TSGetSNES(ts, snes, ierr)) 51 PetscCallA(SNESSetOptionsPrefix(snes, 'ex1_', ierr)) ! should get this from the dm or give it to the dm 52 PetscCallA(SNESGetLineSearch(snes, linesearch, ierr)) 53 PetscCallA(SNESLineSearchSetType(linesearch, SNESLINESEARCHBASIC, ierr)) 54 PetscCallA(TSSetIFunction(ts, PETSC_NULL_VEC, DMPlexLandauIFunction, PETSC_NULL_VEC, ierr)) 55 PetscCallA(TSSetIJacobian(ts, J, J, DMPlexLandauIJacobian, PETSC_NULL_VEC, ierr)) 56 PetscCallA(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER, ierr)) 57 58 PetscCallA(SNESGetKSP(snes, ksp, ierr)) 59 PetscCallA(KSPSetOptionsPrefix(ksp, 'ex1_', ierr)) ! should get this from the dm or give it to the dm 60 PetscCallA(KSPGetPC(ksp, pc, ierr)) 61 PetscCallA(PCSetOptionsPrefix(pc, 'ex1_', ierr)) ! should get this from the dm or give it to the dm 62 63 PetscCallA(TSSetFromOptions(ts, ierr)) 64 PetscCallA(TSSetSolution(ts, X, ierr)) 65 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 66 ! Solve nonlinear system 67 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 68 PetscCallA(TSSolve(ts, X, ierr)) 69 ii = 1 70 PetscCallA(DMPlexLandauPrintNorms(X, ii, ierr)) 71 PetscCallA(TSGetTime(ts, mone, ierr)) 72 PetscCallA(DMSetOutputSequenceNumber(dm, ii, mone, ierr)) 73 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 74 ! remove f_0 75 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 76 scalar = -1. 77 PetscCallA(VecAXPY(X, scalar, X_0, ierr)) 78 PetscCallA(DMPlexLandauDestroyVelocitySpace(dm, ierr)) 79 PetscCallA(TSDestroy(ts, ierr)) 80 PetscCallA(VecDestroy(X, ierr)) 81 PetscCallA(VecDestroy(X_0, ierr)) 82 PetscCallA(PetscFinalize(ierr)) 83end program DMPlexTestLandauInterface 84 85!/*TEST 86! build: 87! requires: defined(PETSC_USING_F90FREEFORM) defined(PETSC_USE_DMLANDAU_2D) 88! 89! test: 90! suffix: 0 91! requires: p4est !complex !kokkos_kernels !cuda 92! 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 -ex1_ts_monitor -ex1_snes_rtol 1.e-14 -ex1_snes_stol 1.e-14 -ex1_snes_monitor -ex1_snes_converged_reason -ex1_ts_type arkimex -ex1_ts_arkimex_type 1bee -ex1_ts_max_snes_failures unlimited -ex1_ts_rtol 1e-1 -ex1_ts_dt 1.e-3 -ex1_ts_max_time 1 -ex1_ts_adapt_clip .5,1.25 -ex1_ts_adapt_scale_solve_failed 0.75 -ex1_ts_adapt_time_step_increase_delay 5 -ex1_ts_max_steps 1 -ex1_pc_type lu -ex1_ksp_type preonly -dm_landau_amr_levels_max 2,1 -dm_landau_device_type cpu -dm_landau_verbose 4 -dm_landau_normalization_grid 1 93! 94!TEST*/ 95