diff --git a/CMakeLists.txt b/CMakeLists.txt
index 5e0068e5d4189093c358bf34c9987706d5815f1a..8d3def61aa2a8264b2b7670d93e9d348c45ea5d3 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -253,6 +253,8 @@ set(LIBGMSHDG_SOURCES src/gmsh_io.cpp
src/kernels_gpu_mi.cpp
src/kernels_cuda.cu
src/maxwell_interface.cpp
+ src/maxwell_cpu.cpp
+ src/maxwell_gpu.cpp
)
add_library(gmshdg SHARED ${LIBGMSHDG_SOURCES})
diff --git a/include/maxwell_interface.h b/include/maxwell_interface.h
index 2481d886eec061ed99bb2d9e921520fd30497b4e..371b57d401e1936e3523754538963d96bdefbbf2 100644
--- a/include/maxwell_interface.h
+++ b/include/maxwell_interface.h
@@ -1,9 +1,12 @@
#pragma once
+#include "gmsh_io.h"
#include "types.h"
#include "gpu_support.hpp"
-struct electromagnetic_field
+namespace maxwell {
+
+struct field
{
vecxd Ex;
vecxd Ey;
@@ -17,7 +20,7 @@ struct electromagnetic_field
void resize(size_t);
};
-struct electromagnetic_material_params
+struct material_params
{
size_t num_dofs;
size_t num_fluxes;
@@ -37,7 +40,7 @@ struct electromagnetic_material_params
vecxd bc_coeffs;
};
-struct electromagnetic_material_params_gpu
+struct material_params_gpu
{
size_t num_dofs;
size_t num_fluxes;
@@ -57,7 +60,7 @@ struct electromagnetic_material_params_gpu
device_vector<double> bc_coeffs;
};
-struct electromagnetic_field_gpu
+struct field_gpu
{
device_vector<double> Ex;
device_vector<double> Ey;
@@ -83,7 +86,50 @@ struct electromagnetic_field_gpu
void zero(void);
void resize(size_t);
raw_ptrs data(void);
- void copyin(const electromagnetic_field&);
- void copyout(electromagnetic_field&);
+ void copyin(const field&);
+ void copyout(field&);
+
+};
+
+struct solver_state
+{
+ std::vector<entity_data_cpu> eds;
+ field emf_curr;
+ field emf_next;
+ material_params matparams;
+
+ double delta_t;
+ double curr_time;
+ size_t curr_timestep;
+
+ field dx;
+ field dy;
+ field dz;
+};
+
+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 matparams;
+
+ double delta_t;
+ double curr_time;
+ size_t curr_timestep;
+
+ field_gpu dx;
+ field_gpu dy;
+ field_gpu dz;
+};
+
+void init_from_model(const model&, solver_state&);
+void apply_operator(solver_state&, const field&, field&);
+void export_to_silo(const model&, const solver_state&, const std::string&);
+
+void init_from_model(const model&, solver_state_gpu&);
+void apply_operator(solver_state_gpu&, const field_gpu&, field_gpu&);
+void export_to_silo(const model&, const field&, const std::string&);
-};
\ No newline at end of file
+} // namespace maxwell
\ No newline at end of file
diff --git a/src/maxwell_interface.cpp b/src/maxwell_interface.cpp
index 76b732ca6389bff7943f8edbbf8a4f1e95c115f0..a26f26ca83777a9b59a2ecb3b16b20fdce420a9a 100644
--- a/src/maxwell_interface.cpp
+++ b/src/maxwell_interface.cpp
@@ -1,7 +1,9 @@
#include "maxwell_interface.h"
+namespace maxwell {
+
void
-electromagnetic_field::resize(size_t p_num_dofs)
+field::resize(size_t p_num_dofs)
{
Ex = vecxd::Zero(p_num_dofs);
Ey = vecxd::Zero(p_num_dofs);
@@ -13,7 +15,7 @@ electromagnetic_field::resize(size_t p_num_dofs)
}
void
-electromagnetic_field_gpu::zero()
+field_gpu::zero()
{
Ex.zero();
Ey.zero();
@@ -24,7 +26,7 @@ electromagnetic_field_gpu::zero()
}
void
-electromagnetic_field_gpu::resize(size_t p_num_dofs)
+field_gpu::resize(size_t p_num_dofs)
{
Ex.resize(p_num_dofs);
Ey.resize(p_num_dofs);
@@ -35,8 +37,8 @@ electromagnetic_field_gpu::resize(size_t p_num_dofs)
num_dofs = p_num_dofs;
}
-electromagnetic_field_gpu::raw_ptrs
-electromagnetic_field_gpu::data(void)
+field_gpu::raw_ptrs
+field_gpu::data(void)
{
raw_ptrs ret;
@@ -52,7 +54,7 @@ electromagnetic_field_gpu::data(void)
}
void
-electromagnetic_field_gpu::copyin(const electromagnetic_field& emf)
+field_gpu::copyin(const field& emf)
{
Ex.copyin(emf.Ex.data(), emf.Ex.size());
Ey.copyin(emf.Ey.data(), emf.Ey.size());
@@ -63,7 +65,7 @@ electromagnetic_field_gpu::copyin(const electromagnetic_field& emf)
}
void
-electromagnetic_field_gpu::copyout(electromagnetic_field& emf)
+field_gpu::copyout(field& emf)
{
if (num_dofs != emf.num_dofs)
emf.resize(num_dofs);
@@ -74,4 +76,6 @@ electromagnetic_field_gpu::copyout(electromagnetic_field& emf)
Hx.copyout(emf.Hx.data());
Hy.copyout(emf.Hy.data());
Hz.copyout(emf.Hz.data());
-}
\ No newline at end of file
+}
+
+} // namespace maxwell
diff --git a/src/maxwell_solver.cpp b/src/maxwell_solver.cpp
index ca3b15b78c7cc7012948a9a623131214b0ca3350..5f65394c1be8422eb38eaf2144d55fe8e3180ac8 100644
--- a/src/maxwell_solver.cpp
+++ b/src/maxwell_solver.cpp
@@ -1,10 +1,6 @@
-#include <iostream>
-
-#include "gmsh_io.h"
-#include "maxwell_interface.h"
-#include "kernels_cpu.h"
-#include "kernels_gpu.h"
+#include "gmsh.h"
#include "silo_output.hpp"
+#include "maxwell_interface.h"
static void
make_geometry(int order, double mesh_h)
@@ -32,267 +28,6 @@ make_geometry(int order, double mesh_h)
gmm::setSize(vp, mesh_h);
}
-struct solver_state
-{
- std::vector<entity_data_cpu> eds;
- electromagnetic_field emf_curr;
- electromagnetic_field emf_next;
- electromagnetic_material_params matparams;
-
- double delta_t;
- double curr_time;
- size_t curr_timestep;
-
- electromagnetic_field dx;
- electromagnetic_field dy;
- electromagnetic_field dz;
-};
-
-struct solver_state_gpu
-{
- std::vector<entity_data_cpu> eds;
- std::vector<entity_data_gpu> edgs;
- electromagnetic_field_gpu emf_curr;
- electromagnetic_field_gpu emf_next;
- electromagnetic_material_params matparams;
-
- double delta_t;
- double curr_time;
- size_t curr_timestep;
-
- electromagnetic_field_gpu dx;
- electromagnetic_field_gpu dy;
- electromagnetic_field_gpu dz;
-};
-
-void init_from_model(const model& mod, solver_state& state)
-{
- state.emf_curr.resize( mod.num_dofs() );
- state.emf_next.resize( mod.num_dofs() );
- state.dx.resize( mod.num_dofs() );
- state.dy.resize( mod.num_dofs() );
- state.dz.resize( mod.num_dofs() );
-
- auto E = [](const point_3d& pt) -> vec3d {
- vec3d ret;
- ret(0) = std::sin(M_PI * pt.y());
- ret(1) = 0;
- ret(2) = 0;
- return ret;
- };
-
- for (auto& e : mod)
- {
- entity_data_cpu ed;
- e.populate_entity_data(ed);
- state.eds.push_back( std::move(ed) );
- e.project(E, state.emf_curr.Ex, state.emf_curr.Ey, state.emf_curr.Ez);
- }
-}
-
-void init_from_model(const model& mod, solver_state_gpu& state)
-{
- state.emf_curr.resize( mod.num_dofs() );
- state.emf_next.resize( mod.num_dofs() );
- state.dx.resize( mod.num_dofs() );
- state.dy.resize( mod.num_dofs() );
- state.dz.resize( mod.num_dofs() );
-
- electromagnetic_field emf_io;
- emf_io.resize( mod.num_dofs() );
-
- auto E = [](const point_3d& pt) -> vec3d {
- vec3d ret;
- ret(0) = std::sin(M_PI * pt.y());
- ret(1) = 0;
- ret(2) = 0;
- return ret;
- };
-
-
- for (auto& e : mod)
- {
- entity_data_cpu ed;
- e.populate_entity_data(ed);
- e.project(E, emf_io.Ex, emf_io.Ey, emf_io.Ez);
- entity_data_gpu edg(ed);
- state.eds.push_back( std::move(ed) );
- state.edgs.push_back( std::move(edg) );
- }
-
- state.emf_curr.copyin(emf_io);
-}
-
-void compute_curls(solver_state& state, const electromagnetic_field& curr,
- electromagnetic_field& next)
-{
- for (const auto& ed : state.eds)
- {
- compute_field_derivatives(ed, curr.Ex, state.dx.Ex, state.dy.Ex, state.dz.Ex);
- compute_field_derivatives(ed, curr.Ey, state.dx.Ey, state.dy.Ey, state.dz.Ey);
- compute_field_derivatives(ed, curr.Ez, state.dx.Ez, state.dy.Ez, state.dz.Ez);
- compute_field_derivatives(ed, curr.Hx, state.dx.Hx, state.dy.Hx, state.dz.Hx);
- compute_field_derivatives(ed, curr.Hy, state.dx.Hy, state.dy.Hy, state.dz.Hy);
- compute_field_derivatives(ed, curr.Hz, state.dx.Hz, state.dy.Hz, state.dz.Hz);
- }
-
- next.Ex = state.dy.Hz - state.dz.Hy;
- next.Ey = state.dz.Hx - state.dx.Hz;
- next.Ez = state.dx.Hy - state.dy.Hx;
- next.Hx = state.dz.Ey - state.dy.Ez;
- next.Hy = state.dx.Ez - state.dz.Ex;
- next.Hz = state.dy.Ex - state.dx.Ey;
-}
-
-void compute_curls(solver_state_gpu& state, const electromagnetic_field_gpu& curr,
- electromagnetic_field_gpu& next)
-{
- auto currp = curr.data();
- auto nextp = next.data();
- auto dxp = state.dx.data();
- auto dyp = state.dy.data();
- auto dzp = state.dz.data();
-
- for (const auto& edg : state.edgs)
- {
- gpu_compute_field_derivatives(edg, currp.Ex, dxp.Ex, dyp.Ex, dzp.Ex, 1.0);
- gpu_compute_field_derivatives(edg, currp.Ey, dxp.Ey, dyp.Ey, dzp.Ey, 1.0);
- gpu_compute_field_derivatives(edg, currp.Ez, dxp.Ez, dyp.Ez, dzp.Ez, 1.0);
- gpu_compute_field_derivatives(edg, currp.Hx, dxp.Hx, dyp.Hx, dzp.Hx, 1.0);
- gpu_compute_field_derivatives(edg, currp.Hy, dxp.Hy, dyp.Hy, dzp.Hy, 1.0);
- gpu_compute_field_derivatives(edg, currp.Hz, dxp.Hz, dyp.Hz, dzp.Hz, 1.0);
- }
-
-
- gpu_subtract(nextp.Ex, dyp.Hz, dzp.Hy, nextp.num_dofs);
- gpu_subtract(nextp.Ey, dzp.Hx, dxp.Hz, nextp.num_dofs);
- gpu_subtract(nextp.Ez, dxp.Hy, dyp.Hx, nextp.num_dofs);
- gpu_subtract(nextp.Hx, dzp.Ey, dyp.Ez, nextp.num_dofs);
- gpu_subtract(nextp.Hy, dxp.Ez, dzp.Ex, nextp.num_dofs);
- gpu_subtract(nextp.Hz, dyp.Ex, dxp.Ey, nextp.num_dofs);
-}
-
-void apply_operator(solver_state& state, const electromagnetic_field& curr,
- electromagnetic_field& next)
-{
- compute_curls(state, curr, next);
-}
-
-void apply_operator(solver_state_gpu& state, const electromagnetic_field_gpu& curr,
- electromagnetic_field_gpu& next)
-{
- compute_curls(state, curr, next);
-}
-
-void
-export_to_silo(const model& mod, const solver_state& state, const std::string& filename)
-{
- silo sdb;
- sdb.create_db(filename);
- sdb.import_mesh_from_gmsh();
- sdb.write_mesh();
-
- std::vector<double> curr_Ex; curr_Ex.resize( mod.num_cells() );
- std::vector<double> curr_Ey; curr_Ey.resize( mod.num_cells() );
- std::vector<double> curr_Ez; curr_Ez.resize( mod.num_cells() );
- std::vector<double> curr_Hx; curr_Hx.resize( mod.num_cells() );
- std::vector<double> curr_Hy; curr_Hy.resize( mod.num_cells() );
- std::vector<double> curr_Hz; curr_Hz.resize( mod.num_cells() );
-
- std::vector<double> next_Ex; next_Ex.resize( mod.num_cells() );
- std::vector<double> next_Ey; next_Ey.resize( mod.num_cells() );
- std::vector<double> next_Ez; next_Ez.resize( mod.num_cells() );
- std::vector<double> next_Hx; next_Hx.resize( mod.num_cells() );
- std::vector<double> next_Hy; next_Hy.resize( mod.num_cells() );
- std::vector<double> next_Hz; next_Hz.resize( mod.num_cells() );
-
- for (size_t iE = 0; iE < mod.num_entities(); iE++)
- {
- auto& e = mod.entity_at(iE);
- for (size_t iT = 0; iT < e.num_cells(); iT++)
- {
- auto& pe = e.cell(iT);
- auto& re = e.cell_refelem(pe);
- vecxd phi_bar = re.basis_functions({1./3., 1./3., 1./3.});
- auto num_bf = re.num_basis_functions();
- auto ofs = e.cell_global_dof_offset(iT);
- auto gi = e.cell_global_index_by_gmsh(iT);
- curr_Ex[gi] = state.emf_curr.Ex.segment(ofs, num_bf).dot(phi_bar);
- curr_Ey[gi] = state.emf_curr.Ey.segment(ofs, num_bf).dot(phi_bar);
- curr_Ez[gi] = state.emf_curr.Ez.segment(ofs, num_bf).dot(phi_bar);
- curr_Hx[gi] = state.emf_curr.Hx.segment(ofs, num_bf).dot(phi_bar);
- curr_Hy[gi] = state.emf_curr.Hy.segment(ofs, num_bf).dot(phi_bar);
- curr_Hz[gi] = state.emf_curr.Hz.segment(ofs, num_bf).dot(phi_bar);
- next_Ex[gi] = state.emf_next.Ex.segment(ofs, num_bf).dot(phi_bar);
- next_Ey[gi] = state.emf_next.Ey.segment(ofs, num_bf).dot(phi_bar);
- next_Ez[gi] = state.emf_next.Ez.segment(ofs, num_bf).dot(phi_bar);
- next_Hx[gi] = state.emf_next.Hx.segment(ofs, num_bf).dot(phi_bar);
- next_Hy[gi] = state.emf_next.Hy.segment(ofs, num_bf).dot(phi_bar);
- next_Hz[gi] = state.emf_next.Hz.segment(ofs, num_bf).dot(phi_bar);
- }
- }
-
- sdb.write_zonal_variable("curr_Ex", curr_Ex);
- sdb.write_zonal_variable("curr_Ey", curr_Ey);
- sdb.write_zonal_variable("curr_Ez", curr_Ez);
- sdb.write_zonal_variable("curr_Hx", curr_Hx);
- sdb.write_zonal_variable("curr_Hy", curr_Hy);
- sdb.write_zonal_variable("curr_Hz", curr_Hz);
- //sdb.write_vector_expression("curr_E", "{curr_Ex, curr_Ey, curr_Ez}");
- //sdb.write_vector_expression("curr_H", "{curr_Hx, curr_Hy, curr_Hz}");
- sdb.write_zonal_variable("next_Ex", next_Ex);
- sdb.write_zonal_variable("next_Ey", next_Ey);
- sdb.write_zonal_variable("next_Ez", next_Ez);
- sdb.write_zonal_variable("next_Hx", next_Hx);
- sdb.write_zonal_variable("next_Hy", next_Hy);
- sdb.write_zonal_variable("next_Hz", next_Hz);
- //sdb.write_vector_expression("next_E", "{next_Ex, next_Ey, next_Ez}");
- //sdb.write_vector_expression("next_H", "{next_Hx, next_Hy, next_Hz}");
-}
-
-void
-export_to_silo(const model& mod, const electromagnetic_field& emf,
- const std::string& filename)
-{
- silo sdb;
- sdb.create_db(filename);
- sdb.import_mesh_from_gmsh();
- sdb.write_mesh();
-
- std::vector<double> Ex; Ex.resize( mod.num_cells() );
- std::vector<double> Ey; Ey.resize( mod.num_cells() );
- std::vector<double> Ez; Ez.resize( mod.num_cells() );
- std::vector<double> Hx; Hx.resize( mod.num_cells() );
- std::vector<double> Hy; Hy.resize( mod.num_cells() );
- std::vector<double> Hz; Hz.resize( mod.num_cells() );
-
- for (size_t iE = 0; iE < mod.num_entities(); iE++)
- {
- auto& e = mod.entity_at(iE);
- for (size_t iT = 0; iT < e.num_cells(); iT++)
- {
- auto& pe = e.cell(iT);
- auto& re = e.cell_refelem(pe);
- vecxd phi_bar = re.basis_functions({1./3., 1./3., 1./3.});
- auto num_bf = re.num_basis_functions();
- auto ofs = e.cell_global_dof_offset(iT);
- auto gi = e.cell_global_index_by_gmsh(iT);
- Ex[gi] = emf.Ex.segment(ofs, num_bf).dot(phi_bar);
- Ey[gi] = emf.Ey.segment(ofs, num_bf).dot(phi_bar);
- Ez[gi] = emf.Ez.segment(ofs, num_bf).dot(phi_bar);
- Hx[gi] = emf.Hx.segment(ofs, num_bf).dot(phi_bar);
- Hy[gi] = emf.Hy.segment(ofs, num_bf).dot(phi_bar);
- Hz[gi] = emf.Hz.segment(ofs, num_bf).dot(phi_bar);
- }
- }
-
- sdb.write_zonal_variable("Ex", Ex);
- sdb.write_zonal_variable("Ey", Ey);
- sdb.write_zonal_variable("Ez", Ez);
- sdb.write_zonal_variable("Hx", Hx);
- sdb.write_zonal_variable("Hy", Hy);
- sdb.write_zonal_variable("Hz", Hz);
-}
int main(int argc, const char *argv[])
{
@@ -301,15 +36,15 @@ int main(int argc, const char *argv[])
make_geometry(1, 0.05);
model mod(1,2);
- solver_state_gpu state;
- init_from_model(mod, state);
-
- apply_operator(state, state.emf_curr, state.emf_next);
+
+ maxwell::solver_state_gpu state;
+ maxwell::init_from_model(mod, state);
+ maxwell::apply_operator(state, state.emf_curr, state.emf_next);
- electromagnetic_field res;
+ maxwell::field res;
state.emf_next.copyout(res);
- export_to_silo(mod, res, "test.silo");
+ maxwell::export_to_silo(mod, res, "test.silo");
//gmsh::finalize();