diff --git a/doc/lua_api.md b/doc/lua_api.md
index 46e30ab76b3f3488e4d887d9bc938071d5013c9b..6ab8f7fc6bf68dd5556d406d2d0935e399f05a4f 100644
--- a/doc/lua_api.md
+++ b/doc/lua_api.md
@@ -144,7 +144,7 @@ None.
### Debug and validation
The solver has some facilities to allow validation and comparison with analytical solutions. Those facilities are accessed via the `debug` table. In the current implementation the `debug` table is undefined by default, so it must be initialized explicitly in the configuration with a statement `debug = {}`. Possible members of the `debug` table are:
-- `analytical_solution(tag, x, y, z, t)`: define the analytical solution of the problem. If defined, the numerical solution is compared with the analytical solution at each timestep in the whole domain. The function has to return 6 values, namely `Ex`, `Ey`, `Ez`, `Hx`, `Hy` and `Hz`.
+- `analytical_solution`: a function `f(tag, x, y, z, t)` defining the analytical solution of the problem. If defined, it can be used to compare the numerical solution with the analytical solution during the timestepping. The function has to return 6 values, namely `Ex`, `Ey`, `Ez`, `Hx`, `Hy` and `Hz`.
- `dump_cell_ranks` (bool): if true, and if the solver has MPI support compiled in, add to the Silo output a variable named `cell_ranks` showing the cell-to-rank mapping.
### Postprocessing
diff --git a/include/libgmshdg/gmsh_io.h b/include/libgmshdg/gmsh_io.h
index 629356b70cbdfa118ecbcc515ba1888864573de1..81d959b50f647796988a042dad898b3a5a9f8453 100644
--- a/include/libgmshdg/gmsh_io.h
+++ b/include/libgmshdg/gmsh_io.h
@@ -29,8 +29,6 @@ std::string quadrature_name(int);
std::string basis_func_name(int);
std::string basis_grad_name(int);
-using face_key = std::array<size_t, 3>;
-
enum class face_type : int
{
NONE = 0,
diff --git a/include/maxwell/maxwell_interface.h b/include/maxwell/maxwell_interface.h
index 00662a72f2b6ee8a79b3fdcb32da73d4091abbcf..0fb3af3d17f51c1aa6d579210f59ff83d7926c1a 100644
--- a/include/maxwell/maxwell_interface.h
+++ b/include/maxwell/maxwell_interface.h
@@ -435,4 +435,42 @@ void decompress_bndsrc(const solver_state_gpu& state, const field_gpu& csrcs,
#endif /* ENABLE_GPU_SOLVER */
+
+
+
+
+struct field_values {
+ double Ex;
+ double Ey;
+ double Ez;
+ double Hx;
+ double Hy;
+ double Hz;
+
+ field_values();
+ field_values(const vec3d&, const vec3d&);
+
+ field_values(double, double, double, double, double, double);
+
+ field_values& operator+=(const field_values&);
+ field_values operator+(const field_values&) const;
+ field_values& operator-=(const field_values&);
+ field_values operator-(const field_values&) const;
+ field_values& operator*=(const field_values&);
+ field_values operator*(const field_values&) const;
+ field_values& operator/=(const field_values&);
+ field_values operator/(const field_values&) const;
+
+ field_values& operator*=(double);
+ field_values operator*(double) const;
+ field_values& operator/=(double);
+ field_values operator/(double) const;
+};
+
+field_values operator*(double, const field_values&);
+std::ostream& operator<<(std::ostream&, const field_values&);
+
+field_values sqrt(const field_values&);
+double hsum(const field_values& m);
+
} // namespace maxwell
diff --git a/include/maxwell/maxwell_postpro.h b/include/maxwell/maxwell_postpro.h
new file mode 100644
index 0000000000000000000000000000000000000000..75e94cf99060336e194caa671bfa9e857760dbe6
--- /dev/null
+++ b/include/maxwell/maxwell_postpro.h
@@ -0,0 +1,11 @@
+#pragma once
+
+#include "maxwell/maxwell_interface.h"
+
+namespace maxwell {
+
+field_values eval_field(const field&, size_t, size_t, const vecxd&);
+field_values compute_error(const model&, const solver_state&, const parameter_loader&);
+field_values compute_energy(const model&, const solver_state&, const parameter_loader&);
+
+} // namespace maxwell
\ No newline at end of file
diff --git a/share/capacitor/capacitor.geo b/share/capacitor/capacitor.geo
deleted file mode 100644
index 9fd46c5cb2c3301f13cdf5a84d5953f30b4b94de..0000000000000000000000000000000000000000
--- a/share/capacitor/capacitor.geo
+++ /dev/null
@@ -1,15 +0,0 @@
-l = 0.02;
-d = 0.001;
-t = 0.0003;
-sh = 2*t+d;
-
-SetFactory("OpenCASCADE");
-Mesh.Algorithm3D = 10;
-
-Box(1) = {-l/2, -l/2, -sh/2, l, l, t};
-Box(2) = {-l/2, -l/2, -sh/2 + t, l, l, d};
-Box(3) = {-l/2, -l/2, -sh/2 + t+d, l, l, t};
-Sphere(4) = {0, 0, 0, l, -Pi/2, Pi/2, 2*Pi};
-Coherence;
-
-MeshSize{ PointsOf{ Volume{1,2,3}; } } = 0.0008;
diff --git a/share/capacitor/capacitor.m b/share/capacitor/capacitor.m
deleted file mode 100644
index f7abfbaf6ce246f4494f8b90ef456222bbfe0ee9..0000000000000000000000000000000000000000
--- a/share/capacitor/capacitor.m
+++ /dev/null
@@ -1,16 +0,0 @@
-l = 0.02;
-d = 0.001;
-epsilon = 10 * 8.85e-12;
-t = 1e-13;
-J = 1;
-
-A = l*l;
-I = J*A;
-
-C = epsilon*A/d;
-Q = I*t;
-
-V = Q/C
-E = V/d
-
-Q*d/(epsilon*l*A)
diff --git a/share/capacitor/params_capacitor.lua b/share/capacitor/params_capacitor.lua
deleted file mode 100644
index bfba7ec06d02a3b2a45a53d167983274389032e1..0000000000000000000000000000000000000000
--- a/share/capacitor/params_capacitor.lua
+++ /dev/null
@@ -1,66 +0,0 @@
-sim.name = "capacitor" -- simulation name
-sim.dt = 1e-13 -- timestep size
-sim.timesteps = 100000 -- num of iterations
-sim.gmsh_model = "capacitor.geo" -- gmsh model filename
-sim.use_gpu = 0 -- 0: cpu, 1: gpu
-sim.approx_order = 1 -- approximation order
-sim.time_integrator = "leapfrog"
-postpro.silo_output_rate = 10 -- rate at which to write silo files
-postpro.cycle_print_rate = 10 -- console print rate
-
-local alu = {1, 3}
-for i,v in ipairs(alu) do
- materials[v] = {}
- materials[v].epsilon = 1
- materials[v].mu = 1
- materials[v].sigma = 3.69e7
-end
-
-local diel = { 2 }
-for i,v in ipairs(diel) do
- materials[v] = {}
- materials[v].epsilon = 10
- materials[v].mu = 1
- materials[v].sigma = 0
-end
-
-local air = { 4 }
-for i,v in ipairs(air) do
- materials[v] = {}
- materials[v].epsilon = 1
- materials[v].mu = 1
- materials[v].sigma = 0
-end
-
--- ** Boundary conditions **
-local absorbing_bcs = {1}
-for i,v in ipairs(absorbing_bcs) do
- bndconds[v] = {}
- bndconds[v].kind = "impedance"
-end
-
-function interp(t, t0, t1, y0, y1)
- return (t-t0)*(y1-y0)/(t1-t0) + y0
-end
-
-local current_shape = { {0, 1}, {10*sim.dt, 1} }
-
-
-function discharge_source(tag, x, y, z, t)
- for ii = 1,(#current_shape-1) do
- local t0 = current_shape[ii][1]
- local c0 = current_shape[ii][2]
- local t1 = current_shape[ii+1][1]
- local c1 = current_shape[ii+1][2]
- if (t >= t0 and t < t1) then
- Ez = interp(t, t0, t1, c0, c1)
- return 0, 0, Ez
- end
- end
- return 0, 0, 0
-end
-
-sources[2] = discharge_source
-
-
-
diff --git a/share/examples/capacitor/capacitor.geo b/share/examples/capacitor/capacitor.geo
new file mode 100644
index 0000000000000000000000000000000000000000..3a2306b212ef89647147b81b13ab669689d57781
--- /dev/null
+++ b/share/examples/capacitor/capacitor.geo
@@ -0,0 +1,14 @@
+R = 0.01; // disk radius
+d = 0.001; // dielectric thickness
+t = 0.0003; // disk thickness
+
+SetFactory("OpenCASCADE");
+Mesh.Algorithm3D = 10;
+
+Cylinder(1) = {0, 0, d/2, 0, 0, t, R};
+Cylinder(2) = {0, 0, -d/2, 0, 0, d, R};
+Cylinder(3) = {0, 0, -d/2-t, 0, 0, t, R};
+Sphere(4) = {0, 0, 0, 2*R, -Pi/2, Pi/2, 2*Pi};
+Coherence;
+
+MeshSize{ PointsOf{ Volume{2}; } } = 0.0003;
diff --git a/share/examples/capacitor/params_capacitor.lua b/share/examples/capacitor/params_capacitor.lua
new file mode 100644
index 0000000000000000000000000000000000000000..b694a17e95827437c7cbf831154c707e779b2f45
--- /dev/null
+++ b/share/examples/capacitor/params_capacitor.lua
@@ -0,0 +1,98 @@
+sim.name = "capacitor" -- simulation name
+sim.dt = 1e-14 -- timestep size
+sim.timesteps = 301 -- num of iterations
+sim.gmsh_model = "capacitor.geo" -- gmsh model filename
+sim.use_gpu = 0 -- 0: cpu, 1: gpu
+sim.approx_order = 1 -- approximation order
+sim.time_integrator = "leapfrog" -- time integration method
+postpro.silo_output_rate = 10 -- rate at which to write silo files
+postpro.cycle_print_rate = 10 -- console print rate
+
+local diel_epsr = 10 -- Permittivity of the capacitor dielectric
+
+-- Aluminum plates material parameters
+local alu = {1, 3}
+for i,v in ipairs(alu) do
+ materials[v] = {}
+ materials[v].epsilon = 1
+ materials[v].mu = 1
+ materials[v].sigma = 3.69e7
+end
+
+-- Dielectric parameters
+local diel = { 2 }
+for i,v in ipairs(diel) do
+ materials[v] = {}
+ materials[v].epsilon = diel_epsr
+ materials[v].mu = 1
+ materials[v].sigma = 0
+end
+
+-- Air parameters
+local air = { 4 }
+for i,v in ipairs(air) do
+ materials[v] = {}
+ materials[v].epsilon = 1
+ materials[v].mu = 1
+ materials[v].sigma = 0
+end
+
+-- ** Boundary conditions **
+local absorbing_bcs = {1}
+for i,v in ipairs(absorbing_bcs) do
+ bndconds[v] = {}
+ bndconds[v].kind = "impedance"
+end
+
+function interp(t, t0, t1, y0, y1)
+ return (t-t0)*(y1-y0)/(t1-t0) + y0
+end
+
+local R = 0.01 -- Capacitor disk radius (must match capacitor.geo)
+local d = 0.001 -- Dielectric thickness (must match capacitor.geo)
+local tV = 1 -- Target capacitor voltage
+local cts = 100 -- Charging current pulse timesteps
+
+local ct = cts*sim.dt
+local eps = const.eps0 * diel_epsr
+local A = math.pi*R*R
+local C = eps*A/d
+local J = tV*eps/(ct*d)
+local I = J*A
+local Q = I*ct
+
+local E = Q/(eps*A) -- Expected capacitor field (should equal tV/d)
+
+local do_IO = (not parallel) or (parallel and parallel.comm_rank == 0)
+
+if ( do_IO ) then
+ print("\x1b[1mExpected capacitor field: " .. E .. " V/m")
+ print("Charge: " .. Q .. " C")
+ print("Current density: " .. J .. " A/m^2\x1b[0m")
+end
+
+-- Define the shape of charging current density in terms of
+-- {time, current density} pairs. Linear interpolation between
+-- pairs is used
+local current_shape = { {0, J}, {ct, J} }
+
+-- Define the function evaluating the current source
+function current_source(tag, x, y, z, t)
+ for ii = 1,(#current_shape-1) do
+ local t0 = current_shape[ii][1]
+ local c0 = current_shape[ii][2]
+ local t1 = current_shape[ii+1][1]
+ local c1 = current_shape[ii+1][2]
+ if (t >= t0 and t < t1) then
+ Ez = interp(t, t0, t1, c0, c1)
+ return 0, 0, Ez
+ end
+ end
+ return 0, 0, 0
+end
+
+-- Apply the current source to entity 2
+sources[2] = current_source
+
+
+
diff --git a/share/microwave_oven/oven.geo b/share/examples/microwave_oven/oven.geo
similarity index 100%
rename from share/microwave_oven/oven.geo
rename to share/examples/microwave_oven/oven.geo
diff --git a/share/microwave_oven/params_oven.lua b/share/examples/microwave_oven/params_oven.lua
similarity index 99%
rename from share/microwave_oven/params_oven.lua
rename to share/examples/microwave_oven/params_oven.lua
index 1f0459ae7b66c85593c304167ebc73175738f7e4..99d88da514f1ded3dbad3136ab13d7952f3ca449 100644
--- a/share/microwave_oven/params_oven.lua
+++ b/share/examples/microwave_oven/params_oven.lua
@@ -33,4 +33,3 @@ bndconds[100].kind = "plane_wave_E"
bndconds[100].source = source
-
diff --git a/share/modeconv/modeconv.geo b/share/examples/modeconv/modeconv.geo
similarity index 99%
rename from share/modeconv/modeconv.geo
rename to share/examples/modeconv/modeconv.geo
index b1f7b5d11266daa3eee4cb4dfe59cca69f9bcd19..2c5c74181e0807c5529b0c64080066f721a595be 100644
--- a/share/modeconv/modeconv.geo
+++ b/share/examples/modeconv/modeconv.geo
@@ -13,7 +13,7 @@ z_base = 0;
wg_xlen = 0.14;
wg_ylen = 0.0229;
-wg_zlen = 0.005;
+wg_zlen = 0.01;
rods_start_x = 0.12;
diff --git a/share/modeconv/modeconv.i1 b/share/examples/modeconv/modeconv.i1
similarity index 100%
rename from share/modeconv/modeconv.i1
rename to share/examples/modeconv/modeconv.i1
diff --git a/share/modeconv/params_modeconv.lua b/share/examples/modeconv/params_modeconv.lua
similarity index 95%
rename from share/modeconv/params_modeconv.lua
rename to share/examples/modeconv/params_modeconv.lua
index 20d3c8e99f1610a9bcb8cb7add0a1013108f56f2..bdf106cdf852de7ddd5ab90eba90f1e580e6ced7 100644
--- a/share/modeconv/params_modeconv.lua
+++ b/share/examples/modeconv/params_modeconv.lua
@@ -7,7 +7,7 @@ sim.timesteps = 50000 -- num of iterations
sim.gmsh_model = "modeconv.geo" -- gmsh model filename
sim.use_gpu = 0 -- 0: cpu, 1: gpu
sim.approx_order = 2 -- approximation order
---sim_geom_order = 1 -- geometric order, only 1 for now
+sim.time_integrator = "rk4"
postpro.silo_output_rate = 100 -- rate at which to write silo files
postpro.cycle_print_rate = 10 -- console print rate
diff --git a/share/yagi/params_yagi.lua b/share/examples/yagi/params_yagi.lua
similarity index 87%
rename from share/yagi/params_yagi.lua
rename to share/examples/yagi/params_yagi.lua
index e3d17fdad4abcd9c9a422b5ef7ce4a50415e67d6..5e545ec6b0d0dafdc54e5699593c265dbc8c1648 100644
--- a/share/yagi/params_yagi.lua
+++ b/share/examples/yagi/params_yagi.lua
@@ -3,11 +3,11 @@
]]--
sim.name = "yagi" -- simulation name
sim.dt = 1e-12 -- timestep size
-sim.timesteps = 2000 -- num of iterations
+sim.timesteps = 20000 -- num of iterations
sim.gmsh_model = "yagi.geo" -- gmsh model filename
-sim.use_gpu = 1 -- 0: cpu, 1: gpu
+sim.use_gpu = 0 -- 0: cpu, 1: gpu
sim.approx_order = 1 -- approximation order
---sim_geom_order = 1 -- geometric order, only 1 for now
+sim.time_integrator = "leapfrog"
postpro.silo_output_rate = 100 -- rate at which to write silo files
postpro.cycle_print_rate = 10 -- console print rate
diff --git a/share/yagi/yagi.geo b/share/examples/yagi/yagi.geo
similarity index 100%
rename from share/yagi/yagi.geo
rename to share/examples/yagi/yagi.geo
diff --git a/share/yagi_rad/params_yagi_rad.lua b/share/examples/yagi_rad/params_yagi_rad.lua
similarity index 74%
rename from share/yagi_rad/params_yagi_rad.lua
rename to share/examples/yagi_rad/params_yagi_rad.lua
index de966b78518e76ccfac3d90223b98860353b2359..2090af836ae781570502060765f3167ca696f10d 100644
--- a/share/yagi_rad/params_yagi_rad.lua
+++ b/share/examples/yagi_rad/params_yagi_rad.lua
@@ -3,15 +3,14 @@
The source is modeled as a volumetric current between the
two halves of the radiator.
--]]
-sim.name = "yagi_rad" -- simulation name
+sim.name = "yagi_rad" -- simulation name
sim.dt = 1e-13 -- timestep size
-sim.timesteps = 50000 -- num of iterations
-sim.gmsh_model = "yagi_rad.geo" -- gmsh model filename
-sim.use_gpu = 1 -- 0: cpu, 1: gpu
+sim.timesteps = 50000 -- num of iterations
+sim.gmsh_model = "yagi_rad.geo" -- gmsh model filename
+sim.use_gpu = 0 -- 0: cpu, 1: gpu
sim.approx_order = 2 -- approximation order
---sim_geom_order = 1 -- geometric order, only 1 for now
sim.time_integrator = "leapfrog"
-postpro.silo_output_rate = 500 -- rate at which to write silo files
+postpro.silo_output_rate = 500 -- rate at which to write silo files
postpro.cycle_print_rate = 10 -- console print rate
postpro["H"].silo_mode = "none"
@@ -24,7 +23,7 @@ for i,v in ipairs(alu) do
materials[v] = {}
materials[v].epsilon = 1
materials[v].mu = 1
- materials[v].sigma = 3.69e7
+ materials[v].sigma = 3.69e7
end
local air = {4, 7}
diff --git a/share/yagi_rad/yagi_rad.geo b/share/examples/yagi_rad/yagi_rad.geo
similarity index 100%
rename from share/yagi_rad/yagi_rad.geo
rename to share/examples/yagi_rad/yagi_rad.geo
diff --git a/src/libgmshdg/gmsh_io.cpp b/src/libgmshdg/gmsh_io.cpp
index 5c805974d627e0af69d989812198d4038d6212bc..23f5fb7225fb0b35689990da4999362e00b1eb97 100644
--- a/src/libgmshdg/gmsh_io.cpp
+++ b/src/libgmshdg/gmsh_io.cpp
@@ -410,6 +410,15 @@ model::map_boundaries(void)
#endif /* USE_MPI */
std::vector<bfk_t> bfk = get_facekey_tag_pairs();
+ for (size_t i = 1; i < bfk.size(); i++)
+ {
+ if ( bfk[i-1].first == bfk[i].first )
+ {
+ std::cout << bfk[i-1].second << " " << bfk[i].second << std::endl;
+ throw 42;
+ }
+ }
+
bnd_descriptors.resize( num_faces() );
size_t fbase = 0;
/* For each entity */
diff --git a/src/maxwell/CMakeLists.txt b/src/maxwell/CMakeLists.txt
index e4f2fba6de7d1bd5f9de3f31bca84809453ddd36..00e2876a280c917563397dcb91ad086d415e16f2 100644
--- a/src/maxwell/CMakeLists.txt
+++ b/src/maxwell/CMakeLists.txt
@@ -2,6 +2,7 @@ set(MAXWELL_SOURCES maxwell_interface.cpp
maxwell_cpu.cpp
maxwell_common.cpp
maxwell_solver.cpp
+ maxwell_postpro.cpp
)
if (OPT_ENABLE_GPU_SOLVER)
diff --git a/src/maxwell/maxwell_interface.cpp b/src/maxwell/maxwell_interface.cpp
index b04e16ea468a643a19585f433b6cb8bc7d57cf53..ca5526a836e5734e6edf37b0b733c1cc54e1eae1 100644
--- a/src/maxwell/maxwell_interface.cpp
+++ b/src/maxwell/maxwell_interface.cpp
@@ -86,6 +86,179 @@ void receive_field(field& f, int src, int tag, size_t offset, size_t length, MPI
}
#endif /* USE_MPI */
+field_values::field_values()
+ : Ex(0.0), Ey(0.0), Ez(0.0),
+ Hx(0.0), Hy(0.0), Hz(0.0)
+{}
+
+field_values::field_values(const vec3d& E, const vec3d& H)
+ : Ex( E(0) ), Ey( E(1) ), Ez( E(2) ),
+ Hx( H(0) ), Hy( H(1) ), Hz( H(2) )
+{}
+
+field_values::field_values(double pEx, double pEy, double pEz,
+ double pHx, double pHy, double pHz)
+ : Ex(pEx), Ey(pEy), Ez(pEz), Hx(pHx), Hy(pHy), Hz(pHz)
+{}
+
+field_values&
+field_values::operator+=(const field_values& other)
+{
+ Ex += other.Ex;
+ Ey += other.Ey;
+ Ez += other.Ez;
+ Hx += other.Hx;
+ Hy += other.Hy;
+ Hz += other.Hz;
+ return *this;
+}
+
+field_values
+field_values::operator+(const field_values& other) const
+{
+ field_values ret = *this;
+ ret += other;
+ return ret;
+}
+
+field_values&
+field_values::operator-=(const field_values& other)
+{
+ Ex -= other.Ex;
+ Ey -= other.Ey;
+ Ez -= other.Ez;
+ Hx -= other.Hx;
+ Hy -= other.Hy;
+ Hz -= other.Hz;
+ return *this;
+}
+
+field_values
+field_values::operator-(const field_values& other) const
+{
+ field_values ret = *this;
+ ret -= other;
+ return ret;
+}
+
+field_values&
+field_values::operator*=(const field_values& other)
+{
+ Ex *= other.Ex;
+ Ey *= other.Ey;
+ Ez *= other.Ez;
+ Hx *= other.Hx;
+ Hy *= other.Hy;
+ Hz *= other.Hz;
+ return *this;
+}
+
+field_values
+field_values::operator*(const field_values& other) const
+{
+ field_values ret = *this;
+ ret *= other;
+ return ret;
+}
+
+field_values&
+field_values::operator/=(const field_values& other)
+{
+ Ex /= other.Ex;
+ Ey /= other.Ey;
+ Ez /= other.Ez;
+ Hx /= other.Hx;
+ Hy /= other.Hy;
+ Hz /= other.Hz;
+ return *this;
+}
+
+field_values
+field_values::operator/(const field_values& other) const
+{
+ field_values ret = *this;
+ ret /= other;
+ return ret;
+}
+
+field_values&
+field_values::operator*=(double other)
+{
+ Ex *= other;
+ Ey *= other;
+ Ez *= other;
+ Hx *= other;
+ Hy *= other;
+ Hz *= other;
+ return *this;
+}
+
+field_values
+field_values::operator*(double other) const
+{
+ field_values ret = *this;
+ ret *= other;
+ return ret;
+}
+
+field_values&
+field_values::operator/=(double other)
+{
+ Ex /= other;
+ Ey /= other;
+ Ez /= other;
+ Hx /= other;
+ Hy /= other;
+ Hz /= other;
+ return *this;
+}
+
+field_values
+field_values::operator/(double other) const
+{
+ field_values ret = *this;
+ ret /= other;
+ return ret;
+}
+
+field_values
+operator*(double d, const field_values& m)
+{
+ return m*d;
+}
+
+std::ostream&
+operator<<(std::ostream& os, const field_values& m)
+{
+ os << m.Ex << " " << m.Ey << " " << m.Ez << " ";
+ os << m.Hx << " " << m.Hy << " " << m.Hz;
+ return os;
+}
+
+field_values
+sqrt(const field_values& m)
+{
+ field_values ret;
+ ret.Ex = std::sqrt(m.Ex);
+ ret.Ey = std::sqrt(m.Ey);
+ ret.Ez = std::sqrt(m.Ez);
+ ret.Hx = std::sqrt(m.Hx);
+ ret.Hy = std::sqrt(m.Hy);
+ ret.Hz = std::sqrt(m.Hz);
+ return ret;
+}
+
+double
+hsum(const field_values& m)
+{
+ return m.Ex + m.Ey + m.Ez + m.Hx + m.Hy + m.Hz;
+}
+
+
+
+
+
+
#ifdef ENABLE_GPU_SOLVER
pinned_field::pinned_field()
@@ -239,7 +412,6 @@ field_gpu::copyout(field& emf) const
}
-
material_params_gpu::raw_ptrs
material_params_gpu::data(void)
{
diff --git a/src/maxwell/maxwell_postpro.cpp b/src/maxwell/maxwell_postpro.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9d9d603bd61e1ac6277c3736389af2d9c9b1d97c
--- /dev/null
+++ b/src/maxwell/maxwell_postpro.cpp
@@ -0,0 +1,116 @@
+
+#ifdef USE_MPI
+#include <mpi.h>
+#endif /* USE_MPI */
+
+#include "maxwell/maxwell_postpro.h"
+
+namespace maxwell {
+
+field_values
+eval_field(const field& f, size_t ofs, size_t cbs, const vecxd& phi)
+{
+ assert(ofs+cbs <= f.num_dofs);
+ field_values ret;
+ ret.Ex = f.Ex.segment(ofs, cbs).dot(phi);
+ ret.Ey = f.Ey.segment(ofs, cbs).dot(phi);
+ ret.Ez = f.Ez.segment(ofs, cbs).dot(phi);
+ ret.Hx = f.Hx.segment(ofs, cbs).dot(phi);
+ ret.Hy = f.Hy.segment(ofs, cbs).dot(phi);
+ ret.Hz = f.Hz.segment(ofs, cbs).dot(phi);
+ return ret;
+}
+
+field_values
+compute_error(const model& mod, const solver_state& state, const parameter_loader& mpl)
+{
+ field_values err;
+
+ if (not mpl.has_analytical_solution())
+ throw std::invalid_argument("analytical solution not defined in the Lua configuration");
+
+ for (auto& e : mod)
+ {
+ for (size_t iT = 0; iT < e.num_cells(); iT++)
+ {
+ const auto& pe = e.cell(iT);
+ const auto& re = e.cell_refelem(pe);
+
+ const auto pqps = pe.integration_points();
+ const auto dets = pe.determinants();
+
+ const auto ofs = e.cell_model_dof_offset(iT);
+ const auto cbs = re.num_basis_functions();
+
+ assert(pqps.size() == dets.size());
+ for (size_t iQp = 0; iQp < pqps.size(); iQp++)
+ {
+ const auto& pqp = pqps[iQp];
+ const vecxd phi = re.basis_functions(iQp);
+
+ field_values Fh = eval_field(state.emf_curr, ofs, cbs, phi);
+
+ auto [E, H] = mpl.eval_analytical_solution(e.material_tag(),
+ pqp.point(), state.curr_time);
+ field_values F(E,H);
+ field_values Fdiff = F - Fh;
+ err += Fdiff*Fdiff*pqp.weight();
+ }
+ }
+ }
+
+#ifdef USE_MPI
+ MPI_Allreduce(MPI_IN_PLACE, &err, 6, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+#endif /* USE_MPI */
+ return sqrt(err);
+}
+
+field_values
+compute_energy(const model& mod, const solver_state& state, const parameter_loader& mpl)
+{
+ field_values energy;
+
+ for (auto& e : mod)
+ {
+ for (size_t iT = 0; iT < e.num_cells(); iT++)
+ {
+ const auto& pe = e.cell(iT);
+ const auto& re = e.cell_refelem(pe);
+
+ const auto pqps = pe.integration_points();
+ const auto dets = pe.determinants();
+
+ const auto ofs = e.cell_model_dof_offset(iT);
+ const auto cbs = re.num_basis_functions();
+
+ assert(pqps.size() == dets.size());
+ const auto bar = pe.barycenter();
+ auto eps = mpl.epsilon( e.material_tag(), bar );
+ auto mu = mpl.mu( e.material_tag(), bar );
+
+ field_values mat(eps, eps, eps, mu, mu, mu);
+
+ for (size_t iQp = 0; iQp < pqps.size(); iQp++)
+ {
+ const auto& pqp = pqps[iQp];
+ const vecxd phi = re.basis_functions(iQp);
+ field_values Fh = eval_field(state.emf_curr, ofs, cbs, phi);
+ energy += 0.5*mat*Fh*Fh*pqp.weight();
+ }
+ }
+ }
+
+#ifdef USE_MPI
+ MPI_Allreduce(MPI_IN_PLACE, &energy, 6, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+#endif /* USE_MPI */
+ return energy;
+}
+
+#ifdef ENABLE_GPU_SOLVER
+static void
+validate(const model&, maxwell::solver_state_gpu&,
+ maxwell::parameter_loader&, std::ofstream&)
+{}
+#endif /* ENABLE_GPU_SOLVER */
+
+} // namespace maxwell
\ No newline at end of file
diff --git a/src/maxwell/maxwell_solver.cpp b/src/maxwell/maxwell_solver.cpp
index b6575125392d33bea1544401355639d4a07768da..179afbcb2d57edc34abeb1336b0eb06e23b92d42 100644
--- a/src/maxwell/maxwell_solver.cpp
+++ b/src/maxwell/maxwell_solver.cpp
@@ -19,41 +19,12 @@
#include "libgmshdg/param_loader.h"
#include "maxwell/maxwell_interface.h"
#include "maxwell/maxwell_common.h"
+#include "maxwell/maxwell_postpro.h"
#include "timecounter.h"
#include "sgr.hpp"
using namespace sgr;
-#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)
{
@@ -196,265 +167,10 @@ reinit_materials(const model& mod, maxwell::solver_state& state, maxwell::parame
}
#endif /* ENABLE_EXP_GEOMETRIC_DIODE */
-struct maxwell_6ple {
- double Ex;
- double Ey;
- double Ez;
- double Hx;
- double Hy;
- double Hz;
-
- maxwell_6ple()
- : Ex(0.0), Ey(0.0), Ez(0.0),
- Hx(0.0), Hy(0.0), Hz(0.0)
- {}
-
- maxwell_6ple(const vec3d& E, const vec3d& H)
- : Ex( E(0) ), Ey( E(1) ), Ez( E(2) ),
- Hx( H(0) ), Hy( H(1) ), Hz( H(2) )
- {}
-
- maxwell_6ple(double pEx, double pEy, double pEz, double pHx, double pHy, double pHz)
- : Ex(pEx), Ey(pEy), Ez(pEz), Hx(pHx), Hy(pHy), Hz(pHz)
- {}
-
- maxwell_6ple& operator+=(const maxwell_6ple& other)
- {
- Ex += other.Ex;
- Ey += other.Ey;
- Ez += other.Ez;
- Hx += other.Hx;
- Hy += other.Hy;
- Hz += other.Hz;
- return *this;
- }
-
- maxwell_6ple operator+(const maxwell_6ple& other) const
- {
- maxwell_6ple ret = *this;
- ret += other;
- return ret;
- }
-
- maxwell_6ple& operator-=(const maxwell_6ple& other)
- {
- Ex -= other.Ex;
- Ey -= other.Ey;
- Ez -= other.Ez;
- Hx -= other.Hx;
- Hy -= other.Hy;
- Hz -= other.Hz;
- return *this;
- }
-
- maxwell_6ple operator-(const maxwell_6ple& other) const
- {
- maxwell_6ple ret = *this;
- ret -= other;
- return ret;
- }
-
- maxwell_6ple& operator*=(const maxwell_6ple& other)
- {
- Ex *= other.Ex;
- Ey *= other.Ey;
- Ez *= other.Ez;
- Hx *= other.Hx;
- Hy *= other.Hy;
- Hz *= other.Hz;
- return *this;
- }
-
- maxwell_6ple operator*(const maxwell_6ple& other) const
- {
- maxwell_6ple ret = *this;
- ret *= other;
- return ret;
- }
-
- maxwell_6ple& operator*=(double other)
- {
- Ex *= other;
- Ey *= other;
- Ez *= other;
- Hx *= other;
- Hy *= other;
- Hz *= other;
- return *this;
- }
-
- maxwell_6ple operator*(double other) const
- {
- maxwell_6ple ret = *this;
- ret *= other;
- return ret;
- }
-
- maxwell_6ple& operator/=(const maxwell_6ple& other)
- {
- Ex /= other.Ex;
- Ey /= other.Ey;
- Ez /= other.Ez;
- Hx /= other.Hx;
- Hy /= other.Hy;
- Hz /= other.Hz;
- return *this;
- }
-
- maxwell_6ple operator/(const maxwell_6ple& other) const
- {
- maxwell_6ple ret = *this;
- ret /= other;
- return ret;
- }
-};
-
-static
-maxwell_6ple operator*(double d, const maxwell_6ple& m)
-{
- return m*d;
-}
-
-static std::ostream&
-operator<<(std::ostream& os, const maxwell_6ple& m)
-{
- os << m.Ex << " " << m.Ey << " " << m.Ez << " ";
- os << m.Hx << " " << m.Hy << " " << m.Hz;
- return os;
-}
-
-static maxwell_6ple
-sqrt(const maxwell_6ple& m)
-{
- maxwell_6ple ret;
- ret.Ex = std::sqrt(m.Ex);
- ret.Ey = std::sqrt(m.Ey);
- ret.Ez = std::sqrt(m.Ez);
- ret.Hx = std::sqrt(m.Hx);
- ret.Hy = std::sqrt(m.Hy);
- ret.Hz = std::sqrt(m.Hz);
- return ret;
-}
-
-static double
-hsum(const maxwell_6ple& m)
-{
- return m.Ex + m.Ey + m.Ez + m.Hx + m.Hy + m.Hz;
-}
-
-static maxwell_6ple
-eval_field(const maxwell::field& f, size_t ofs, size_t cbs,
- const vecxd& phi)
-{
- assert(ofs+cbs <= f.num_dofs);
- maxwell_6ple ret;
- ret.Ex = f.Ex.segment(ofs, cbs).dot(phi);
- ret.Ey = f.Ey.segment(ofs, cbs).dot(phi);
- ret.Ez = f.Ez.segment(ofs, cbs).dot(phi);
- ret.Hx = f.Hx.segment(ofs, cbs).dot(phi);
- ret.Hy = f.Hy.segment(ofs, cbs).dot(phi);
- ret.Hz = f.Hz.segment(ofs, cbs).dot(phi);
- return ret;
-}
-
-static maxwell_6ple
-compute_error(const model& mod, maxwell::solver_state& state,
- maxwell::parameter_loader& mpl)
-{
- maxwell_6ple err;
-
- if (not mpl.has_analytical_solution())
- throw std::invalid_argument("analytical solution not defined in the Lua configuration");
-
- for (auto& e : mod)
- {
- for (size_t iT = 0; iT < e.num_cells(); iT++)
- {
- const auto& pe = e.cell(iT);
- const auto& re = e.cell_refelem(pe);
-
- const auto pqps = pe.integration_points();
- const auto dets = pe.determinants();
-
- const auto ofs = e.cell_model_dof_offset(iT);
- const auto cbs = re.num_basis_functions();
-
- assert(pqps.size() == dets.size());
- const auto bar = pe.barycenter();
- auto eps = mpl.epsilon( e.material_tag(), bar );
- auto mu = mpl.mu( e.material_tag(), bar );
-
- for (size_t iQp = 0; iQp < pqps.size(); iQp++)
- {
- const auto& pqp = pqps[iQp];
- const vecxd phi = re.basis_functions(iQp);
-
- maxwell_6ple Fh = eval_field(state.emf_curr, ofs, cbs, phi);
-
- auto [E, H] = mpl.eval_analytical_solution(e.material_tag(),
- pqp.point(), state.curr_time);
- maxwell_6ple F(E,H);
- maxwell_6ple Fdiff = F - Fh;
- err += Fdiff*Fdiff*pqp.weight();
- }
- }
- }
-
-#ifdef USE_MPI
- MPI_Allreduce(MPI_IN_PLACE, &err, 6, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
-#endif /* USE_MPI */
- return sqrt(err);
-}
-
-static maxwell_6ple
-compute_energy(const model& mod, const maxwell::solver_state& state,
- const maxwell::parameter_loader& mpl)
-{
- maxwell_6ple energy;
-
- for (auto& e : mod)
- {
- for (size_t iT = 0; iT < e.num_cells(); iT++)
- {
- const auto& pe = e.cell(iT);
- const auto& re = e.cell_refelem(pe);
-
- const auto pqps = pe.integration_points();
- const auto dets = pe.determinants();
-
- const auto ofs = e.cell_model_dof_offset(iT);
- const auto cbs = re.num_basis_functions();
-
- assert(pqps.size() == dets.size());
- const auto bar = pe.barycenter();
- auto eps = mpl.epsilon( e.material_tag(), bar );
- auto mu = mpl.mu( e.material_tag(), bar );
-
- maxwell_6ple mat(eps, eps, eps, mu, mu, mu);
-
- for (size_t iQp = 0; iQp < pqps.size(); iQp++)
- {
- const auto& pqp = pqps[iQp];
- const vecxd phi = re.basis_functions(iQp);
- maxwell_6ple Fh = eval_field(state.emf_curr, ofs, cbs, phi);
- energy += 0.5*mat*Fh*Fh*pqp.weight();
- }
- }
- }
-
-#ifdef USE_MPI
- MPI_Allreduce(MPI_IN_PLACE, &energy, 6, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
-#endif /* USE_MPI */
- return energy;
-}
-
-#ifdef ENABLE_GPU_SOLVER
-static void
-validate(const model&, maxwell::solver_state_gpu&,
- maxwell::parameter_loader&, std::ofstream&)
-{}
-#endif /* ENABLE_GPU_SOLVER */
+/****************************************************************************/
+/* Solver driver
+ */
template<typename State>
void solver_mainloop(const model& mod, State& state, maxwell::parameter_loader& mpl)
{
@@ -533,19 +249,25 @@ void solver_mainloop(const model& mod, State& state, maxwell::parameter_loader&
std::cout << showcursor << std::endl;
}
+/****************************************************************************/
+/* Register Lua usertypes valid for both CPU and GPU
+ */
static void
register_lua_usertypes(maxwell::parameter_loader& mpl)
{
+ using namespace maxwell;
+
sol::state& lua = mpl.lua_state();
- sol::usertype<maxwell_6ple> maxwell_6ple_ut = lua.new_usertype<maxwell_6ple>("maxwell_6ple",
- sol::constructors<maxwell_6ple()>()
- );
- maxwell_6ple_ut["Ex"] = &maxwell_6ple::Ex;
- maxwell_6ple_ut["Ey"] = &maxwell_6ple::Ey;
- maxwell_6ple_ut["Ez"] = &maxwell_6ple::Ez;
- maxwell_6ple_ut["Hx"] = &maxwell_6ple::Hx;
- maxwell_6ple_ut["Hy"] = &maxwell_6ple::Hy;
- maxwell_6ple_ut["Hz"] = &maxwell_6ple::Hz;
+ sol::usertype<field_values> field_values_ut =
+ lua.new_usertype<field_values>("field_values",
+ sol::constructors<field_values()>()
+ );
+ field_values_ut["Ex"] = &field_values::Ex;
+ field_values_ut["Ey"] = &field_values::Ey;
+ field_values_ut["Ez"] = &field_values::Ez;
+ field_values_ut["Hx"] = &field_values::Hx;
+ field_values_ut["Hy"] = &field_values::Hy;
+ field_values_ut["Hz"] = &field_values::Hz;
#ifdef ENABLE_EXP_GEOMETRIC_DIODE
sol::usertype<point_3d> point_3d_ut = lua.new_usertype<point_3d>("point",
@@ -559,6 +281,9 @@ register_lua_usertypes(maxwell::parameter_loader& mpl)
#endif /* ENABLE_EXP_GEOMETRIC_DIODE */
}
+/****************************************************************************/
+/* Register Lua usertypes valid only for CPU
+ */
static void
register_lua_usertypes_bystate(maxwell::parameter_loader& mpl,
model& mod, maxwell::solver_state& state)
@@ -583,6 +308,9 @@ register_lua_usertypes_bystate(maxwell::parameter_loader& mpl,
#endif /*ENABLE_EXP_GEOMETRIC_DIODE */
}
+/****************************************************************************/
+/* Register Lua usertypes valid only for GPU
+ */
#ifdef ENABLE_GPU_SOLVER
static void
register_lua_usertypes_bystate(maxwell::parameter_loader& mpl,
@@ -590,6 +318,9 @@ register_lua_usertypes_bystate(maxwell::parameter_loader& mpl,
{}
#endif
+/****************************************************************************/
+/* Main
+ */
int main(int argc, char *argv[])
{
#ifdef USE_MPI