Actual source code: ex1f90.F90

      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

     28:   PetscCallA(PetscInitialize(ierr))

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

     59:   PetscCallA(SNESGetKSP(snes,ksp,ierr))
     60:   PetscCallA(KSPSetOptionsPrefix(ksp, 'ex1_', ierr)) ! should get this from the dm or give it to the dm
     61:   PetscCallA(KSPGetPC(ksp,pc,ierr))
     62:   PetscCallA(PCSetOptionsPrefix(pc, 'ex1_', ierr)) ! should get this from the dm or give it to the dm

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

     86: !/*TEST
     87: !  build:
     88: !    requires: defined(PETSC_USING_F90FREEFORM) defined(PETSC_USE_DMLANDAU_2D)
     89: !
     90: !  test:
     91: !    suffix: 0
     92: !    requires: p4est !complex !kokkos_kernels !cuda
     93: !    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
     94: !
     95: !TEST*/