1: ! test phase space (Maxwellian) mesh construction (serial)
  2: !
  3: ! Contributed by Mark Adams
  4: #include <petsc/finclude/petscts.h>
  5: #include <petsc/finclude/petscdmplex.h>
  6: program DMPlexTestLandauInterface
  7:   use petscts
  8:   use petscdmplex
  9:   implicit none

 11:   external DMPlexLandauIFunction
 12:   external DMPlexLandauIJacobian
 13:   DM dm
 14:   PetscInt, parameter :: dim = 2
 15:   PetscErrorCode ierr
 16:   TS ts
 17:   Vec X, X_0
 18:   Mat J
 19:   SNES snes
 20:   KSP ksp
 21:   PC pc
 22:   SNESLineSearch linesearch
 23:   PetscReal :: time
 24:   PetscScalar, parameter :: scalar = -1.0

 26:   PetscCallA(PetscInitialize(ierr))

 28:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 29:   !  Create mesh (DM), read in parameters, create and add f_0 (X)
 30:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 31:   PetscCallA(DMPlexLandauCreateVelocitySpace(PETSC_COMM_WORLD, dim, '', X, J, dm, ierr))
 32:   PetscCallA(DMSetUp(dm, ierr))
 33:   PetscCallA(VecDuplicate(X, X_0, ierr))
 34:   PetscCallA(VecCopy(X, X_0, ierr))
 35:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 36:   !  View
 37:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 38:   PetscCallA(DMPlexLandauPrintNorms(X, 0_PETSC_INT_KIND, ierr))
 39:   PetscCallA(DMSetOutputSequenceNumber(dm, 0_PETSC_INT_KIND, 0.0_PETSC_REAL_KIND, ierr))
 40:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 41:   !    Create timestepping solver context
 42:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 43:   PetscCallA(TSCreate(PETSC_COMM_SELF, ts, ierr))
 44:   PetscCallA(TSSetOptionsPrefix(ts, 'ex1_', ierr)) ! should get this from the dm or give it to the dm
 45:   PetscCallA(TSSetDM(ts, dm, ierr))
 46:   PetscCallA(TSGetSNES(ts, snes, ierr))
 47:   PetscCallA(SNESSetOptionsPrefix(snes, 'ex1_', ierr)) ! should get this from the dm or give it to the dm
 48:   PetscCallA(SNESGetLineSearch(snes, linesearch, ierr))
 49:   PetscCallA(SNESLineSearchSetType(linesearch, SNESLINESEARCHBASIC, ierr))
 50:   PetscCallA(TSSetIFunction(ts, PETSC_NULL_VEC, DMPlexLandauIFunction, PETSC_NULL_VEC, ierr))
 51:   PetscCallA(TSSetIJacobian(ts, J, J, DMPlexLandauIJacobian, PETSC_NULL_VEC, ierr))
 52:   PetscCallA(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER, ierr))

 54:   PetscCallA(SNESGetKSP(snes, ksp, ierr))
 55:   PetscCallA(KSPSetOptionsPrefix(ksp, 'ex1_', ierr)) ! should get this from the dm or give it to the dm
 56:   PetscCallA(KSPGetPC(ksp, pc, ierr))
 57:   PetscCallA(PCSetOptionsPrefix(pc, 'ex1_', ierr)) ! should get this from the dm or give it to the dm

 59:   PetscCallA(TSSetFromOptions(ts, ierr))
 60:   PetscCallA(TSSetSolution(ts, X, ierr))
 61:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 62:   !  Solve nonlinear system
 63:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 64:   PetscCallA(TSSolve(ts, X, ierr))
 65:   PetscCallA(DMPlexLandauPrintNorms(X, 1_PETSC_INT_KIND, ierr))
 66:   PetscCallA(TSGetTime(ts, time, ierr))
 67:   PetscCallA(DMSetOutputSequenceNumber(dm, 1_PETSC_INT_KIND, time, ierr))
 68:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 69:   !  remove f_0
 70:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 71:   PetscCallA(VecAXPY(X, scalar, X_0, ierr))
 72:   PetscCallA(DMPlexLandauDestroyVelocitySpace(dm, ierr))
 73:   PetscCallA(TSDestroy(ts, ierr))
 74:   PetscCallA(VecDestroy(X, ierr))
 75:   PetscCallA(VecDestroy(X_0, ierr))
 76:   PetscCallA(PetscFinalize(ierr))
 77: end program DMPlexTestLandauInterface

 79: !/*TEST
 80: !
 81: !  build:
 82: !    requires: defined(PETSC_USE_DMLANDAU_2D)
 83: !
 84: !  test:
 85: !    suffix: 0
 86: !    requires: p4est !complex !kokkos_kernels !cuda
 87: !    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_time_step 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
 88: !
 89: !TEST*/