1 #if !defined(PETSCLANDAU_H) 2 #define PETSCLANDAU_H 3 4 #include <petscdmplex.h> /*I "petscdmplex.h" I*/ 5 #include <petscts.h> 6 7 PETSC_EXTERN PetscErrorCode LandauPrintNorms(Vec, PetscInt); 8 PETSC_EXTERN PetscErrorCode LandauCreateVelocitySpace(MPI_Comm,PetscInt,const char[],Vec*,Mat*,DM*); 9 PETSC_EXTERN PetscErrorCode LandauDestroyVelocitySpace(DM*); 10 PETSC_EXTERN PetscErrorCode LandauAddMaxwellians(DM, Vec, PetscReal, PetscReal[], PetscReal[], void *); 11 PETSC_EXTERN PetscErrorCode LandauCreateMassMatrix(DM dm, Mat *Amat); 12 PETSC_EXTERN PetscErrorCode LandauIFunction(TS, PetscReal,Vec,Vec,Vec,void *); 13 PETSC_EXTERN PetscErrorCode LandauIJacobian(TS, PetscReal,Vec,Vec,PetscReal,Mat,Mat,void *); 14 15 /* the Fokker-Planck-Landau context */ 16 #if !defined(LANDAU_DIM) 17 #define LANDAU_DIM 2 18 #endif 19 20 #if !defined(LANDAU_MAX_SPECIES) 21 #if LANDAU_DIM==2 22 #define LANDAU_MAX_SPECIES 10 23 #else 24 #define LANDAU_MAX_SPECIES 3 25 #endif 26 #endif 27 28 #if !defined(LANDAU_MAX_Q) 29 #if defined(LANDAU_MAX_NQ) 30 #error"LANDAU_MAX_NQ but not LANDAU_MAX_Q. Use -DLANDAU_MAX_Q=4 for Q3 elements" 31 #endif 32 #if LANDAU_DIM==2 33 #define LANDAU_MAX_Q 5 34 #define LANDAU_MAX_Q_FACE LANDAU_MAX_Q 35 #else 36 #define LANDAU_MAX_Q 3 37 #define LANDAU_MAX_Q_FACE (LANDAU_MAX_Q*LANDAU_MAX_Q) 38 #endif 39 #else 40 #undef LANDAU_MAX_NQ 41 #if LANDAU_DIM==2 42 #define LANDAU_MAX_NQ (LANDAU_MAX_Q*LANDAU_MAX_Q) 43 #define LANDAU_MAX_Q_FACE LANDAU_MAX_Q 44 #else 45 #define LANDAU_MAX_NQ (LANDAU_MAX_Q*LANDAU_MAX_Q*LANDAU_MAX_Q) 46 #define LANDAU_MAX_Q_FACE (LANDAU_MAX_Q*LANDAU_MAX_Q) 47 #endif 48 #endif 49 50 #if LANDAU_DIM==2 51 #define LANDAU_MAX_NQ (LANDAU_MAX_Q*LANDAU_MAX_Q) 52 #else 53 #define LANDAU_MAX_NQ (LANDAU_MAX_Q*LANDAU_MAX_Q*LANDAU_MAX_Q) 54 #endif 55 56 typedef enum {LANDAU_CUDA, LANDAU_KOKKOS, LANDAU_CPU} LandauDeviceType; 57 typedef struct { 58 PetscBool interpolate; /* Generate intermediate mesh elements */ 59 PetscBool gpu_assembly; 60 PetscFE fe[LANDAU_MAX_SPECIES]; 61 /* geometry */ 62 PetscReal i_radius; 63 PetscReal e_radius; 64 PetscInt num_sections; 65 PetscReal radius; 66 PetscReal re_radius; /* radius of refinement along v_perp=0, z>0 */ 67 PetscReal vperp0_radius1; /* radius of refinement along v_perp=0 */ 68 PetscReal vperp0_radius2; /* radius of refinement along v_perp=0 after origin AMR refinement */ 69 PetscBool sphere; 70 PetscBool inflate; 71 PetscInt numRERefine; /* refinement along v_perp=0, z > 0 */ 72 PetscInt nZRefine1; /* origin refinement after v_perp=0 refinement */ 73 PetscInt nZRefine2; /* origin refinement after origin AMR refinement */ 74 PetscInt maxRefIts; /* normal AMR - refine from origin */ 75 PetscInt postAMRRefine; /* uniform refinement of AMR */ 76 /* discretization - AMR */ 77 PetscErrorCode (*errorIndicator)(PetscInt, PetscReal, PetscReal [], PetscInt, const PetscInt[], const PetscScalar[], const PetscScalar[], PetscReal *, void *); 78 PetscReal refineTol[LANDAU_MAX_SPECIES]; 79 PetscReal coarsenTol[LANDAU_MAX_SPECIES]; 80 /* physics */ 81 PetscReal thermal_temps[LANDAU_MAX_SPECIES]; 82 PetscReal masses[LANDAU_MAX_SPECIES]; /* mass of each species */ 83 PetscReal charges[LANDAU_MAX_SPECIES]; /* charge of each species */ 84 PetscReal n[LANDAU_MAX_SPECIES]; /* number density of each species */ 85 PetscReal m_0; /* reference mass */ 86 PetscReal v_0; /* reference velocity */ 87 PetscReal n_0; /* reference number density */ 88 PetscReal t_0; /* reference time */ 89 PetscReal Ez; 90 PetscReal epsilon0; 91 PetscReal k; 92 PetscReal lnLam; 93 PetscReal electronShift; /* for tests */ 94 PetscInt num_species; 95 /* diagnostics */ 96 PetscInt verbose; 97 PetscLogEvent events[20]; 98 DM dmv; 99 /* cache */ 100 Mat J; 101 Mat M; 102 Vec X; 103 PetscReal normJ; /* used to see if function changed */ 104 /* derived type */ 105 void *data; 106 PetscBool aux_bool; /* helper */ 107 /* computing */ 108 LandauDeviceType deviceType; 109 PetscInt subThreadBlockSize; 110 } LandauCtx; 111 112 typedef int LandauIdx; 113 typedef struct { 114 PetscReal scale; 115 LandauIdx gid; // Lanadu matrix index (<10,000) 116 } pointInterpolationP4est; 117 typedef struct _lP4estVertexMaps { 118 LandauIdx (*gIdx)[LANDAU_MAX_SPECIES][LANDAU_MAX_NQ]; // #elems * LANDAU_MAX_NQ (spoof for max , Nb) on device, 119 LandauIdx num_elements; 120 LandauIdx num_reduced; 121 LandauIdx num_face; // (Q or Q^2 for 3D) 122 LandauDeviceType deviceType; 123 PetscInt Nf; 124 PetscInt Nq; 125 pointInterpolationP4est (*c_maps)[LANDAU_MAX_Q_FACE]; 126 struct _lP4estVertexMaps*data; 127 void *vp1,*vp2,*vp3; 128 } P4estVertexMaps; 129 130 typedef PetscReal LandauIPReal; 131 typedef struct { 132 LandauIPReal *w; 133 LandauIPReal *x; 134 LandauIPReal *y; 135 LandauIPReal *z; 136 LandauIPReal *coefs; 137 int dim_,ns_,nip_; 138 } LandauIPData; 139 140 PETSC_EXTERN PetscErrorCode LandauCreateColoring(Mat, DM, PetscContainer *); 141 PETSC_EXTERN PetscErrorCode LandauFormJacobian_Internal(Vec, Mat, const PetscInt, void *); 142 PETSC_EXTERN int LandauGetIPDataSize(const LandauIPData *const); 143 #if defined(PETSC_HAVE_CUDA) 144 PETSC_EXTERN PetscErrorCode LandauCUDAJacobian(DM, const PetscInt, const PetscReal [], const PetscReal [], const PetscReal[], const PetscReal[], 145 const LandauIPData *const, const PetscReal [], const PetscLogEvent[], Mat); 146 PETSC_EXTERN PetscErrorCode LandauCUDACreateMatMaps(P4estVertexMaps *, pointInterpolationP4est (*)[LANDAU_MAX_Q_FACE], PetscInt, PetscInt); 147 PETSC_EXTERN PetscErrorCode LandauCUDADestroyMatMaps(P4estVertexMaps *); 148 149 #endif 150 #if defined(PETSC_HAVE_KOKKOS) 151 /* TODO: this won't work if PETSc is built with C++ */ 152 #if !defined(__cplusplus) 153 PETSC_EXTERN PetscErrorCode LandauKokkosJacobian(DM, const PetscInt, const PetscReal [], const PetscReal [], const PetscReal[], const PetscReal[], 154 const LandauIPData *const, const PetscReal [],const PetscInt, const PetscLogEvent[], Mat); 155 PETSC_EXTERN PetscErrorCode LandauKokkosCreateMatMaps(P4estVertexMaps *, pointInterpolationP4est (*)[LANDAU_MAX_Q_FACE], PetscInt, PetscInt); 156 PETSC_EXTERN PetscErrorCode LandauKokkosDestroyMatMaps(P4estVertexMaps *); 157 #endif 158 #endif 159 160 #endif /* PETSCLANDAU_H */ 161