DM_wrapper.hpp

Go to the documentation of this file.

#ifndef DM_WRAPPER_HPP
#define DM_WRAPPER_HPP

#define VELOCITY_SPACE_DIMENSION 2

#ifdef NO_PETSC
// To avoid having to check for PETSc elsewhere define an empty struct if no PETSc
struct DMWrapper{
};
#else

#include "sml.hpp"
#include <Kokkos_Core.hpp>
#include "space_settings.hpp"
#include <petscdm.h>
#include <petscversion.h>
#include <petscsys.h>
#include <petscdmplex.h>
#if PETSC_VERSION_GE(3,14,0)
#include <petscds.h>
#include "petscdmswarm.h"
#include "petscksp.h"
#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; 
};

// Wrapper for PETSc DM object to facilitate cleanup in destructor
struct DMWrapper{
    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(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)
    {
#if PETSC_VERSION_LT(3,14,0)
        // To be able to compile with older PETSc versions
        printf("\n\n\nERROR: initialize_pseudo_inverse_interpolation (should never reach here)\n");
        return 0;
#else
        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;

        char           order_str[50];
        sprintf(order_str, "%d", order);
        ierr = PetscOptionsSetValue(NULL, "-petscspace_degree", order_str);CHKERRQ(ierr);

        /* create mesh */
        PetscFE        fe;
        PetscInt       field = 0, dim = VELOCITY_SPACE_DIMENSION;
        PetscBool      interpolate = PETSC_TRUE;

        for (int id = 0; id < dm_len; id++){
            ierr = DMPlexCreateBoxMesh(PETSC_COMM_SELF, dim, PETSC_FALSE, cells, lo, hi, NULL, interpolate, &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 = PetscFECreateLagrange(PETSC_COMM_SELF, dim, Nspecies, PETSC_FALSE, order, PETSC_DETERMINE, &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(0);
#endif
//PETSC_VERSION_LT(3,14,0)
    } //end initialize_pseudo_inverse_interpolation()

    PetscErrorCode destroy_DM(){
        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(0);
    }
};
#endif
//NO_PETSC

#endif