Actual source code: ex1f90.F90

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

     65:   call SNESGetKSP(snes,ksp,ierr);CHKERRA(ierr)
     66:   call KSPSetOptionsPrefix(ksp, 'ex1f90_', ierr);CHKERRA(ierr) ! should get this from the dm or give it to the dm
     67:   call KSPGetPC(ksp,pc,ierr);CHKERRA(ierr)
     68:   call PCSetOptionsPrefix(pc, 'ex1f90_', ierr);CHKERRA(ierr) ! should get this from the dm or give it to the dm

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

     92: !/*TEST
     93: !  build:
     94: !    requires: define(PETSC_USING_F90FREEFORM)
     95: !
     96: !  test:
     97: !    suffix: 0
     98: !    requires: p4est !complex
     99: !    args: -petscspace_degree 3 -petscspace_poly_tensor 1 -dm_landau_type p4est -info :dm,tsadapt -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 -ex1f90_ts_monitor -ex1f90_snes_rtol 1.e-6 -ex1f90_snes_monitor -ex1f90_snes_converged_reason -ex1f90_ts_type arkimex -ex1f90_ts_arkimex_type 1bee -ex1f90_ts_max_snes_failures -1 -ex1f90_ts_rtol 1e-4 -ex1f90_ts_dt 1.e-6 -ex1f90_ts_max_time 1 -ex1f90_ts_adapt_clip .5,1.25 -ex1f90_ts_adapt_scale_solve_failed 0.75 -ex1f90_ts_adapt_time_step_increase_delay 5 -ex1f90_ts_max_steps 1 -ex1f90_pc_type lu -ex1f90_ksp_type preonly -dm_landau_amr_levels_max 7 -dm_landau_domain_radius 5 -dm_landau_amr_re_levels 0 -dm_landau_re_radius 1 -dm_landau_amr_z_refine1 1 -dm_landau_amr_z_refine2 0 -dm_landau_amr_post_refine 0 -dm_landau_z_radius1 .1 -dm_landau_z_radius2 .1 -dm_refine 1 -dm_landau_device_type cpu
    100: !
    101: !TEST*/