Actual source code: petsclandau.h

  1: #if !defined(PETSCLANDAU_H)
  2: #define PETSCLANDAU_H

  4: #include <petscdmplex.h>
  5: #include <petscts.h>

  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[], PetscInt, 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 *);

 15: /* the Fokker-Planck-Landau context */
 16: #if !defined(LANDAU_DIM)
 17: #define LANDAU_DIM 2
 18: #endif

 20: #if !defined(LANDAU_MAX_SPECIES)
 21: #if LANDAU_DIM==2
 22: #define LANDAU_MAX_SPECIES 10
 23: #define LANDAU_MAX_GRIDS 3
 24: #else
 25: #define LANDAU_MAX_SPECIES 3
 26: #define LANDAU_MAX_GRIDS 2
 27: #endif
 28: #else
 29: #define LANDAU_MAX_GRIDS 3
 30: #endif

 32: #if !defined(LANDAU_MAX_Q)
 33: #if defined(LANDAU_MAX_NQ)
 34: #error"LANDAU_MAX_NQ but not LANDAU_MAX_Q. Use -DLANDAU_MAX_Q=4 for Q3 elements"
 35: #endif
 36: #if LANDAU_DIM==2
 37: #define LANDAU_MAX_Q 5
 38: #else
 39: #define LANDAU_MAX_Q 3
 40: #endif
 41: #else
 42: #undef LANDAU_MAX_NQ
 43: #endif

 45: #if LANDAU_DIM==2
 46: #define LANDAU_MAX_Q_FACE LANDAU_MAX_Q
 47: #define LANDAU_MAX_NQ (LANDAU_MAX_Q*LANDAU_MAX_Q)
 48: #else
 49: #define LANDAU_MAX_Q_FACE (LANDAU_MAX_Q*LANDAU_MAX_Q)
 50: #define LANDAU_MAX_NQ (LANDAU_MAX_Q*LANDAU_MAX_Q*LANDAU_MAX_Q)
 51: #endif

 53: typedef enum {LANDAU_CUDA, LANDAU_KOKKOS, LANDAU_CPU} LandauDeviceType;

 55: // static data
 56: typedef struct {
 57:   void  *invJ;  // nip*dim*dim
 58:   void  *D;     // nq*nb*dim
 59:   void  *B;     // nq*nb
 60:   void  *alpha; // ns
 61:   void  *beta;  // ns
 62:   void  *invMass; // ns
 63:   void  *w; // nip
 64:   void  *x; // nip
 65:   void  *y; // nip
 66:   void  *z; // nip
 67:   void  *Eq_m; // ns - dynamic
 68:   void  *f; //  nip*Nf - dynamic (IP)
 69:   void  *dfdx; // nip*Nf - dynamic (IP)
 70:   void  *dfdy; // nip*Nf - dynamic (IP)
 71:   void  *dfdz; // nip*Nf - dynamic (IP)
 72:   int   dim_,ns_,nip_,nq_,nb_;
 73:   void  *NCells;
 74:   void  *species_offset;
 75:   void  *mat_offset;
 76:   void  *elem_offset;
 77:   void  *ip_offset;
 78:   void  *ipf_offset;
 79:   void  *ipfdf_data;
 80:   void  *maps;
 81: } LandauStaticData;

 83: typedef struct {
 84:   PetscBool      interpolate;                  /* Generate intermediate mesh elements */
 85:   PetscBool      gpu_assembly;
 86:   MPI_Comm       comm; /* global communicator to use for errors and diagnostics */
 87:   double         times[1];
 88:   PetscBool      initialized;
 89:   PetscBool      use_matrix_mass;
 90:   /* FE */
 91:   PetscFE        fe[LANDAU_MAX_SPECIES];
 92:   /* geometry  */
 93:   PetscReal      radius[LANDAU_MAX_GRIDS];
 94:   PetscReal      re_radius;           /* RE: radius of refinement along v_perp=0, z>0 */
 95:   PetscReal      vperp0_radius1;      /* RE: radius of refinement along v_perp=0 */
 96:   PetscReal      vperp0_radius2;      /* RE: radius of refinement along v_perp=0 after origin AMR refinement */
 97:   PetscBool      sphere; // not used
 98:   PetscBool      inflate; // not used
 99:   PetscReal      i_radius[LANDAU_MAX_GRIDS]; // not used
100:   PetscReal      e_radius; // not used
101:   PetscInt       num_sections; // not used
102:   /* AMR */
103:   PetscBool      use_p4est;
104:   PetscInt       numRERefine;        /* RE: refinement along v_perp=0, z > 0 */
105:   PetscInt       nZRefine1;          /* RE: origin refinement after v_perp=0 refinement */
106:   PetscInt       nZRefine2;          /* RE: origin refinement after origin AMR refinement */
107:   PetscInt       numAMRRefine[LANDAU_MAX_GRIDS];  /* normal AMR - refine from origin */
108:   PetscInt       postAMRRefine[LANDAU_MAX_GRIDS]; /* uniform refinement of AMR */
109:   /* relativistic */
110:   PetscBool      use_energy_tensor_trick;
111:   PetscBool      use_relativistic_corrections;
112:   /* physics */
113:   PetscReal      thermal_temps[LANDAU_MAX_SPECIES];
114:   PetscReal      masses[LANDAU_MAX_SPECIES];  /* mass of each species  */
115:   PetscReal      charges[LANDAU_MAX_SPECIES]; /* charge of each species  */
116:   PetscReal      n[LANDAU_MAX_SPECIES];       /* number density of each species  */
117:   PetscReal      m_0;      /* reference mass */
118:   PetscReal      v_0;      /* reference velocity */
119:   PetscReal      n_0;      /* reference number density */
120:   PetscReal      t_0;      /* reference time */
121:   PetscReal      Ez;
122:   PetscReal      epsilon0;
123:   PetscReal      k;
124:   PetscReal      lnLam;
125:   PetscReal      electronShift;
126:   PetscInt       num_species;
127:   PetscInt       num_grids;
128:   PetscInt       species_offset[LANDAU_MAX_GRIDS+1];
129:   PetscInt       mat_offset[LANDAU_MAX_GRIDS+1];
130:   /* cache */
131:   Mat            J;
132:   Mat            M;
133:   Vec            X;
134:   /* derived type */
135:   void          *data;
136:   PetscBool      aux_bool;  /* helper */
137:   /* computing */
138:   LandauDeviceType deviceType;
139:   PetscInt       subThreadBlockSize;
140:   PetscInt       numConcurrency; /* number of SMs in Cuda to use */
141:   DM             pack;
142:   DM             plex[LANDAU_MAX_GRIDS];
143:   LandauStaticData SData_d; /* static geometric data on device, this holds host pointers */
144:   /* diagnostics */
145:   PetscInt       verbose;
146:   PetscLogEvent  events[20];
147: } LandauCtx;

149: typedef int LandauIdx;
150: typedef struct {
151:   PetscReal scale;
152:   LandauIdx gid;   // Lanadu matrix index (<10,000)
153: } pointInterpolationP4est;
154: typedef struct _lP4estVertexMaps {
155:   LandauIdx                (*gIdx)[LANDAU_MAX_SPECIES][LANDAU_MAX_NQ]; // #elems *  LANDAU_MAX_NQ (spoof for max , Nb) on device,
156:   LandauIdx                num_elements;
157:   LandauIdx                num_reduced;
158:   LandauIdx                num_face;  // (Q or Q^2 for 3D)
159:   LandauDeviceType         deviceType;
160:   PetscInt                 Nf;
161:   PetscInt                 Nq;
162:   pointInterpolationP4est (*c_maps)[LANDAU_MAX_Q_FACE];
163:   struct _lP4estVertexMaps*d_self;
164:   void                    *vp1,*vp2,*vp3;
165:   PetscInt                numgrids;
166: } P4estVertexMaps;

168: PETSC_EXTERN PetscErrorCode LandauCreateColoring(Mat, DM, PetscContainer *);
169: #if defined(PETSC_HAVE_CUDA)
170: PETSC_EXTERN PetscErrorCode LandauCUDAJacobian(DM[], const PetscInt, const PetscInt, const PetscInt[], PetscReal[], PetscScalar[], const PetscInt, const PetscScalar[],
171:                                                const LandauStaticData *, const PetscInt, const PetscReal, const PetscLogEvent[],const PetscInt[], const PetscInt[], Mat[], Mat);
172: PETSC_EXTERN PetscErrorCode LandauCUDACreateMatMaps(P4estVertexMaps *, pointInterpolationP4est (*)[LANDAU_MAX_Q_FACE], PetscInt[], PetscInt, PetscInt);
173: PETSC_EXTERN PetscErrorCode LandauCUDADestroyMatMaps(P4estVertexMaps *, PetscInt);
174: PETSC_EXTERN PetscErrorCode LandauCUDAStaticDataSet(DM, const PetscInt, const PetscInt, PetscInt[], PetscInt[], PetscInt[], PetscReal [], PetscReal [], PetscReal[], PetscReal[], PetscReal[], PetscReal[], PetscReal[],
175:                                                     PetscReal[], LandauStaticData *);
176: PETSC_EXTERN PetscErrorCode LandauCUDAStaticDataClear(LandauStaticData *);
177: #endif
178: #if defined(PETSC_HAVE_KOKKOS)
179: PETSC_EXTERN PetscErrorCode LandauKokkosJacobian(DM[], const PetscInt, const PetscInt, const PetscInt[], PetscReal[], PetscScalar[],  const PetscInt, const PetscScalar[],
180:                                                  const LandauStaticData *, const PetscInt, const PetscReal, const PetscLogEvent[], const PetscInt[], const PetscInt[], Mat[], Mat);
181: PETSC_EXTERN PetscErrorCode LandauKokkosCreateMatMaps(P4estVertexMaps *, pointInterpolationP4est (*)[LANDAU_MAX_Q_FACE], PetscInt[], PetscInt, PetscInt);
182: PETSC_EXTERN PetscErrorCode LandauKokkosDestroyMatMaps(P4estVertexMaps *, PetscInt);
183: PETSC_EXTERN PetscErrorCode LandauKokkosStaticDataSet(DM, const PetscInt, const PetscInt, PetscInt[], PetscInt[], PetscInt[], PetscReal [], PetscReal [], PetscReal[], PetscReal[], PetscReal[], PetscReal[], PetscReal[],
184:                                                       PetscReal[], LandauStaticData *);
185: PETSC_EXTERN PetscErrorCode LandauKokkosStaticDataClear(LandauStaticData*);
186: #endif

188: #endif /* PETSCLANDAU_H */