1: ! test phase space (Maxwellian) mesh construction (serial)
  2: !
  3: !
  4: !
  5: ! Contributed by Mark Adams
  6: program DMPlexTestLandauInterface
  7:   use petscts
  8:   use petscdmplex
  9: #include <petsc/finclude/petscts.h>
 10: #include <petsc/finclude/petscdmplex.h>
 11:   implicit none
 12:   external DMPlexLandauIFunction
 13:   external DMPlexLandauIJacobian
 14:   DM             dm
 15:   PetscInt       dim
 16:   PetscInt       ii
 17:   PetscErrorCode ierr
 18:   TS             ts
 19:   Vec            X,X_0
 20:   Mat            J
 21:   SNES           snes
 22:   KSP            ksp
 23:   PC             pc
 24:   SNESLineSearch linesearch
 25:   PetscReal      mone
 26:   PetscScalar    scalar
 27:   call PetscInitialize(PETSC_NULL_CHARACTER, ierr)
 28:   if (ierr .ne. 0) then
 29:      print*,'Unable to initialize PETSc'
 30:      stop
 31:   endif
 32:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 33:   !  Create mesh (DM), read in parameters, create and add f_0 (X)
 34:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 35:   dim = 2
 36:   call DMPlexLandauCreateVelocitySpace(PETSC_COMM_WORLD, dim, '', X, J, dm, ierr);CHKERRA(ierr)
 37:   call DMSetUp(dm,ierr);CHKERRA(ierr)
 38:   call VecDuplicate(X,X_0,ierr);CHKERRA(ierr)
 39:   call VecCopy(X,X_0,ierr)
 40:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 41:   !  View
 42:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 43:   ii = 0
 44:   call DMPlexLandauPrintNorms(X,ii,ierr);CHKERRA(ierr)
 45:   mone = 0;
 46:   call DMSetOutputSequenceNumber(dm, ii, mone, ierr);CHKERRA(ierr);
 47:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 48:   !    Create timestepping solver context
 49:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 50:   call TSCreate(PETSC_COMM_SELF,ts,ierr);CHKERRA(ierr)
 51:   call TSSetOptionsPrefix(ts, 'ex1_', ierr);CHKERRA(ierr) ! should get this from the dm or give it to the dm
 52:   call TSSetDM(ts,dm,ierr);CHKERRA(ierr)
 53:   call TSGetSNES(ts,snes,ierr);CHKERRA(ierr)
 54:   call SNESSetOptionsPrefix(snes, 'ex1_', ierr);CHKERRA(ierr) ! should get this from the dm or give it to the dm
 55:   call SNESGetLineSearch(snes,linesearch,ierr);CHKERRA(ierr)
 56:   call SNESLineSearchSetType(linesearch,SNESLINESEARCHBASIC,ierr);CHKERRA(ierr)
 57:   call TSSetIFunction(ts,PETSC_NULL_VEC,DMPlexLandauIFunction,PETSC_NULL_VEC,ierr);CHKERRA(ierr)
 58:   call TSSetIJacobian(ts,J,J,DMPlexLandauIJacobian,PETSC_NULL_VEC,ierr);CHKERRA(ierr)
 59:   call TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER,ierr);CHKERRA(ierr)

 61:   call SNESGetKSP(snes,ksp,ierr);CHKERRA(ierr)
 62:   call KSPSetOptionsPrefix(ksp, 'ex1_', ierr);CHKERRA(ierr) ! should get this from the dm or give it to the dm
 63:   call KSPGetPC(ksp,pc,ierr);CHKERRA(ierr)
 64:   call PCSetOptionsPrefix(pc, 'ex1_', ierr);CHKERRA(ierr) ! should get this from the dm or give it to the dm

 66:   call TSSetFromOptions(ts,ierr);CHKERRA(ierr)
 67:   call TSSetSolution(ts,X,ierr);CHKERRA(ierr)
 68:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 69:   !  Solve nonlinear system
 70:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 71:   call TSSolve(ts,X,ierr);CHKERRA(ierr)
 72:   ii = 1
 73:   call DMPlexLandauPrintNorms(X,ii,ierr);CHKERRA(ierr)
 74:   call TSGetTime(ts, mone, ierr);CHKERRA(ierr);
 75:   call DMSetOutputSequenceNumber(dm, ii, mone, ierr);CHKERRA(ierr);
 76:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 77:   !  remove f_0
 78:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 79:   scalar = -1.
 80:   call VecAXPY(X,scalar,X_0,ierr);CHKERRA(ierr)
 81:   call DMPlexLandauDestroyVelocitySpace(dm, ierr);CHKERRA(ierr)
 82:   call TSDestroy(ts, ierr);CHKERRA(ierr)
 83:   call VecDestroy(X, ierr);CHKERRA(ierr)
 84:   call VecDestroy(X_0, ierr);CHKERRA(ierr)
 85:   call PetscFinalize(ierr)
 86: end program DMPlexTestLandauInterface

 88: !/*TEST
 89: !  build:
 90: !    requires: defined(PETSC_USING_F90FREEFORM)
 91: !
 92: !  test:
 93: !    suffix: 0
 94: !    requires: p4est !complex  !kokkos_kernels !cuda
 95: !    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 -1 -ex1_ts_rtol 1e-1 -ex1_ts_dt 1.e-1 -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
 96: !
 97: !TEST*/