maxwell_interface.h

#pragma once

/* This file is included by the solver code and by the kernels.
 * The solver code needs all the parameter_loader stuff, which
 * ends up including sol2. Sol2 is C++17, which nvcc does not
 * support. Since the kernels need only the data structure
 * definitions, we hide all the other stuff when compiling
 * with nvcc. */

#include "gmsh_io.h"
#include "types.h"

#ifdef ENABLE_GPU_SOLVER
#include "gpu_support.hpp"
#endif /* ENABLE_GPU_SOLVER */

#ifndef HIDE_THIS_FROM_NVCC
#include "param_loader.h"
#endif

namespace maxwell {

#ifndef HIDE_THIS_FROM_NVCC
enum class boundary_condition
{
    UNSPECIFIED,
    PEC,
    PMC,
    IMPEDANCE,
    PLANE_WAVE_E,
    PLANE_WAVE_H,
    E_FIELD,
    H_FIELD,
    SURFACE_CURRENT
};

enum class interface_condition
{
    UNSPECIFIED,
    E_FIELD,
    H_FIELD,
    SURFACE_CURRENT
};

class parameter_loader : public parameter_loader_base
{
    void init(void);

    bool    validate_materials(const std::string&, int) const;
    bool    validate_boundary_condition(const model&, int) const;
    bool    validate_interface_condition(const model&, int) const;
    bool    validate_conditions(const model&) const;

public:
    parameter_loader();

    bool    validate_model_params(const model&) const;
    double  epsilon(int, const point_3d&) const;
    double  mu(int, const point_3d&) const;
    double  sigma(int, const point_3d&) const;

    bool    initial_Efield_defined(void) const;
    vec3d   initial_Efield(const point_3d&) const;

    bool    initial_Hfield_defined(void) const;
    vec3d   initial_Hfield(const point_3d&) const;

    vec3d   eval_boundary_source(int, const point_3d&, double) const;
    vec3d   eval_interface_source(int, const point_3d&, double) const;

    boundary_condition  boundary_type(int tag) const;
    interface_condition interface_type(int tag) const;
    field_export_mode postpro_fieldExportMode(const std::string&) const;
};
#endif /* HIDE_THIS_FROM_NVCC */

/************************************************
 * FIELDS
 */
struct field
{
    vecxd   Ex;
    vecxd   Ey;
    vecxd   Ez;
    vecxd   Hx;
    vecxd   Hy;
    vecxd   Hz;

    size_t  num_dofs;

    void    zero(void);
    void    resize(size_t);
};

#ifdef ENABLE_GPU_SOLVER
struct pinned_field
{
    double  *Ex;
    double  *Ey;
    double  *Ez;
    double  *Hx;
    double  *Hy;
    double  *Hz;

    size_t  num_dofs;

    pinned_field();
    ~pinned_field();
    pinned_field(const pinned_field&) = delete;
    pinned_field& operator=(const pinned_field&) = delete;
    
    void    zero(void);
    void    resize(size_t);
};
#endif /* ENABLE_GPU_SOLVER */

#ifdef ENABLE_GPU_SOLVER
struct field_gpu
{
    device_vector<double>   Ex;
    device_vector<double>   Ey;
    device_vector<double>   Ez;
    device_vector<double>   Hx;
    device_vector<double>   Hy;
    device_vector<double>   Hz;

    size_t num_dofs;

    template<typename T>
    struct priv_raw_ptrs
    {
        size_t  num_dofs;

        T  *Ex;
        T  *Ey;
        T  *Ez;
        T  *Hx;
        T  *Hy;
        T  *Hz;
    };

    using raw_ptrs = priv_raw_ptrs<double>;
    using const_raw_ptrs = priv_raw_ptrs<const double>;

    void            zero(void);
    void            resize(size_t);
    raw_ptrs        data(void);
    const_raw_ptrs  data(void) const;
    void            copyin(const field&);
    void            copyin(const field&, const stream&);
    void            copyin(const pinned_field&, const stream&);
    void            copyout(field&) const;

};
#endif /* ENABLE_GPU_SOLVER */

/************************************************
 * MATERIAL PARAMETERS
 */
struct material_params
{
    size_t  num_dofs;
    size_t  num_fluxes;

    vecxd   inv_epsilon;
    vecxd   inv_mu;
    vecxd   sigma;
    vecxd   sigma_over_epsilon;

    /* Lax-Milgram flux coefficients */
    vecxd   aE;
    vecxd   bE;
    vecxd   aH;
    vecxd   bH;

    /* Boundary conditions coefficients */
    vecxd   bc_coeffs;
};

#ifdef ENABLE_GPU_SOLVER
struct material_params_gpu
{
    size_t  num_dofs;
    size_t  num_fluxes;

    device_vector<double>   inv_epsilon;
    device_vector<double>   inv_mu;
    device_vector<double>   sigma;
    device_vector<double>   sigma_over_epsilon;

    /* Lax-Milgram flux coefficients */
    device_vector<double>   aE;
    device_vector<double>   bE;
    device_vector<double>   aH;
    device_vector<double>   bH;

    /* Boundary conditions coefficients */
    device_vector<double>   bc_coeffs;

    template<typename T>
    struct priv_raw_ptrs
    {
        size_t num_dofs;
        size_t num_fluxes;

        T *inv_epsilon;
        T *inv_mu;
        T *sigma;
        T *sigma_over_epsilon;
        T *aE;
        T *bE;
        T *aH;
        T *bH;
        T *bc_coeffs; 
    };

    using raw_ptrs = priv_raw_ptrs<double>;
    using const_raw_ptrs = priv_raw_ptrs<const double>;

    raw_ptrs        data(void);
    const_raw_ptrs  data(void) const;
    void            copyin(const material_params&);
};
#endif

/************************************************
 * SOLVER STATE
 */
struct solver_state
{
    std::vector<entity_data_cpu>        eds;
    field                               emf_curr;
    field                               emf_next;
    material_params                     matparams;

    field                               jumps;
    field                               fluxes;

    double                              delta_t;
    double                              curr_time;
    size_t                              curr_timestep;

    field                               dx;
    field                               dy;
    field                               dz;

    field                               k1;
    field                               k2;
    field                               k3;
    field                               k4;
    field                               tmp;

    field                               bndsrcs;
};

#ifdef ENABLE_GPU_SOLVER
struct solver_state_gpu
{
    std::vector<entity_data_cpu>        eds;
    std::vector<entity_data_gpu>        edgs;
    field_gpu                           emf_curr;
    field_gpu                           emf_next;
    material_params_gpu                 matparams;

    field_gpu                           jumps;
    field_gpu                           fluxes;

    double                              delta_t;
    double                              curr_time;
    size_t                              curr_timestep;

    field_gpu                           dx;
    field_gpu                           dy;
    field_gpu                           dz;

    field_gpu                           k1;
    field_gpu                           k2;
    field_gpu                           k3;
    field_gpu                           k4;
    field_gpu                           tmp;

    field_gpu                           bndsrcs;
    field_gpu                           bndsrcs_buf;

    field                               bndsrcs_decomp_cpu;
    pinned_field                        bndsrcs_cpu;

    std::vector<size_t>                 bndsrcs_decomp_table_cpu;
    device_vector<size_t>               bndsrcs_decomp_table;

    stream                              memcpy_stream;
    stream                              compute_stream;
};
#endif /* ENABLE_GPU_SOLVER */

#ifndef HIDE_THIS_FROM_NVCC
void init_from_model(const model&, solver_state&);
void init_matparams(const model&, solver_state&, const parameter_loader&);
//void apply_operator(solver_state&, const field&, field&);
void export_fields_to_silo(const model&, const solver_state&, const parameter_loader&);
void timestep(solver_state&, time_integrator_type);
void prepare_sources(const model&, maxwell::solver_state&, maxwell::parameter_loader&);
void do_sources(const model&, maxwell::solver_state&, maxwell::parameter_loader&);
void swap(maxwell::solver_state&);

#ifdef ENABLE_GPU_SOLVER
void init_from_model(const model&, solver_state_gpu&);
void init_matparams(const model&, solver_state_gpu&, const parameter_loader&);
//void apply_operator(solver_state_gpu&, const field_gpu&, field_gpu&);
void export_fields_to_silo(const model&, const solver_state_gpu&, const parameter_loader&);
void timestep(solver_state_gpu&, time_integrator_type);
void prepare_sources(const model&, maxwell::solver_state_gpu&, maxwell::parameter_loader&);
void do_sources(const model&, maxwell::solver_state_gpu&, maxwell::parameter_loader&);
void swap(maxwell::solver_state_gpu&);
#endif /* ENABLE_GPU_SOLVER */

void init_boundary_conditions(const model&, const parameter_loader&, vecxd&);
void init_lax_milgram(const model&, const parameter_loader&, vecxd&, vecxd&, vecxd&, vecxd&);
//void eval_boundary_sources(const model&, const parameter_loader&, solver_state&, field&);
//void eval_interface_sources(const model&, const parameter_loader&, solver_state&, field&);
#endif /* HIDE_THIS_FROM_NVCC */

template<typename Function>
void
init_E_field(const model& mod, solver_state& state, const Function& E)
{
    for (auto& e : mod)
        e.project(E, state.emf_curr.Ex, state.emf_curr.Ey, state.emf_curr.Ez);
}

template<typename Function>
void
init_H_field(const model& mod, solver_state& state, const Function& H)
{
    for (auto& e : mod)
        e.project(H, state.emf_curr.Hx, state.emf_curr.Hy, state.emf_curr.Hz);
}

#ifdef ENABLE_GPU_SOLVER
template<typename Function>
void
init_E_field(const model& mod, solver_state_gpu& state, const Function& E)
{
    vecxd Ex = vecxd::Zero( mod.num_dofs() );
    vecxd Ey = vecxd::Zero( mod.num_dofs() );
    vecxd Ez = vecxd::Zero( mod.num_dofs() );

    for (auto& e : mod)
        e.project(E, Ex, Ey, Ez);
    
    state.emf_curr.Ex.copyin(Ex.data(), Ex.size());
    state.emf_curr.Ey.copyin(Ey.data(), Ey.size());
    state.emf_curr.Ez.copyin(Ez.data(), Ez.size());
}

template<typename Function>
void
init_H_field(const model& mod, solver_state_gpu& state, const Function& H)
{
    vecxd Hx = vecxd::Zero( mod.num_dofs() );
    vecxd Hy = vecxd::Zero( mod.num_dofs() );
    vecxd Hz = vecxd::Zero( mod.num_dofs() );

    for (auto& e : mod)
        e.project(H, Hx, Hy, Hz);
    
    state.emf_curr.Hx.copyin(Hx.data(), Hx.size());
    state.emf_curr.Hy.copyin(Hy.data(), Hy.size());
    state.emf_curr.Hz.copyin(Hz.data(), Hz.size());
}


#endif /* ENABLE_GPU_SOLVER */

void gpu_compute_jumps(const entity_data_gpu&, const field_gpu&, field_gpu&,
    double *, gpuStream_t stream = 0);

void gpu_compute_jumps_E(const entity_data_gpu&, const field_gpu&, field_gpu&,
    double *, gpuStream_t stream = 0);

void gpu_compute_jumps_H(const entity_data_gpu&, const field_gpu&, field_gpu&,
    double *, gpuStream_t stream = 0);

void gpu_compute_fluxes(const entity_data_gpu&, const field_gpu&, field_gpu&,
    const field_gpu&, const material_params_gpu&, gpuStream_t stream = 0);

void gpu_compute_fluxes_E(const entity_data_gpu&, const field_gpu&, field_gpu&,
    const field_gpu&, const material_params_gpu&, gpuStream_t stream = 0);

void gpu_compute_fluxes_H(const entity_data_gpu&, const field_gpu&, field_gpu&,
    const field_gpu&, const material_params_gpu&, gpuStream_t stream = 0);

void decompress_bndsrc(const solver_state_gpu& state, const field_gpu& csrcs,
    field_gpu& srcs);

} // namespace maxwell