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

Volumetric sources.

parent 4b2bf8ef
No related branches found
No related tags found
No related merge requests found
...@@ -65,6 +65,7 @@ public: ...@@ -65,6 +65,7 @@ public:
//matxd project(const vector_function&) const; //matxd project(const vector_function&) const;
void project(const scalar_function&, vecxd&) const; void project(const scalar_function&, vecxd&) const;
void project(const vector_function&, vecxd&, vecxd&, vecxd&) const; void project(const vector_function&, vecxd&, vecxd&, vecxd&) const;
void project(const vector_function&, double *, double *, double *) const;
matxd project_on_face(size_t iF, const vector_function&) const; matxd project_on_face(size_t iF, const vector_function&) const;
......
...@@ -199,6 +199,74 @@ public: ...@@ -199,6 +199,74 @@ public:
} }
}; };
template<typename T>
class pinned_vector
{
T * m_ptr;
size_t m_size;
public:
pinned_vector()
: m_ptr(nullptr), m_size(0)
{}
pinned_vector(size_t size)
: m_ptr(nullptr), m_size(0)
{
(void)checkGPU( gpuMallocHost((void**)&m_ptr, size*sizeof(T)) );
memset(m_ptr, '\0', size*sizeof(double));
m_size = size;
}
pinned_vector(const pinned_vector&) = delete;
pinned_vector(pinned_vector&&) = delete;
pinned_vector& operator=(const pinned_vector&) = delete;
pinned_vector& operator=(pinned_vector&&) = delete;
void resize(size_t size)
{
if (size == m_size)
return;
if (m_ptr != nullptr)
(void) checkGPU(gpuFreeHost(m_ptr));
(void)checkGPU( gpuMallocHost((void**)&m_ptr, size*sizeof(T)) );
memset(m_ptr, '\0', size*sizeof(double));
m_size = size;
}
void zero(void)
{
memset(m_ptr, '\0', m_size*sizeof(double));
}
size_t size(void) const
{
return m_size;
}
T* data(void)
{
return m_ptr;
}
const T* data(void) const
{
return m_ptr;
}
T operator[](size_t pos) const
{
return m_ptr[pos];
}
T& operator[](size_t pos)
{
return m_ptr[pos];
}
};
template<typename T> template<typename T>
class device_vector class device_vector
{ {
...@@ -255,6 +323,15 @@ public: ...@@ -255,6 +323,15 @@ public:
(void)checkGPU( gpuMemcpyAsync(m_devptr, src, size*sizeof(T), gpuMemcpyHostToDevice, st) ); (void)checkGPU( gpuMemcpyAsync(m_devptr, src, size*sizeof(T), gpuMemcpyHostToDevice, st) );
} }
void copyin(const T *src, size_t start, size_t num, const stream& st)
{
assert(start < m_size);
assert(start+num <= m_size);
const T *c_begin = src + start;
T *d_begin = m_devptr + start;
(void)checkGPU( gpuMemcpyAsync(d_begin, c_begin, num*sizeof(T), gpuMemcpyHostToDevice, st) );
}
void copyout(T *dst) const void copyout(T *dst) const
{ {
if (m_devptr == nullptr) if (m_devptr == nullptr)
......
...@@ -148,9 +148,12 @@ gpu_compute_field_derivatives(const entity_data_gpu& edg, const double *F, ...@@ -148,9 +148,12 @@ gpu_compute_field_derivatives(const entity_data_gpu& edg, const double *F,
void gpu_compute_flux_lifting(entity_data_gpu&, const double *, double *, void gpu_compute_flux_lifting(entity_data_gpu&, const double *, double *,
gpuStream_t stream = 0); gpuStream_t stream = 0);
void gpu_subtract(double *, const double *, const double *, size_t, void gpu_curl(double *, const double *, const double *, size_t,
gpuStream_t stream = 0); gpuStream_t stream = 0);
void gpu_curl(double *, const double *, const double *, const double *,
size_t, gpuStream_t stream = 0);
void void
gpu_compute_euler_update(double *out, const double *y, const double *k, gpu_compute_euler_update(double *out, const double *y, const double *k,
const double *reaction, const double *material, size_t num_dofs, double dt, const double *reaction, const double *material, size_t num_dofs, double dt,
......
...@@ -65,6 +65,8 @@ public: ...@@ -65,6 +65,8 @@ public:
bool initial_Hfield_defined(void) const; bool initial_Hfield_defined(void) const;
vec3d initial_Hfield(const point_3d&) const; vec3d initial_Hfield(const point_3d&) const;
bool volume_has_source(int) const;
vec3d eval_volume_source(int, const point_3d&, double) const;
vec3d eval_boundary_source(int, const point_3d&, double) const; vec3d eval_boundary_source(int, const point_3d&, double) const;
vec3d eval_interface_source(int, const point_3d&, double) const; vec3d eval_interface_source(int, const point_3d&, double) const;
...@@ -250,6 +252,9 @@ struct solver_state ...@@ -250,6 +252,9 @@ struct solver_state
field k4; field k4;
field tmp; field tmp;
vecxd Jx_src;
vecxd Jy_src;
vecxd Jz_src;
field bndsrcs; field bndsrcs;
}; };
...@@ -279,6 +284,17 @@ struct solver_state_gpu ...@@ -279,6 +284,17 @@ struct solver_state_gpu
field_gpu k4; field_gpu k4;
field_gpu tmp; field_gpu tmp;
pinned_vector<double> Jx_src;
pinned_vector<double> Jy_src;
pinned_vector<double> Jz_src;
device_vector<double> Jx_src_gpu;
device_vector<double> Jy_src_gpu;
device_vector<double> Jz_src_gpu;
device_vector<double> Jx_src_gpu_buf;
device_vector<double> Jy_src_gpu_buf;
device_vector<double> Jz_src_gpu_buf;
field_gpu bndsrcs; field_gpu bndsrcs;
field_gpu bndsrcs_buf; field_gpu bndsrcs_buf;
......
...@@ -166,6 +166,51 @@ entity::project(const vector_function& function, vecxd& pfx, vecxd& pfy, vecxd& ...@@ -166,6 +166,51 @@ entity::project(const vector_function& function, vecxd& pfx, vecxd& pfy, vecxd&
} }
} }
void
entity::project(const vector_function& function, double *pfx, double *pfy, double *pfz) const
{
auto num_bf = reference_cells[0].num_basis_functions();
for (size_t iT = 0; iT < physical_cells.size(); iT++)
{
const auto& pe = physical_cells[iT];
assert( pe.orientation() < reference_cells.size() );
const auto& re = reference_cells[pe.orientation()];
const auto pqps = pe.integration_points();
const auto dets = pe.determinants();
assert(pqps.size() == dets.size());
matxd mass = matxd::Zero(num_bf, num_bf);
vecxd rhs_x = vecxd::Zero(num_bf);
vecxd rhs_y = vecxd::Zero(num_bf);
vecxd rhs_z = vecxd::Zero(num_bf);
for (size_t iQp = 0; iQp < pqps.size(); iQp++)
{
const auto& pqp = pqps[iQp];
const vecxd phi = re.basis_functions(iQp);
mass += pqp.weight() * phi * phi.transpose();
vec3d fval = function(pqp.point());
rhs_x += pqp.weight() * fval(0) * phi;
rhs_y += pqp.weight() * fval(1) * phi;
rhs_z += pqp.weight() * fval(2) * phi;
}
Eigen::LDLT<matxd> mass_ldlt;
mass_ldlt.compute(mass);
vecxd p_x = mass_ldlt.solve(rhs_x);
vecxd p_y = mass_ldlt.solve(rhs_y);
vecxd p_z = mass_ldlt.solve(rhs_z);
for (size_t i = 0; i < num_bf; i++)
{
pfx[m_dof_base + iT*num_bf + i] = p_x[i];
pfy[m_dof_base + iT*num_bf + i] = p_y[i];
pfz[m_dof_base + iT*num_bf + i] = p_z[i];
}
}
}
matxd matxd
entity::project_on_face(size_t iF, const vector_function& function) const entity::project_on_face(size_t iF, const vector_function& function) const
{ {
......
...@@ -292,7 +292,7 @@ gpu_compute_flux_lifting(entity_data_gpu& edg, const double *f, double *out, ...@@ -292,7 +292,7 @@ gpu_compute_flux_lifting(entity_data_gpu& edg, const double *f, double *out,
} }
void __global__ void __global__
gpu_subtract_kernel(double * __restrict__ dst, const double * __restrict__ add, gpu_curl_kernel(double * __restrict__ dst, const double * __restrict__ add,
const double * __restrict__ sub, size_t num_dofs) const double * __restrict__ sub, size_t num_dofs)
{ {
int32_t dof = blockIdx.x * blockDim.x + threadIdx.x; int32_t dof = blockIdx.x * blockDim.x + threadIdx.x;
...@@ -301,7 +301,7 @@ gpu_subtract_kernel(double * __restrict__ dst, const double * __restrict__ add, ...@@ -301,7 +301,7 @@ gpu_subtract_kernel(double * __restrict__ dst, const double * __restrict__ add,
} }
void void
gpu_subtract(double *dst, const double *add, const double *sub, size_t num_dofs, gpu_curl(double *dst, const double *add, const double *sub, size_t num_dofs,
cudaStream_t stream) cudaStream_t stream)
{ {
static const size_t SUBTRACT_THREADS = 256; static const size_t SUBTRACT_THREADS = 256;
...@@ -309,7 +309,29 @@ gpu_subtract(double *dst, const double *add, const double *sub, size_t num_dofs, ...@@ -309,7 +309,29 @@ gpu_subtract(double *dst, const double *add, const double *sub, size_t num_dofs,
if (num_dofs % SUBTRACT_THREADS != 0) if (num_dofs % SUBTRACT_THREADS != 0)
gs += 1; gs += 1;
gpu_subtract_kernel<<<gs, SUBTRACT_THREADS, 0, stream>>>(dst, add, sub, num_dofs); gpu_curl_kernel<<<gs, SUBTRACT_THREADS, 0, stream>>>(dst, add, sub, num_dofs);
}
void __global__
gpu_curl_kernel(double * __restrict__ dst, const double * __restrict__ add,
const double * __restrict__ sub, const double * __restrict__ source,
size_t num_dofs)
{
int32_t dof = blockIdx.x * blockDim.x + threadIdx.x;
if (dof < num_dofs)
dst[dof] = add[dof] - sub[dof] - source[dof];
}
void
gpu_curl(double *dst, const double *add, const double *sub, const double *source,
size_t num_dofs, cudaStream_t stream)
{
static const size_t SUBTRACT_THREADS = 256;
auto gs = num_dofs/SUBTRACT_THREADS;
if (num_dofs % SUBTRACT_THREADS != 0)
gs += 1;
gpu_curl_kernel<<<gs, SUBTRACT_THREADS, 0, stream>>>(dst, add, sub, source, num_dofs);
} }
......
...@@ -71,6 +71,9 @@ init_from_model(const model& mod, solver_state& state) ...@@ -71,6 +71,9 @@ init_from_model(const model& mod, solver_state& state)
state.jumps.resize( mod.num_fluxes() ); state.jumps.resize( mod.num_fluxes() );
state.fluxes.resize( mod.num_fluxes() ); state.fluxes.resize( mod.num_fluxes() );
state.Jx_src = vecxd::Zero( mod.num_dofs() );
state.Jy_src = vecxd::Zero( mod.num_dofs() );
state.Jz_src = vecxd::Zero( mod.num_dofs() );
state.bndsrcs.resize( mod.num_fluxes() ); state.bndsrcs.resize( mod.num_fluxes() );
for (auto& e : mod) for (auto& e : mod)
...@@ -97,9 +100,9 @@ compute_curls(solver_state& state, const field& curr, field& next) ...@@ -97,9 +100,9 @@ compute_curls(solver_state& state, const field& curr, field& next)
compute_field_derivatives(ed, curr.Hz, state.dx.Hz, state.dy.Hz, state.dz.Hz); compute_field_derivatives(ed, curr.Hz, state.dx.Hz, state.dy.Hz, state.dz.Hz);
} }
next.Ex = state.dy.Hz - state.dz.Hy; next.Ex = state.dy.Hz - state.dz.Hy - state.Jx_src;
next.Ey = state.dz.Hx - state.dx.Hz; next.Ey = state.dz.Hx - state.dx.Hz - state.Jy_src;
next.Ez = state.dx.Hy - state.dy.Hx; next.Ez = state.dx.Hy - state.dy.Hx - state.Jz_src;
next.Hx = state.dz.Ey - state.dy.Ez; next.Hx = state.dz.Ey - state.dy.Ez;
next.Hy = state.dx.Ez - state.dz.Ex; next.Hy = state.dx.Ez - state.dz.Ex;
next.Hz = state.dy.Ex - state.dx.Ey; next.Hz = state.dy.Ex - state.dx.Ey;
...@@ -130,9 +133,9 @@ compute_curls_H(solver_state& state, const field& curr, field& next) ...@@ -130,9 +133,9 @@ compute_curls_H(solver_state& state, const field& curr, field& next)
compute_field_derivatives(ed, curr.Hz, state.dx.Hz, state.dy.Hz, state.dz.Hz); compute_field_derivatives(ed, curr.Hz, state.dx.Hz, state.dy.Hz, state.dz.Hz);
} }
next.Ex = state.dy.Hz - state.dz.Hy; next.Ex = state.dy.Hz - state.dz.Hy - state.Jx_src;
next.Ey = state.dz.Hx - state.dx.Hz; next.Ey = state.dz.Hx - state.dx.Hz - state.Jy_src;
next.Ez = state.dx.Hy - state.dy.Hx; next.Ez = state.dx.Hy - state.dy.Hx - state.Jz_src;
} }
...@@ -536,20 +539,42 @@ timestep(solver_state& state, time_integrator_type ti) ...@@ -536,20 +539,42 @@ timestep(solver_state& state, time_integrator_type ti)
} }
void void
do_sources(const model& mod, maxwell::solver_state& state, eval_volume_sources(const model& mod, const parameter_loader& mpl, solver_state& state)
maxwell::parameter_loader& mpl) {
for (auto& e : mod)
{
auto etag = e.gmsh_tag();
if ( not mpl.volume_has_source(etag) )
continue;
auto fJ = [&](const point_3d& pt) -> vec3d {
return mpl.eval_volume_source(etag, pt, state.curr_time);
};
e.project(fJ, state.Jx_src, state.Jy_src, state.Jz_src);
}
}
void
do_sources(const model& mod, solver_state& state, parameter_loader& mpl)
{ {
if ( mpl.source_has_changed_state() ) if ( mpl.source_has_changed_state() )
{ {
state.Jx_src = vecxd::Zero( state.Jx_src.size() );
state.Jy_src = vecxd::Zero( state.Jy_src.size() );
state.Jz_src = vecxd::Zero( state.Jz_src.size() );
state.bndsrcs.zero(); state.bndsrcs.zero();
mpl.source_was_cleared(); mpl.source_was_cleared();
} }
if ( mpl.volume_sources_enabled() )
eval_volume_sources(mod, mpl, state);
if ( mpl.boundary_sources_enabled() ) if ( mpl.boundary_sources_enabled() )
maxwell::eval_boundary_sources(mod, mpl, state, state.bndsrcs); eval_boundary_sources(mod, mpl, state, state.bndsrcs);
if ( mpl.interface_sources_enabled() ) if ( mpl.interface_sources_enabled() )
maxwell::eval_interface_sources(mod, mpl, state, state.bndsrcs); eval_interface_sources(mod, mpl, state, state.bndsrcs);
} }
void void
......
...@@ -79,6 +79,16 @@ void init_from_model(const model& mod, solver_state_gpu& state) ...@@ -79,6 +79,16 @@ void init_from_model(const model& mod, solver_state_gpu& state)
state.k4.resize( mod.num_dofs() ); state.k4.resize( mod.num_dofs() );
state.tmp.resize( mod.num_dofs() ); state.tmp.resize( mod.num_dofs() );
state.Jx_src.resize( mod.num_dofs() );
state.Jy_src.resize( mod.num_dofs() );
state.Jz_src.resize( mod.num_dofs() );
state.Jx_src_gpu.resize( mod.num_dofs() );
state.Jy_src_gpu.resize( mod.num_dofs() );
state.Jz_src_gpu.resize( mod.num_dofs() );
state.Jx_src_gpu_buf.resize( mod.num_dofs() );
state.Jy_src_gpu_buf.resize( mod.num_dofs() );
state.Jz_src_gpu_buf.resize( mod.num_dofs() );
state.bndsrcs.resize( mod.num_fluxes() ); state.bndsrcs.resize( mod.num_fluxes() );
state.bndsrcs_decomp_cpu.resize( mod.num_fluxes() ); state.bndsrcs_decomp_cpu.resize( mod.num_fluxes() );
...@@ -158,13 +168,16 @@ compute_curls(solver_state_gpu& state, const field_gpu& curr, field_gpu& next) ...@@ -158,13 +168,16 @@ compute_curls(solver_state_gpu& state, const field_gpu& curr, field_gpu& next)
gpu_compute_field_derivatives(edg, currp.Hz, dxp.Hz, dyp.Hz, dzp.Hz, 1.0, state.compute_stream); gpu_compute_field_derivatives(edg, currp.Hz, dxp.Hz, dyp.Hz, dzp.Hz, 1.0, state.compute_stream);
} }
auto Jx = state.Jx_src_gpu.data();
auto Jy = state.Jy_src_gpu.data();
auto Jz = state.Jz_src_gpu.data();
gpu_subtract(nextp.Ex, dyp.Hz, dzp.Hy, nextp.num_dofs, state.compute_stream); gpu_curl(nextp.Ex, dyp.Hz, dzp.Hy, Jx, nextp.num_dofs, state.compute_stream);
gpu_subtract(nextp.Ey, dzp.Hx, dxp.Hz, nextp.num_dofs, state.compute_stream); gpu_curl(nextp.Ey, dzp.Hx, dxp.Hz, Jy, nextp.num_dofs, state.compute_stream);
gpu_subtract(nextp.Ez, dxp.Hy, dyp.Hx, nextp.num_dofs, state.compute_stream); gpu_curl(nextp.Ez, dxp.Hy, dyp.Hx, Jz, nextp.num_dofs, state.compute_stream);
gpu_subtract(nextp.Hx, dzp.Ey, dyp.Ez, nextp.num_dofs, state.compute_stream); gpu_curl(nextp.Hx, dzp.Ey, dyp.Ez, nextp.num_dofs, state.compute_stream);
gpu_subtract(nextp.Hy, dxp.Ez, dzp.Ex, nextp.num_dofs, state.compute_stream); gpu_curl(nextp.Hy, dxp.Ez, dzp.Ex, nextp.num_dofs, state.compute_stream);
gpu_subtract(nextp.Hz, dyp.Ex, dxp.Ey, nextp.num_dofs, state.compute_stream); gpu_curl(nextp.Hz, dyp.Ex, dxp.Ey, nextp.num_dofs, state.compute_stream);
} }
static void static void
...@@ -184,9 +197,9 @@ compute_curls_E(solver_state_gpu& state, const field_gpu& curr, field_gpu& next) ...@@ -184,9 +197,9 @@ compute_curls_E(solver_state_gpu& state, const field_gpu& curr, field_gpu& next)
} }
gpu_subtract(nextp.Hx, dzp.Ey, dyp.Ez, nextp.num_dofs, state.compute_stream); gpu_curl(nextp.Hx, dzp.Ey, dyp.Ez, nextp.num_dofs, state.compute_stream);
gpu_subtract(nextp.Hy, dxp.Ez, dzp.Ex, nextp.num_dofs, state.compute_stream); gpu_curl(nextp.Hy, dxp.Ez, dzp.Ex, nextp.num_dofs, state.compute_stream);
gpu_subtract(nextp.Hz, dyp.Ex, dxp.Ey, nextp.num_dofs, state.compute_stream); gpu_curl(nextp.Hz, dyp.Ex, dxp.Ey, nextp.num_dofs, state.compute_stream);
} }
static void static void
...@@ -205,10 +218,13 @@ compute_curls_H(solver_state_gpu& state, const field_gpu& curr, field_gpu& next) ...@@ -205,10 +218,13 @@ compute_curls_H(solver_state_gpu& state, const field_gpu& curr, field_gpu& next)
gpu_compute_field_derivatives(edg, currp.Hz, dxp.Hz, dyp.Hz, dzp.Hz, 1.0, state.compute_stream); gpu_compute_field_derivatives(edg, currp.Hz, dxp.Hz, dyp.Hz, dzp.Hz, 1.0, state.compute_stream);
} }
auto Jx = state.Jx_src_gpu.data();
auto Jy = state.Jy_src_gpu.data();
auto Jz = state.Jz_src_gpu.data();
gpu_subtract(nextp.Ex, dyp.Hz, dzp.Hy, nextp.num_dofs, state.compute_stream); gpu_curl(nextp.Ex, dyp.Hz, dzp.Hy, Jx, nextp.num_dofs, state.compute_stream);
gpu_subtract(nextp.Ey, dzp.Hx, dxp.Hz, nextp.num_dofs, state.compute_stream); gpu_curl(nextp.Ey, dzp.Hx, dxp.Hz, Jy, nextp.num_dofs, state.compute_stream);
gpu_subtract(nextp.Ez, dxp.Hy, dyp.Hx, nextp.num_dofs, state.compute_stream); gpu_curl(nextp.Ez, dxp.Hy, dyp.Hx, Jz, nextp.num_dofs, state.compute_stream);
} }
static void static void
...@@ -426,12 +442,47 @@ prepare_sources(const model& mod, maxwell::solver_state_gpu& state, ...@@ -426,12 +442,47 @@ prepare_sources(const model& mod, maxwell::solver_state_gpu& state,
decompress_bndsrc(state, state.bndsrcs_buf, state.bndsrcs); decompress_bndsrc(state, state.bndsrcs_buf, state.bndsrcs);
} }
static void
copyin_volume_sources(const entity& e, solver_state_gpu& state)
{
auto ofs = e.dof_base();
auto num = e.num_dofs();
state.Jx_src_gpu_buf.copyin( state.Jx_src.data(), ofs, num, state.memcpy_stream );
state.Jy_src_gpu_buf.copyin( state.Jy_src.data(), ofs, num, state.memcpy_stream );
state.Jz_src_gpu_buf.copyin( state.Jz_src.data(), ofs, num, state.memcpy_stream );
}
static void
eval_volume_sources(const model& mod, const parameter_loader& mpl, solver_state_gpu& state)
{
for (auto& e : mod)
{
auto etag = e.gmsh_tag();
if ( not mpl.volume_has_source(etag) )
continue;
auto fJ = [&](const point_3d& pt) -> vec3d {
return mpl.eval_volume_source(etag, pt, state.curr_time);
};
e.project(fJ, state.Jx_src.data(), state.Jy_src.data(), state.Jz_src.data());
copyin_volume_sources(e, state);
}
}
void void
do_sources(const model& mod, maxwell::solver_state_gpu& state, do_sources(const model& mod, maxwell::solver_state_gpu& state,
maxwell::parameter_loader& mpl) maxwell::parameter_loader& mpl)
{ {
if ( mpl.source_has_changed_state() ) if ( mpl.source_has_changed_state() )
{ {
state.Jx_src.zero();
state.Jy_src.zero();
state.Jz_src.zero();
state.Jx_src_gpu.zero();
state.Jy_src_gpu.zero();
state.Jz_src_gpu.zero();
state.bndsrcs_decomp_cpu.zero(); state.bndsrcs_decomp_cpu.zero();
state.bndsrcs_cpu.zero(); state.bndsrcs_cpu.zero();
mpl.source_was_cleared(); mpl.source_was_cleared();
...@@ -441,11 +492,15 @@ do_sources(const model& mod, maxwell::solver_state_gpu& state, ...@@ -441,11 +492,15 @@ do_sources(const model& mod, maxwell::solver_state_gpu& state,
auto ie = mpl.interface_sources_enabled(); auto ie = mpl.interface_sources_enabled();
auto ve = mpl.volume_sources_enabled(); auto ve = mpl.volume_sources_enabled();
if ( ve )
eval_volume_sources(mod, mpl, state);
if ( be ) if ( be )
maxwell::eval_boundary_sources(mod, mpl, state, state.bndsrcs_decomp_cpu); eval_boundary_sources(mod, mpl, state, state.bndsrcs_decomp_cpu);
if ( ie ) if ( ie )
maxwell::eval_interface_sources(mod, mpl, state, state.bndsrcs_decomp_cpu); eval_interface_sources(mod, mpl, state, state.bndsrcs_decomp_cpu);
if ( be or ie or ve ) if ( be or ie or ve )
{ {
...@@ -458,6 +513,9 @@ void ...@@ -458,6 +513,9 @@ void
swap(maxwell::solver_state_gpu& state) swap(maxwell::solver_state_gpu& state)
{ {
checkGPU( gpuDeviceSynchronize() ); checkGPU( gpuDeviceSynchronize() );
std::swap(state.Jx_src_gpu, state.Jx_src_gpu_buf);
std::swap(state.Jy_src_gpu, state.Jy_src_gpu_buf);
std::swap(state.Jz_src_gpu, state.Jz_src_gpu_buf);
decompress_bndsrc(state, state.bndsrcs_buf, state.bndsrcs); decompress_bndsrc(state, state.bndsrcs_buf, state.bndsrcs);
std::swap(state.emf_curr, state.emf_next); std::swap(state.emf_curr, state.emf_next);
} }
......
...@@ -237,7 +237,7 @@ material_params_gpu::copyin(const material_params& mp) ...@@ -237,7 +237,7 @@ material_params_gpu::copyin(const material_params& mp)
#define SEC_VSRCS "sources"
#define SEC_BCONDS "bndconds" #define SEC_BCONDS "bndconds"
#define BC_KIND "kind" #define BC_KIND "kind"
#define BC_SOURCE "source" #define BC_SOURCE "source"
...@@ -277,6 +277,7 @@ parameter_loader::parameter_loader() ...@@ -277,6 +277,7 @@ parameter_loader::parameter_loader()
lua["postpro"]["H"]["silo_mode"] = "nodal"; lua["postpro"]["H"]["silo_mode"] = "nodal";
lua["postpro"]["J"] = lua.create_table(); lua["postpro"]["J"] = lua.create_table();
lua["postpro"]["J"]["silo_mode"] = "nodal"; lua["postpro"]["J"]["silo_mode"] = "nodal";
lua[SEC_VSRCS] = lua.create_table();
} }
bool bool
...@@ -559,6 +560,22 @@ parameter_loader::interface_type(int tag) const ...@@ -559,6 +560,22 @@ parameter_loader::interface_type(int tag) const
return interface_condition::UNSPECIFIED; return interface_condition::UNSPECIFIED;
} }
bool
parameter_loader::volume_has_source(int tag) const
{
auto src = lua[SEC_VSRCS][tag];
return src.valid();
}
vec3d
parameter_loader::eval_volume_source(int tag, const point_3d& pt, double t) const
{
vec3d ret;
sol::tie(ret(0), ret(1), ret(2)) =
lua[SEC_VSRCS][tag](tag, pt.x(), pt.y(), pt.z(), t);
return ret;
}
vec3d vec3d
parameter_loader::eval_boundary_source(int tag, const point_3d& pt, double t) const parameter_loader::eval_boundary_source(int tag, const point_3d& pt, double t) const
{ {
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment