XGCa
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
DM_wrapper.hpp
Go to the documentation of this file.
1 #ifndef DM_WRAPPER_HPP
2 #define DM_WRAPPER_HPP
3 
4 #define VELOCITY_SPACE_DIMENSION 2
5 
6 #ifdef NO_PETSC
7 // To avoid having to check for PETSc elsewhere define an empty struct if no PETSc
8 struct DMWrapper{
9 };
10 #else
11 
12 #include "sml.hpp"
13 #include <Kokkos_Core.hpp>
14 #include "space_settings.hpp"
15 #include <petscdm.h>
16 #include <petscversion.h>
17 #include <petscsys.h>
18 #include <petscdmplex.h>
19 #if PETSC_VERSION_GE(3,14,0)
20 #include <petscds.h>
21 #include "petscdmswarm.h"
22 #include "petscksp.h"
23 #endif
24 
25 // For some reason Kokkos cannot handle Views of PETSc objects, but it works when putting them in a wrapper
26 // Pseudo-inverse mesh PETSc objects
27 struct DMObjects{
28  DM dm;
29  Mat MM;
30 };
31 
32 // Wrapper for PETSc DM object to facilitate cleanup in destructor
33 struct DMWrapper{
34  bool is_owner;
35  int dm_len;
36  Kokkos::View<DMObjects*, Kokkos::LayoutRight, HostType> obj;
37 
39  : is_owner(false),
40  dm_len(0){
41  Kokkos::resize(obj, 0);
42  }
43 
44  ~DMWrapper(){ // destructor
45  if(is_owner){
46  for (int id = 0; id < dm_len; id++){
47  DMDestroy(&obj(id).dm);
48  MatDestroy(&obj(id).MM);
49  }
50  }
51  }
52 
53  DMWrapper(const DMWrapper& original) // copy constructor
54  : is_owner(false),
55  dm_len(original.dm_len){
56  Kokkos::resize(obj, dm_len);
57  for (int id = 0; id < dm_len; id++){
58  obj(id).dm = original.obj(id).dm;
59  obj(id).MM = original.obj(id).MM;
60  }
61  }
62 
63  DMWrapper& operator=(const DMWrapper& original) { // copy assignment
64  if (this == &original){
65  return *this;
66  }
67  is_owner = false;
68  dm_len = original.dm_len;
69  Kokkos::resize(obj, dm_len);
70  for (int id = 0; id < dm_len; id++){
71  obj(id).dm = original.obj(id).dm;
72  obj(id).MM = original.obj(id).MM;
73  }
74  return *this;
75  }
76 
94  PetscErrorCode initialize_pseudo_inverse_interpolation(int nthreads,
95  PetscInt order, const PetscInt Nspecies,
96  const PetscReal vPara_min, const PetscReal vPara_max,
97  const PetscInt num_vPara, const PetscReal vPerp_min,
98  const PetscReal vPerp_max, const PetscInt num_vPerp)
99  {
100 #if PETSC_VERSION_LT(3,14,0)
101  // To be able to compile with older PETSc versions
102  printf("\n\n\nERROR: initialize_pseudo_inverse_interpolation (should never reach here)\n");
103  return 0;
104 #else
105  is_owner = true;
106  dm_len = nthreads;
107  Kokkos::resize(obj, dm_len);
108 
109  PetscErrorCode ierr;
110 
111  PetscInt cells[2];
112  PetscReal lo[2], hi[2];
113  cells[0] = num_vPara - 1; //2 grid points give 1 cell, etc
114  cells[1] = num_vPerp - 1;
115  lo[0] = vPara_min;
116  lo[1] = vPerp_min;
117  hi[0] = vPara_max;
118  hi[1] = vPerp_max;
119 
120  PetscFunctionBeginUser;
121 
122  char order_str[50];
123  sprintf(order_str, "%d", order);
124  ierr = PetscOptionsSetValue(NULL, "-petscspace_degree", order_str);CHKERRQ(ierr);
125 
126  /* create mesh */
127  PetscFE fe;
128  PetscInt field = 0, dim = VELOCITY_SPACE_DIMENSION;
129  PetscBool interpolate = PETSC_TRUE;
130 
131  for (int id = 0; id < dm_len; id++){
132  ierr = DMPlexCreateBoxMesh(PETSC_COMM_SELF, dim, PETSC_FALSE, cells, lo, hi, NULL, interpolate, &obj(id).dm);CHKERRQ(ierr);
133  ierr = PetscObjectSetName((PetscObject)obj(id).dm, "Potential Grid");CHKERRQ(ierr);
134  if (id == dm_len - 1) {ierr = DMViewFromOptions(obj(id).dm, NULL, "-dm_view");CHKERRQ(ierr);}
135  ierr = DMSetFromOptions(obj(id).dm);CHKERRQ(ierr);
136  ierr = PetscFECreateDefault(PETSC_COMM_SELF, dim, Nspecies, PETSC_FALSE, "", PETSC_DECIDE, &fe);CHKERRQ(ierr);
137  //ierr = PetscFECreateLagrange(PETSC_COMM_SELF, dim, Nspecies, PETSC_FALSE, order, PETSC_DETERMINE, &fe);CHKERRQ(ierr);
138  ierr = PetscFESetFromOptions(fe);CHKERRQ(ierr);
139  ierr = PetscObjectSetName((PetscObject)fe, "fe");CHKERRQ(ierr);
140  ierr = DMSetField(obj(id).dm, field, NULL, (PetscObject)fe);CHKERRQ(ierr);
141  ierr = DMCreateDS(obj(id).dm);CHKERRQ(ierr);
142  ierr = PetscFEDestroy(&fe);CHKERRQ(ierr);
143 
144  ierr = DMCreateMassMatrix(obj(id).dm, obj(id).dm, &obj(id).MM);CHKERRQ(ierr);
145  }
146 
147  PetscFunctionReturn(0);
148 #endif
149 //PETSC_VERSION_LT(3,14,0)
150  } //end initialize_pseudo_inverse_interpolation()
151 
158  PetscErrorCode destroy_DM(){
159  PetscErrorCode ierr;
160  PetscFunctionBeginUser;
161  if(is_owner){
162  for (int id = 0; id < dm_len; id++){
163  if(obj(id).dm != NULL){
164  ierr = DMDestroy(&obj(id).dm);CHKERRQ(ierr);
165  obj(id).dm = NULL;
166  }
167  if(obj(id).MM != NULL){
168  ierr = MatDestroy(&obj(id).MM);CHKERRQ(ierr);
169  obj(id).MM = NULL;
170  }
171  }
172  }
173  PetscFunctionReturn(0);
174  }
175 };
176 #endif
177 //NO_PETSC
178 
179 #endif
PetscErrorCode initialize_pseudo_inverse_interpolation(int nthreads, PetscInt order, const PetscInt Nspecies, const PetscReal vPara_min, const PetscReal vPara_max, const PetscInt num_vPara, const PetscReal vPerp_min, const PetscReal vPerp_max, const PetscInt num_vPerp)
Definition: DM_wrapper.hpp:94
DMWrapper()
Definition: DM_wrapper.hpp:38
DM dm
PETSc DM mesh object representing the velocity grid when pseudo-inverse is used.
Definition: DM_wrapper.hpp:28
#define VELOCITY_SPACE_DIMENSION
Definition: DM_wrapper.hpp:4
bool is_owner
Definition: DM_wrapper.hpp:34
~DMWrapper()
Definition: DM_wrapper.hpp:44
Definition: DM_wrapper.hpp:33
DMWrapper(const DMWrapper &original)
Definition: DM_wrapper.hpp:53
DMWrapper & operator=(const DMWrapper &original)
Definition: DM_wrapper.hpp:63
Kokkos::View< DMObjects *, Kokkos::LayoutRight, HostType > obj
Stores #threads DMObjects.
Definition: DM_wrapper.hpp:36
PetscErrorCode destroy_DM()
Definition: DM_wrapper.hpp:158
int dm_len
Number of PETSc mesh objects to store.
Definition: DM_wrapper.hpp:35
Definition: DM_wrapper.hpp:27
Mat MM
PETSc mass matrix object for the velocity grid when pseudo-inverse is used.
Definition: DM_wrapper.hpp:29