maxwell_solver.cpp

#include <filesystem>

#ifdef _OPENMP
#include <omp.h>
#endif

#include "gmsh.h"
#include "silo_output.hpp"
#include "maxwell_interface.h"
#include "param_loader.h"

#if 0
static void
make_geometry(int order, double mesh_h)
{
    gm::add("difftest");
    /*
    std::vector<std::pair<int,int>> objects;

    objects.push_back(
        std::make_pair(3, gmo::addBox(0.0, 0.0, 0.0, 1.0, 1.0, 0.05) )
        );
    objects.push_back(
        std::make_pair(3, gmo::addBox(0.0, 0.0, 0.05, 1.0, 1.0, 0.05) )
        );

    std::vector<std::pair<int, int>> tools;
    gmsh::vectorpair odt;
    std::vector<gmsh::vectorpair> odtm;
    gmo::fragment(objects, tools, odt, odtm);
    */
    gmo::addBox(0.0, 0.0, 0.0, 1.0, 1.0, 0.1);
    
    gmo::synchronize();

    gvp_t vp;
    gm::getEntities(vp);
    gmm::setSize(vp, mesh_h);
}
#endif

template<typename State>
void initialize_solver(const model& mod, State& state, const maxwell::parameter_loader& mpl)
{
    maxwell::init_from_model(mod, state);
    
    if ( mpl.initial_Efield_defined() )
    {
        auto E = [&](const point_3d& pt) -> vec3d { return mpl.initial_Efield(pt); };
        maxwell::init_E_field(mod, state, E);
    }

    if ( mpl.initial_Hfield_defined() )
    {
        auto H = [&](const point_3d& pt) -> vec3d { return mpl.initial_Hfield(pt); };
        maxwell::init_H_field(mod, state, H);
    }

    init_matparams(mod, state, mpl);
}


template<typename State>
void test_it(const model& mod, State& state, maxwell::parameter_loader& mpl)
{
    initialize_solver(mod, state, mpl);

    state.delta_t = mpl.sim_dt();
    auto num_timesteps = mpl.sim_timesteps();
    auto silo_output_rate = mpl.postpro_siloOutputRate();