DM_wrapper.hpp

Go to the documentation of this file.

#ifndef DM_WRAPPER_HPP
#define DM_WRAPPER_HPP
 
#define VELOCITY_SPACE_DIMENSION 2
 
#ifndef NO_PETSC
 
#include <string>
#include <Kokkos_Core.hpp>
#include "space_settings.hpp"
#include <petscdm.h>
#include <petscversion.h>
#include <petscsys.h>
#include <petscdmplex.h>
#include <petscds.h>
#include "petscdmswarm.h"
#include "petscksp.h"
 
// For older PETSc versions
#if !defined(PETSC_USE_STRICT_PETSCERRORCODE) && !defined(PETSC_SUCCESS)
#define PETSC_SUCCESS 0
#endif
 
// For some reason Kokkos cannot handle Views of PETSc objects, but it works when putting them in a wrapper
// Pseudo-inverse mesh PETSc objects
struct DMObjects{
    DM  dm; 
    Mat MM; 
};
 
#endif
 
// Wrapper for PETSc DM object to facilitate cleanup in destructor
struct DMWrapper{
 
#ifdef NO_PETSC
    DMWrapper(){}
#else
    bool is_owner;
    int dm_len; 
    Kokkos::View<DMObjects*, Kokkos::LayoutRight, HostType> obj; 
 
    DMWrapper()
        : is_owner(false),
          dm_len(0){
        Kokkos::resize(obj, 0);
    }
 
    ~DMWrapper(){ // destructor
        if(is_owner){
            for (int id = 0; id < dm_len; id++){
                DMDestroy(&obj(id).dm);
                MatDestroy(&obj(id).MM);
            }
        }
    }
 
    DMWrapper(const DMWrapper& original) // copy constructor
        : is_owner(false),
          dm_len(original.dm_len){
        Kokkos::resize(obj, dm_len);
        for (int id = 0; id < dm_len; id++){
            obj(id).dm = original.obj(id).dm;
            obj(id).MM = original.obj(id).MM;
        }
    }
 
    DMWrapper& operator=(const DMWrapper& original) { // copy assignment
        if (this == &original){
            return *this;
        }
        is_owner = false;
        dm_len = original.dm_len;
        Kokkos::resize(obj, dm_len);
        for (int id = 0; id < dm_len; id++){
            obj(id).dm = original.obj(id).dm;
            obj(id).MM = original.obj(id).MM;
        }
        return *this;
    }
 
    PetscErrorCode initialize_pseudo_inverse_interpolation_petsc(int nthreads,
                             PetscInt order, const PetscInt Nspecies,
                             const PetscReal vPara_min, const PetscReal vPara_max,
                             const PetscInt nvp, const PetscReal vPerp_min,
                             const PetscReal vPerp_max, const PetscInt nmu)
    {
        int element_order_shift = (order > 0) ? order : 1;
        const PetscInt num_vPara = 2*nvp/element_order_shift + 1;
        const PetscInt num_vPerp = nmu/element_order_shift + 1;
 
        is_owner = true;
        dm_len = nthreads;
        Kokkos::resize(obj, dm_len);
 
        PetscErrorCode ierr;
 
        PetscInt       cells[2];
        PetscReal      lo[2], hi[2];
        cells[0] = num_vPara - 1; //2 grid points give 1 cell, etc
        cells[1] = num_vPerp - 1;
        lo[0] = vPara_min;
        lo[1] = vPerp_min;
        hi[0] = vPara_max;
        hi[1] = vPerp_max;
 
        PetscFunctionBeginUser;
 
        std::string    order_str = std::to_string(order);
        ierr = PetscOptionsSetValue(NULL, "-petscspace_degree", order_str.c_str());CHKERRQ(ierr);
 
        /* create mesh */
        PetscFE        fe;
        PetscInt       field = 0, dim = VELOCITY_SPACE_DIMENSION;
        DMBoundaryType periodicity[2] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE};
        PetscBool      interpolate = PETSC_TRUE;
 
        for (int id = 0; id < dm_len; id++){
            PetscInt petsc_zero=0.0;
            ierr = DMPlexCreateBoxMesh(PETSC_COMM_SELF, dim, PETSC_FALSE, cells, lo, hi, periodicity, interpolate, petsc_zero, PETSC_FALSE, &obj(id).dm);CHKERRQ(ierr);
            ierr = PetscObjectSetName((PetscObject)obj(id).dm, "Potential Grid");CHKERRQ(ierr);
            if (id == dm_len - 1) {ierr = DMViewFromOptions(obj(id).dm, NULL, "-dm_view");CHKERRQ(ierr);}
            ierr = DMSetFromOptions(obj(id).dm);CHKERRQ(ierr);
            ierr = PetscFECreateDefault(PETSC_COMM_SELF, dim, Nspecies, PETSC_FALSE, "", PETSC_DECIDE, &fe);CHKERRQ(ierr);
            ierr = PetscFESetFromOptions(fe);CHKERRQ(ierr);
            ierr = PetscObjectSetName((PetscObject)fe, "fe");CHKERRQ(ierr);
            ierr = DMSetField(obj(id).dm, field, NULL, (PetscObject)fe);CHKERRQ(ierr);
            ierr = DMCreateDS(obj(id).dm);CHKERRQ(ierr);
            ierr = PetscFEDestroy(&fe);CHKERRQ(ierr);
 
            ierr = DMCreateMassMatrix(obj(id).dm, obj(id).dm, &obj(id).MM);CHKERRQ(ierr);
        }
 
        PetscFunctionReturn(PETSC_SUCCESS);
    } //end initialize_pseudo_inverse_interpolation()
 
    PetscErrorCode destroy_DM_petsc(){
        PetscErrorCode ierr;
        PetscFunctionBeginUser;
        if(is_owner){
            for (int id = 0; id < dm_len; id++){
                if(obj(id).dm != NULL){
                    ierr = DMDestroy(&obj(id).dm);CHKERRQ(ierr);
                    obj(id).dm = NULL;
                }
                if(obj(id).MM != NULL){
                    ierr = MatDestroy(&obj(id).MM);CHKERRQ(ierr);
                    obj(id).MM = NULL;
                }
            }
        }
        PetscFunctionReturn(PETSC_SUCCESS);
    }
#endif
 
    void initialize_pseudo_inverse_interpolation(int nthreads,
                             int order, const int Nspecies,
                             const double vPara_min, const double vPara_max,
                             const int nvp, const double vPerp_min,
                             const double vPerp_max, const int nmu){
#ifndef NO_PETSC
        int ierr = initialize_pseudo_inverse_interpolation_petsc(nthreads, order, Nspecies, vPara_min, vPara_max, nvp, vPerp_min, vPerp_max, nmu);
        CHKERRABORT(PETSC_COMM_SELF,ierr);
#endif
    }
 
    void destroy_DM(){
#ifndef NO_PETSC
        int ierr = destroy_DM_petsc();
        CHKERRABORT(PETSC_COMM_SELF,ierr);
#endif
    }
};
 
#endif