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();
    auto cycle_print_rate = mpl.postpro_cyclePrintRate();

#ifdef _OPENMP
    omp_set_num_threads(4);
#endif

    prepare_sources(mod, state, mpl);
    for(size_t i = 0; i < num_timesteps; i++)
    {
        mpl.call_timestep_callback(i);
        timestep(state);
        do_sources(mod, state, mpl);
        std::stringstream ss;
        ss << mpl.sim_name() << "/timestep_" << i << ".silo";

        if (i%silo_output_rate == 0)
            maxwell::export_to_silo(mod, state, ss.str());

        if (i%cycle_print_rate == 0)
        {
            std::cout << "Cycle " << i << ": t = " << state.curr_time << " s";
            std::cout << std::endl;
        }
        swap(state);
    }
}

int main(int argc, const char *argv[])
{
    gmsh::initialize();
    //make_geometry(1, 0.1);

    gmsh::option::setNumber("General.Terminal", 0);

    if (argc != 2)
    {
        std::cout << "Please specify configuration file" << std::endl;
        return 1;
    }

    maxwell::parameter_loader mpl;
    if ( not mpl.load_file( argv[1] ) )
    {
        std::cout << "Configuration problem, exiting" << std::endl;
        return 1;
    }

    std::filesystem::create_directory( mpl.sim_name() );;

    gmsh::open( mpl.sim_gmshmodel() );

    model mod( mpl.sim_geomorder(), mpl.sim_approxorder() );
    if ( not mpl.validate_model_params(mod) )
    {
        std::cout << "Configuration problem, exiting" << std::endl;
        return 1;
    }

    if ( mpl.sim_usegpu() )
    {
#ifdef ENABLE_GPU_SOLVER
        maxwell::solver_state_gpu state_g;
        test_it(mod, state_g, mpl);
#else
        std::cout << "GPU solver not compiled. Exiting." << std::endl;
#endif /* ENABLE_GPU_SOLVER */
    }
    else
    {
        maxwell::solver_state state_c;
        test_it(mod, state_c, mpl);
    }

    //gmsh::finalize();

    return 0;
}