diff --git a/doc/changelog.md b/doc/changelog.md
new file mode 100644
index 0000000000000000000000000000000000000000..3994b2bce078975c91d8ed8188880bef46b0d107
--- /dev/null
+++ b/doc/changelog.md
@@ -0,0 +1,13 @@
+# Changelog
+
+## Version v0.1
+Initial version.
+
+Available features:
+- Boundary sources and interface sources both on CPU and GPU. In GPU computations, sources are computed asynchronously on the CPU and then uploaded, and this can be suboptimal if not problematic (see `lua_api.md` for the details). The solver implements `"plane_wave_E"` and `"surface_current"` as boundary sources. It implements `"surface_current"` as interface sources.
+
+Unavailable features:
+- Geometric order is limited to 1. Approximation order is limited to 6.
+- Volumetric sources are not available yet.
+- In the Maxwell solver, the RK4 time integrator limits the practical values of the conducibility `sigma`. If the absolute value of `1 - sigma*sim.dt/(epsilon*const.eps0)` is greater than 1, the computation fails. This will be solved by leapfrog time integration in a subsequent version.
+- Repetitive sources on GPU.
diff --git a/doc/compile.md b/doc/compile.md
new file mode 100644
index 0000000000000000000000000000000000000000..78d7d8ead5c1ce3c0c827e6e162ef998591b20b8
--- /dev/null
+++ b/doc/compile.md
@@ -0,0 +1,21 @@
+# Compile instructions
+
+- You need a C++20 compiler. This code is tested on GCC 10 and Clang 11. GCC 9 is known to *not* work. The code does not use (yet) fancy C++20 features, but it relies on the correct support of the spaceship operator.
+
+- The code does not have many dependencies. It essentially requires Eigen, Lua and SILO. On Debian they are installed with `apt install libeigen3-dev libsilo-dev libsiloh5-0 liblua5.3-dev`.
+
+- The code uses CMake to handle the builds.
+
+The whole procedure should be
+```
+git clone --recursive https://gitlab.onelab.info/mcicuttin/gmsh_gpu_dg.git
+cd gmsh_gpu_dg
+mkdir build
+cd build
+ccmake ..
+make
+```
+
+If you don't compile tests, you will get an executable called `maxwell_solver`. It takes a single parameter, the path of the Lua configuration script (see in `share` for some examples).
+
+To visualize the results you need VisIt (https://wci.llnl.gov/simulation/computer-codes/visit).
\ No newline at end of file
diff --git a/doc/lua_api.md b/doc/lua_api.md
new file mode 100644
index 0000000000000000000000000000000000000000..41ed455647ee915497a62d0bf68a3d86189a52ae
--- /dev/null
+++ b/doc/lua_api.md
@@ -0,0 +1,116 @@
+# GMSH DG solver Lua API
+
+The solver employs the Lua programming language for its **configuration**. Lua was chosen first of all for simplicity: it is extremely lightweight and carries almost no dependencies, and this allow to keep the solver small and compact. Secondly, Lua was chosen to *deliberately* limit the possibilities of what the user can do in the configuration files. If configurations become full-fledged programs, it means that something is missing in the solver core or that the solver is being used in the wrong way.
+
+This file documents the Lua API available in the solver. API has a *general* part and a *problem-specific* part.
+
+The *general* part has to do with configuration not related to a specific problem (i.e. timestep, geometry file), whereas the *problem-specific* part configures all the parameter that make sense only on a given problem (i.e. materials and sources). This separation is reflected also in how the configuration is handled internally.
+
+## API version
+This document describes the API available on the version v0.1 of the solver.
+
+## General interface
+
+### Simulation variables
+
+- `sim.name` (string): name of the simulation. Used also as the oname of the output directory.
+- `sim.dt` (real): timestep duration.
+- `sim.timesteps` (integer): number of timesteps to do.
+- `sim.gmsh_model` (string): name of the file containing the GMSH model.
+- `sim.use_gpu` (0/1): enable/disable GPU usage.
+- `sim.approx_order` (integer): method approximation order.
+- `sim.geom_order` (integer): geometry order. Only order 1 supported for now.
+
+### Postprocessing general variables
+
+- `postpro.silo_output_rate` (integer):
+- `postpro.cycle_print_rate` (integer):
+
+### Callable functions
+- `enable_boundary_sources()`: takes a boolean parameter that specifies if boundary sources should be enabled or not
+- `enable_interface_sources()`: takes a boolean parameter that specifies if the sources applied on internal interfaces should be enabled or not
+- `enable_volume_sources()`: takes a boolean parameter that specifies if volumetric sources should be enabled or not
+
+### Callbacks
+None.
+
+## Maxwell solver interface
+
+### Materials
+Materials are specified populating the `materials` table exposed by the solver. The material parameters are the relative electric permittivity `epsilon`, the relative magnetic permeability `mu` and the conducibility `sigma`. The empty space parameters are stored in `const.eps0` and `const.mu0`.
+
+It is possible either to provide a function describing the material properties in the whole domain or specify the materials subdomain by subdomain. In the first case, appropriate callbacks have to be installed in the material table, whereas in the second case the appropriate variables must be set.
+
+If, for example, the model has two subdomains with tags `1` and `2`, the electric permeability can be specified as
+
+```
+function epsilon_callback(tag, x, y, z)
+ if tag == 1 then
+ return 1
+ return 2
+end
+
+materials.epsilon = epsilon_callback
+```
+or by setting the subdomain-specific variables as
+```
+materials[1] = {}
+materials[1].epsilon = 1
+materials[2] = {}
+materials[2].epsilon = 2
+```
+Note that in this second case the subdomain-specific tables must be initialized with `{}` before setting any variable.
+
+The two approaches can be mixed and the subdomain-specific variables take priority on the whole-domain function.
+
+The solver expects that the callbacks for materials have 4 parameters:
+- `tag`: the current subdomain tag
+- `x`, `y` and `z`: the current coordinates at which the material needs to be evaluated. Note that currently the materials are assumed to be piecewise constant and that the material is evaluated in the barycenter of the element.
+
+### Boundary conditions
+Boundary conditions are specified for each surface with tag `tag` via the table `bndconds`. If a surface does not have an entry in `bndconds` it is assumed to be PEC.
+
+For example, to specify an impedance boundary condition on the surface with tag `4`, one proceeds as following:
+```
+bndconds[4] = {}
+bndconds[4].kind = "impedance"
+```
+
+There are also nonhomogeneous boundary conditions, as the plane wave condition. A plane wave described in terms of electric field is configured as
+```
+function pw_callback(tag, x, y, z, t)
+ --- Computation of the field components
+ return Ex, Ey, Ez
+end
+
+bndconds[4] = {}
+bndconds[4].kind = "plane_wave_E"
+bndconds[4].source = pw_callback
+```
+The available homogeneous conditions are:
+- `"pec"`: Perfect Electric Conductor.
+- `"pmc"`: Perfect Magnetic Conductor.
+- `"impedance"`: Impedance condition, actual impedance value is computed using the material parameters of the adjacent element.
+
+The available nonhomogeneous conditions are:
+- `"E_field"`: electric field. Not implemented yet.
+- `"H_field"`: magnetic field. Not implemented yet.
+- `"surface_current"`: surface current.
+- `"plane_wave_E"`: plane wave specified by the electric field. Normal incidence is assumed.
+- `"plane_wave_H"`: plane wave specified by the magnetic field. Not implemented yet.
+
+### Interface conditions
+The interface conditions follow the exact same logic as boundary conditions, except that they are applied on internal interfaces. The table in this case is `ifaceconds`.
+
+Currently, only a surface current condition is available:
+- `"surface_current"`: surface current.
+### Callable functions
+None.
+
+### Callbacks
+- `on_timestep(ts)`: if defined, the solver calls this function at every timestep. The solver passes the current timestep number in the parameter `ts`. The function is called **before** the solver state is advanced.
+- `electric_initial_condition(x, y, z)`: if defined, the solver calls this function at the beginning of the simulation in order to initialize the electric field. The parameters `x`, `y` and `z` represent the point at which the field needs to be evaluated. The function has to return a triple specifying the three field components, i.e. `return Ex, Ey, Ez`.
+- `magnetic_initial_condition(x, y, z)`: if defined, the solver calls this function at the beginning of the simulation in order to initialize the magnetic field. The parameters `x`, `y` and `z` represent the point at which the field needs to be evaluated. The function has to return a triple specifying the three field components, i.e. `return Hx, Hy, Hz`.
+
+### Considerations on the evaluation of sources and the computation on GPU
+Currently sources for the timestep `t+1` are evaluated asynchronously on the CPU while the GPU is computing the timestep `t`. When the GPU is done, the updated sources are uploaded from the CPU to the GPU and the cycle restarts. This usually works relatively well on big domains with high approximation order, however on small domains at low approximation order does not give any speedup. For this reason, if the sources do not need to be evaluated along the whole simulation, it is suggested to turn them off at the appropriate timestep by using the `on_timestep()` callback and the appropriate functions described below. A future version of the solver will implement repetitive boundary sources, that will be precomputed on the CPU and uploaded on the GPU only once.
\ No newline at end of file
diff --git a/include/maxwell_common.h b/include/maxwell_common.h
index 4de4333057f219a8a356187d5f743ea9c86c0076..a9d0b16386962b8a914b13434149a3dc940da0e5 100644
--- a/include/maxwell_common.h
+++ b/include/maxwell_common.h
@@ -88,5 +88,72 @@ eval_boundary_sources(const model& mod, const parameter_loader& mpl,
}
}
+/* Field must be either 'field' or 'pinned_field'. Too early to use a concept. */
+template<typename State, typename Field>
+void
+eval_interface_sources(const model& mod, const parameter_loader& mpl,
+ State& state, Field& srcfield)
+{
+ auto& bds = mod.boundary_descriptors();
+ assert( srcfield.num_dofs == mod.num_fluxes() );
+
+ size_t face_num_base = 0;
+ for (auto& e : mod)
+ {
+ for (size_t iF = 0; iF < e.num_faces(); iF++)
+ {
+ vec3d n = state.eds[e.number()].normals.row(iF);
+
+ auto gfnum = face_num_base + iF;
+ auto bd = bds[gfnum];
+ if (bd.type == face_type::NONE or bd.type == face_type::BOUNDARY)
+ continue;
+
+ switch ( mpl.interface_type(bd.gmsh_entity) )
+ {
+ case interface_condition::UNSPECIFIED:
+ break;
+
+ case interface_condition::E_FIELD: {
+ /* This is BS actually... */
+ auto fE = [&](const point_3d& pt) -> vec3d {
+ return mpl.eval_interface_source(bd.gmsh_entity, pt, state.curr_time);
+ };
+
+ matxd prjE = e.project_on_face(iF, fE);
+ auto num_bf = prjE.rows();
+ for (size_t i = 0; i < num_bf; i++)
+ {
+ srcfield.Ex[gfnum*num_bf + i] = prjE(i,0);
+ srcfield.Ey[gfnum*num_bf + i] = prjE(i,1);
+ srcfield.Ez[gfnum*num_bf + i] = prjE(i,2);
+ }
+ } break;
+
+ case interface_condition::SURFACE_CURRENT: {
+ auto fJ = [&](const point_3d& pt) -> vec3d {
+ vec3d J = mpl.eval_interface_source(bd.gmsh_entity, pt, state.curr_time);
+ return n.cross(J);
+ };
+
+ matxd prjJ = e.project_on_face(iF, fJ);
+ auto num_bf = prjJ.rows();
+ for (size_t i = 0; i < num_bf; i++)
+ {
+ srcfield.Hx[gfnum*num_bf + i] = prjJ(i,0);
+ srcfield.Hy[gfnum*num_bf + i] = prjJ(i,1);
+ srcfield.Hz[gfnum*num_bf + i] = prjJ(i,2);
+ }
+ } break;
+
+ case interface_condition::H_FIELD:
+ break;
+ }
+ }
+
+ face_num_base += e.num_faces();
+ }
+}
+
} // namespace maxwell
diff --git a/include/maxwell_interface.h b/include/maxwell_interface.h
index 35defba4a86db21a279fecfe47accb61b30f1c0a..92db3e188b8474ed8afb8b86868e587419b71a36 100644
--- a/include/maxwell_interface.h
+++ b/include/maxwell_interface.h
@@ -310,7 +310,7 @@ void swap(maxwell::solver_state_gpu&);
void init_boundary_conditions(const model&, const parameter_loader&, vecxd&);
void init_lax_milgram(const model&, const parameter_loader&, vecxd&, vecxd&, vecxd&, vecxd&);
//void eval_boundary_sources(const model&, const parameter_loader&, solver_state&, field&);
-void eval_interface_sources(const model&, const parameter_loader&, solver_state&, field&);
+//void eval_interface_sources(const model&, const parameter_loader&, solver_state&, field&);
#endif /* HIDE_THIS_FROM_NVCC */
template<typename Function>
diff --git a/share/binde_simplified/binde_simplified.geo b/share/binde_simplified/binde_simplified.geo
new file mode 100644
index 0000000000000000000000000000000000000000..fd1b270d91bbf427951e1ca30508ad7b66502856
--- /dev/null
+++ b/share/binde_simplified/binde_simplified.geo
@@ -0,0 +1,10 @@
+SetFactory("OpenCASCADE");
+Box(1) = {0, 0, 0, 3, 3, 3};
+Box(2) = {1+0.3, 1+0.45, 1+0.45, 0.4, 0.1, 0.1};
+Coherence;
+
+// Air box
+MeshSize {17, 18, 19, 20, 21, 22, 23, 24} = 0.2;
+
+// Metal block
+MeshSize {9, 10, 11, 12, 13, 14, 15, 16} = 0.05;
diff --git a/share/binde_simplified/params_binde_simplified.lua b/share/binde_simplified/params_binde_simplified.lua
new file mode 100644
index 0000000000000000000000000000000000000000..e47b43b1bf38d7b5921650266900119a6da5541a
--- /dev/null
+++ b/share/binde_simplified/params_binde_simplified.lua
@@ -0,0 +1,74 @@
+sim.name = "binde_simplified" -- simulation name
+sim.dt = 1e-11 -- timestep size
+sim.timesteps = 10000 -- num of iterations
+sim.gmsh_model = "binde_simplified.geo" -- gmsh model filename
+sim.use_gpu = 1 -- 0: cpu, 1: gpu
+sim.approx_order = 1 -- approximation order
+--sim_geom_order = 1 -- geometric order, only 1 for now
+postpro.silo_output_rate = 10 -- rate at which to write silo files
+postpro.cycle_print_rate = 10 -- console print rate
+
+-- ** Materials **
+-- Aluminum
+local alu = {2}
+for i,v in ipairs(alu) do
+ materials[v] = {}
+ materials[v].epsilon = 1
+ materials[v].mu = 1
+ materials[v].sigma = 1
+
+ local coeff = sim.dt*materials[v].sigma/(materials[v].epsilon*const.eps0)
+ if math.abs(1-coeff) > 1 then
+ print("WARNING: sigma too high on domain " .. v .. ", this simulation is going to explode")
+ end
+end
+
+local air = {3}
+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 Z_bnds = {13, 14, 15, 16, 17, 18}
+for i,v in ipairs(Z_bnds) 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, 0}, {9.5e-9, 0.03}, {1.1e-8, 1},
+ {1.4e-8, 0.45}, {1.6e-8, 0.45}, {2.7e-8, 0.6},
+ {9e-8, 0.01}, {1e-7, 0} }
+
+function interface_source_7(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
+
+local discharge_surf = 7
+ifaceconds[discharge_surf] = {}
+ifaceconds[discharge_surf].kind = "surface_current"
+ifaceconds[discharge_surf].source = interface_source_7
+
+function on_timestep(ts)
+ --local time = ts*sim.dt
+ --local curr = interface_source_7(0,0,0,0,time)
+ --print("Time: " .. time .. ", current: " .. curr)
+end
+
+
diff --git a/share/modeconv/modeconv.geo b/share/modeconv/modeconv.geo
new file mode 100644
index 0000000000000000000000000000000000000000..b1f7b5d11266daa3eee4cb4dfe59cca69f9bcd19
--- /dev/null
+++ b/share/modeconv/modeconv.geo
@@ -0,0 +1,78 @@
+// Mode converter GMSH model
+
+SetFactory("OpenCASCADE");
+
+Mesh.Algorithm = 1;
+Mesh.Algorithm3D = 1;
+
+cl = 0.01;
+
+x_base = 0;
+y_base = 0;
+z_base = 0;
+
+wg_xlen = 0.14;
+wg_ylen = 0.0229;
+wg_zlen = 0.005;
+
+rods_start_x = 0.12;
+
+rod_spacing = 0.00954;
+
+rod1_d = 0.01145; rod1_r = 0.000380; rod1_x = x_base + rods_start_x - rod_spacing*0; rod1_y = y_base + rod1_d; rod1_z = z_base;
+rod2_d = 0.01145; rod2_r = 0.000550; rod2_x = x_base + rods_start_x - rod_spacing*1; rod2_y = y_base + rod2_d; rod2_z = z_base;
+rod3_d = 0.01145; rod3_r = 0.000550; rod3_x = x_base + rods_start_x - rod_spacing*2; rod3_y = y_base + rod3_d; rod3_z = z_base;
+rod4_d = 0.01145; rod4_r = 0.000550; rod4_x = x_base + rods_start_x - rod_spacing*3; rod4_y = y_base + rod4_d; rod4_z = z_base;
+rod5_d = 0.01045; rod5_r = 0.000515; rod5_x = x_base + rods_start_x - rod_spacing*4; rod5_y = y_base + rod5_d; rod5_z = z_base;
+rod6_d = 0.00900; rod6_r = 0.000570; rod6_x = x_base + rods_start_x - rod_spacing*5; rod6_y = y_base + rod6_d; rod6_z = z_base;
+rod7_d = 0.00700; rod7_r = 0.000640; rod7_x = x_base + rods_start_x - rod_spacing*6; rod7_y = y_base + rod7_d; rod7_z = z_base;
+rod8_d = 0.00500; rod8_r = 0.000690; rod8_x = x_base + rods_start_x - rod_spacing*7; rod8_y = y_base + rod8_d; rod8_z = z_base;
+rod9_d = 0.00300; rod9_r = 0.000760; rod9_x = x_base + rods_start_x - rod_spacing*8; rod9_y = y_base + rod9_d; rod9_z = z_base;
+rod10_d = 0.00125; rod10_r = 0.000920; rod10_x = x_base + rods_start_x - rod_spacing*9; rod10_y = y_base + rod10_d; rod9_z = z_base;
+rod11_d = 0.02165; rod11_r = 0.000960; rod11_x = x_base + rods_start_x - rod_spacing*9; rod11_y = y_base + rod11_d; rod9_z = z_base;
+
+cl_rod = 0.0005;
+xb = 100; rx = rod1_x; ry = rod1_y; rz = z_base; rr = rod1_r; Include "modeconv.i1";
+xb = 200; rx = rod2_x; ry = rod2_y; rz = z_base; rr = rod2_r; Include "modeconv.i1";
+xb = 300; rx = rod3_x; ry = rod3_y; rz = z_base; rr = rod3_r; Include "modeconv.i1";
+xb = 400; rx = rod4_x; ry = rod4_y; rz = z_base; rr = rod4_r; Include "modeconv.i1";
+xb = 500; rx = rod5_x; ry = rod5_y; rz = z_base; rr = rod5_r; Include "modeconv.i1";
+xb = 600; rx = rod6_x; ry = rod6_y; rz = z_base; rr = rod6_r; Include "modeconv.i1";
+xb = 700; rx = rod7_x; ry = rod7_y; rz = z_base; rr = rod7_r; Include "modeconv.i1";
+xb = 800; rx = rod8_x; ry = rod8_y; rz = z_base; rr = rod8_r; Include "modeconv.i1";
+xb = 900; rx = rod9_x; ry = rod9_y; rz = z_base; rr = rod9_r; Include "modeconv.i1";
+xb = 1000; rx = rod10_x; ry = rod10_y; rz = z_base; rr = rod10_r; Include "modeconv.i1";
+xb = 1100; rx = rod11_x; ry = rod11_y; rz = z_base; rr = rod11_r; Include "modeconv.i1";
+
+Point(1) = { x_base, y_base, z_base, cl/3 };
+Point(2) = { x_base + wg_xlen, y_base, z_base, cl };
+Point(3) = { x_base + wg_xlen, y_base + wg_ylen, z_base, cl };
+Point(4) = { x_base, y_base + wg_ylen, z_base, cl/3 };
+
+Line(1) = {1,2};
+Line(2) = {2,3};
+Line(3) = {3,4};
+Line(4) = {4,1};
+
+Curve Loop(5) = {1,2,3,4};
+
+Plane Surface(6) = {5};
+
+Extrude {0, 0, wg_zlen} {
+ Surface { 6, 121, 221, 321, 421, 521, 621, 721, 821, 921, 1021, 1121 };
+ //Surface {6};
+}
+
+BooleanDifference(100) = { Volume{1}; Delete; } { Volume{2}; };
+BooleanDifference(101) = { Volume{100}; Delete; } { Volume{3}; };
+BooleanDifference(102) = { Volume{101}; Delete; } { Volume{4}; };
+BooleanDifference(103) = { Volume{102}; Delete; } { Volume{5}; };
+BooleanDifference(104) = { Volume{103}; Delete; } { Volume{6}; };
+BooleanDifference(105) = { Volume{104}; Delete; } { Volume{7}; };
+BooleanDifference(106) = { Volume{105}; Delete; } { Volume{8}; };
+BooleanDifference(107) = { Volume{106}; Delete; } { Volume{8}; };
+BooleanDifference(108) = { Volume{107}; Delete; } { Volume{9}; };
+BooleanDifference(109) = { Volume{108}; Delete; } { Volume{10}; };
+BooleanDifference(110) = { Volume{109}; Delete; } { Volume{11}; };
+BooleanDifference(111) = { Volume{110}; Delete; } { Volume{12}; };
+
diff --git a/share/modeconv/modeconv.i1 b/share/modeconv/modeconv.i1
new file mode 100644
index 0000000000000000000000000000000000000000..21101d3298d4a6661261ca33417e49def081afe6
--- /dev/null
+++ b/share/modeconv/modeconv.i1
@@ -0,0 +1,16 @@
+// Mode converter GMSH model
+
+Point(xb + 01) = { rx, ry, rz, cl_rod };
+Point(xb + 02) = { rx + rr, ry, rz, cl_rod };
+Point(xb + 03) = { rx, ry + rr, rz, cl_rod };
+Point(xb + 04) = { rx - rr, ry, rz, cl_rod };
+Point(xb + 05) = { rx, ry - rr, rz, cl_rod };
+
+Circle(xb + 01) = {xb + 02, xb + 01, xb + 03};
+Circle(xb + 02) = {xb + 03, xb + 01, xb + 04};
+Circle(xb + 03) = {xb + 04, xb + 01, xb + 05};
+Circle(xb + 04) = {xb + 05, xb + 01, xb + 02};
+
+Curve Loop(xb + 01) = {xb + 01, xb + 02, xb + 03, xb + 04};
+Plane Surface(xb + 21) = {xb + 01};
+
diff --git a/share/modeconv/params_modeconv.lua b/share/modeconv/params_modeconv.lua
new file mode 100644
index 0000000000000000000000000000000000000000..67b816ce7d2421c1096bbb031090bc47d4d97c57
--- /dev/null
+++ b/share/modeconv/params_modeconv.lua
@@ -0,0 +1,49 @@
+sim.name = "modeconv" -- simulation name
+sim.dt = 1e-14 -- timestep size
+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
+postpro.silo_output_rate = 100 -- rate at which to write silo files
+postpro.cycle_print_rate = 10 -- console print rate
+
+-- ** Materials **
+-- Plastic rods
+local rods = {2,3,4,5,6,7,8,9,10,11,12}
+for i,v in ipairs(rods) do
+ materials[v] = {}
+ materials[v].epsilon = 24
+ materials[v].mu = 1
+ materials[v].sigma = 0
+end
+-- Air
+local air = 111
+materials[air] = {}
+materials[air].epsilon = 1
+materials[air].mu = 1
+materials[air].sigma = 0
+
+-- ** Boundary conditions **
+-- These conditions are not correct for this model, is a test.
+-- waveguide end => impedance
+local Z_bnds = {1184}
+for i,v in ipairs(Z_bnds) do
+ bndconds[v] = {}
+ bndconds[v].kind = "impedance"
+end
+
+-- waveguide start => plane wave source
+local freq = 1.55e10
+function boundary_source_1183(tag, x, y, z, t)
+ local Ex = 0
+ local Ey = 0
+ local Ez = math.sin(2*math.pi*t*freq)
+ return Ex, Ey, Ez
+end
+
+bndconds[1183] = {}
+bndconds[1183].kind = "plane_wave_E"
+bndconds[1183].source = boundary_source_1183
+
+
diff --git a/share/yagi/params_yagi.lua b/share/yagi/params_yagi.lua
new file mode 100644
index 0000000000000000000000000000000000000000..4d73f80e198595d559bdbd73c33462ead7cac9dc
--- /dev/null
+++ b/share/yagi/params_yagi.lua
@@ -0,0 +1,59 @@
+sim.name = "yagi" -- simulation name
+sim.dt = 1e-12 -- timestep size
+sim.timesteps = 2000 -- num of iterations
+sim.gmsh_model = "yagi.geo" -- gmsh model filename
+sim.use_gpu = 1 -- 0: cpu, 1: gpu
+sim.approx_order = 1 -- approximation order
+--sim_geom_order = 1 -- geometric order, only 1 for now
+postpro.silo_output_rate = 100 -- rate at which to write silo files
+postpro.cycle_print_rate = 10 -- console print rate
+
+-- ** Materials **
+-- Aluminum
+local alu = {2, 3, 4}
+for i,v in ipairs(alu) do
+ materials[v] = {}
+ materials[v].epsilon = 1
+ materials[v].mu = 1
+ materials[v].sigma = 17
+
+ local coeff = sim.dt*materials[v].sigma/(materials[v].epsilon*const.eps0)
+ if math.abs(1-coeff) > 1 then
+ print("WARNING: sigma too high on domain " .. v .. ", this simulation is going to explode")
+ end
+end
+
+local air = {5}
+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 Z_bnds = {16}
+for i,v in ipairs(Z_bnds) do
+ bndconds[v] = {}
+ bndconds[v].kind = "impedance"
+end
+
+local pmc_bnds = {17, 19}
+for i,v in ipairs(pmc_bnds) do
+ bndconds[v] = {}
+ bndconds[v].kind = "pmc"
+end
+
+local freq = 3e8
+function interface_source_7(tag, x, y, z, t)
+ local Ex = 0
+ local Ey = 0
+ local Ez = math.sin(2*math.pi*t*freq)
+ return Ex, Ey, Ez
+end
+
+bndconds[21] = {}
+bndconds[21].kind = "plane_wave_E"
+bndconds[21].source = interface_source_7
+
+
diff --git a/share/yagi/yagi.geo b/share/yagi/yagi.geo
new file mode 100644
index 0000000000000000000000000000000000000000..2cdb917cdef2cdec74add8e30d0754abb428cdb5
--- /dev/null
+++ b/share/yagi/yagi.geo
@@ -0,0 +1,18 @@
+SetFactory("OpenCASCADE");
+
+MHz = 300;
+
+R = 0.5*(300/MHz);
+A = 0.44*(300/MHz);
+D = 0.43*(300/MHz);
+sp_RA = 0.25*(300/MHz);
+sp_AD = 0.13*(300/MHz);
+
+Box(1) = {0, 0, 0, 1, 1, 1};
+Cylinder(2) = {0.5-sp_RA, 0.5, 0.5-R/2, 0, 0, R, 0.01, 2*Pi};
+Cylinder(3) = {0.5, 0.5, 0.5-A/2, 0, 0, A, 0.01, 2*Pi};
+Cylinder(4) = {0.5+sp_AD, 0.5, 0.5-D/2, 0, 0, D, 0.01, 2*Pi};
+Coherence;
+
+MeshSize {9,10,11,12,13,14} = 0.01;
+//MeshSize {15,16,17,18,19,20,21,22} = 0.1;
diff --git a/src/maxwell_cpu.cpp b/src/maxwell_cpu.cpp
index 8da503207ae6f71b8e08987848996a92fe575b0e..5099f17417b3494f7a238f10a349a5f473b974a4 100644
--- a/src/maxwell_cpu.cpp
+++ b/src/maxwell_cpu.cpp
@@ -51,64 +51,7 @@ init_matparams(const model& mod, solver_state& state,
-void
-eval_interface_sources(const model& mod, const parameter_loader& mpl,
- solver_state& state, field& srcfield)
-{
- auto& bds = mod.boundary_descriptors();
- assert( srcfield.num_dofs == mod.num_fluxes() );
-
- size_t face_num_base = 0;
- for (auto& e : mod)
- {
- for (size_t iF = 0; iF < e.num_faces(); iF++)
- {
- vec3d n = state.eds[e.number()].normals.row(iF);
-
- auto gfnum = face_num_base + iF;
- auto bd = bds[gfnum];
- if (bd.type == face_type::NONE or bd.type == face_type::BOUNDARY)
- continue;
- switch ( mpl.interface_type(bd.gmsh_entity) )
- {
- case interface_condition::UNSPECIFIED:
- break;
-
- case interface_condition::E_FIELD: {
- /* This is BS actually... */
- auto fE = [&](const point_3d& pt) -> vec3d {
- return mpl.eval_interface_source(bd.gmsh_entity, pt, state.curr_time);
- };
-
- matxd prjE = e.project_on_face(iF, fE);
- auto num_bf = prjE.rows();
- srcfield.Ex.segment(gfnum*num_bf, num_bf) = prjE.col(0);
- srcfield.Ey.segment(gfnum*num_bf, num_bf) = prjE.col(1);
- srcfield.Ez.segment(gfnum*num_bf, num_bf) = prjE.col(2);
- } break;
-
- case interface_condition::SURFACE_CURRENT: {
- auto fJ = [&](const point_3d& pt) -> vec3d {
- vec3d J = mpl.eval_interface_source(bd.gmsh_entity, pt, state.curr_time);
- return n.cross(J);
- };
-
- matxd prjJ = e.project_on_face(iF, fJ);
- auto num_bf = prjJ.rows();
- srcfield.Hx.segment(gfnum*num_bf, num_bf) = prjJ.col(0);
- srcfield.Hy.segment(gfnum*num_bf, num_bf) = prjJ.col(1);
- srcfield.Hz.segment(gfnum*num_bf, num_bf) = prjJ.col(2);
- } break;
-
- case interface_condition::H_FIELD:
- break;
- }
- }
-
- face_num_base += e.num_faces();
- }
-}
void
init_from_model(const model& mod, solver_state& state)
diff --git a/src/maxwell_gpu.cpp b/src/maxwell_gpu.cpp
index 589dbce355cb03b396ea6b00be046a811c485939..a9c79c7365b86d0b93380c84d44169341cf5a1e0 100644
--- a/src/maxwell_gpu.cpp
+++ b/src/maxwell_gpu.cpp
@@ -241,6 +241,9 @@ prepare_sources(const model& mod, maxwell::solver_state_gpu& state,
if ( mpl.boundary_sources_enabled() )
maxwell::eval_boundary_sources(mod, mpl, state, state.bndsrcs_cpu);
+ if ( mpl.interface_sources_enabled() )
+ maxwell::eval_interface_sources(mod, mpl, state, state.bndsrcs_cpu);
+
state.bndsrcs_buf.copyin(state.bndsrcs_cpu, state.memcpy_stream);
state.memcpy_stream.wait();
std::swap(state.bndsrcs, state.bndsrcs_buf);
@@ -263,6 +266,9 @@ do_sources(const model& mod, maxwell::solver_state_gpu& state,
if ( be )
maxwell::eval_boundary_sources(mod, mpl, state, state.bndsrcs_cpu);
+ if ( ie )
+ maxwell::eval_interface_sources(mod, mpl, state, state.bndsrcs_cpu);
+
if ( be or ie or ve )
state.bndsrcs_buf.copyin(state.bndsrcs_cpu, state.memcpy_stream);
}
diff --git a/src/maxwell_solver.cpp b/src/maxwell_solver.cpp
index 2aab09f72faffb75dc80cbf508ba1bfeb52d7e29..a8b323d94027e8f23f23646e3507013fd12c460c 100644
--- a/src/maxwell_solver.cpp
+++ b/src/maxwell_solver.cpp
@@ -126,20 +126,20 @@ int main(int argc, const char *argv[])
return 1;
}
-#ifdef ENABLE_GPU_SOLVER
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
{
-#endif /* ENABLE_GPU_SOLVER */
maxwell::solver_state state_c;
test_it(mod, state_c, mpl);
-#ifdef ENABLE_GPU_SOLVER
}
-#endif
//gmsh::finalize();