Skip to content
Snippets Groups Projects
Commit b4fcecb0 authored by Matteo Cicuttin's avatar Matteo Cicuttin
Browse files

Put maxwell stuff in namespace.

parent a4cb746a
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
CMakeLists.txt
+ 2
− 0
View file @ b4fcecb0
......@@ -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})
......
include/maxwell_interface.h
+ 53
− 7
View file @ b4fcecb0
#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
src/maxwell_interface.cpp
+ 12
− 8
View file @ b4fcecb0
#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
src/maxwell_solver.cpp
+ 8
− 273
View file @ b4fcecb0
#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();
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment