Leapfrog integrator. Still experimental, don't use.

doc/changelog.md

@@ -18,3 +18,5 @@ Unavailable features:
 - Switched back to C++17 due to various problems on Eigen and Mac OS X.
 - Use system GMSH if installed and only after look at the hardcoded path (a55af8ab).
 - Improved export to SILO. E, H, and J fields are now exportable as nodal or zonal variables. It is also possible to disable the export of a specific field, see `lua_api.md` (8a831181).
+- Preliminary implementation of leapfrog integration on CPU
+

doc/lua_api.md

@@ -20,6 +20,7 @@ This document describes the API available on the version v0.1 of the solver.
 - `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.
+-  sim.time_integrator (string): `"rk4"`, `"leapfrog"` or `"euler"`. Euler is only for test purporses.
 
 ### Postprocessing general variables
 

include/maxwell_common.h

@@ -57,8 +57,8 @@ eval_boundary_sources(const model& mod, const parameter_loader& mpl,
                         return mpl.eval_boundary_source(bd.gmsh_entity, pt, state.curr_time);
                     };
                     auto fH = [&](const point_3d& pt) -> vec3d {
-                        vec3d E = mpl.eval_boundary_source(bd.gmsh_entity, pt, state.curr_time);
-                        return n.cross(E)/Z;
+                        vec3d H = mpl.eval_boundary_source(bd.gmsh_entity, pt, state.curr_time);
+                        return n.cross(H)/Z;
                     };
                     matxd prjE = e.project_on_face(iF, fE);
                     auto num_bf = prjE.rows();

include/maxwell_interface.h

@@ -292,7 +292,7 @@ void init_from_model(const model&, solver_state&);
 void init_matparams(const model&, solver_state&, const parameter_loader&);
 //void apply_operator(solver_state&, const field&, field&);
 void export_fields_to_silo(const model&, const solver_state&, const parameter_loader&);
-void timestep(solver_state&);
+void timestep(solver_state&, time_integrator_type);
 void prepare_sources(const model&, maxwell::solver_state&, maxwell::parameter_loader&);
 void do_sources(const model&, maxwell::solver_state&, maxwell::parameter_loader&);
 void swap(maxwell::solver_state&);
@@ -302,7 +302,7 @@ void init_from_model(const model&, solver_state_gpu&);
 void init_matparams(const model&, solver_state_gpu&, const parameter_loader&);
 //void apply_operator(solver_state_gpu&, const field_gpu&, field_gpu&);
 void export_fields_to_silo(const model&, const solver_state_gpu&, const parameter_loader&);
-void timestep(solver_state_gpu&);
+void timestep(solver_state_gpu&, time_integrator_type);
 void prepare_sources(const model&, maxwell::solver_state_gpu&, maxwell::parameter_loader&);
 void do_sources(const model&, maxwell::solver_state_gpu&, maxwell::parameter_loader&);
 void swap(maxwell::solver_state_gpu&);

include/param_loader.h

@@ -7,6 +7,18 @@
 
 #include "gmsh_io.h"
 
+enum class field_export_mode {
+    NONE,
+    NODAL,
+    ZONAL
+};
+
+enum class time_integrator_type {
+    RK4,
+    LEAPFROG,
+    EULER
+};
+
 class parameter_loader_base
 {
     bool        m_bnd_sources_active;
@@ -33,6 +45,7 @@ public:
     size_t          sim_timesteps(void) const;
     std::string     sim_gmshmodel(void) const;
     bool            sim_usegpu(void) const;
+    time_integrator_type    sim_timeIntegrator(void) const;
     size_t          postpro_siloOutputRate(void) const;
     size_t          postpro_cyclePrintRate(void) const;
 
@@ -50,9 +63,5 @@ public:
     void            call_timestep_callback(size_t);
 };
 
-enum class field_export_mode {
-    NONE,
-    NODAL,
-    ZONAL
-};
+
 

src/maxwell_cpu.cpp

@@ -105,6 +105,37 @@ compute_curls(solver_state& state, const field& curr, field& next)
     next.Hz = state.dy.Ex - state.dx.Ey;
 }
 
+static void
+compute_curls_E(solver_state& state, const field& curr, 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);
+    }
+
+    next.Hx = state.dz.Ey - state.dy.Ez;
+    next.Hy = state.dx.Ez - state.dz.Ex;
+    next.Hz = state.dy.Ex - state.dx.Ey;
+}
+
+static void
+compute_curls_H(solver_state& state, const field& curr, field& next)
+{
+    for (const auto& ed : state.eds)
+    {
+        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;
+}
+
+
 static void
 compute_field_jumps(solver_state& state, const field& in)
 {
@@ -143,6 +174,70 @@ compute_field_jumps(solver_state& state, const field& in)
     }
 }
 
+static void
+compute_field_jumps_E(solver_state& state, const field& in)
+{
+    for (auto& ed : state.eds)
+    {
+        size_t num_all_fluxes = num_elems_all_orientations(ed) * ed.num_faces_per_elem * ed.num_fluxes;
+ 
+        #pragma omp parallel for
+        for (size_t i = 0; i < num_all_fluxes; i++)
+        {
+            auto lofs = i;
+            auto gofs = ed.flux_base + i;
+            if (ed.fluxdofs_neigh[lofs] != NOT_PRESENT)
+            {
+                auto pos_mine  = ed.fluxdofs_mine[lofs];
+                auto pos_neigh = ed.fluxdofs_neigh[lofs];
+                state.jumps.Ex[gofs] = in.Ex[pos_mine] - in.Ex[pos_neigh];
+                state.jumps.Ey[gofs] = in.Ey[pos_mine] - in.Ey[pos_neigh];
+                state.jumps.Ez[gofs] = in.Ez[pos_mine] - in.Ez[pos_neigh];
+            }
+            else
+            {
+                auto bc = state.matparams.bc_coeffs[gofs];
+                auto pos_mine  = ed.fluxdofs_mine[lofs];
+                state.jumps.Ex[gofs] = bc * in.Ex[pos_mine];
+                state.jumps.Ey[gofs] = bc * in.Ey[pos_mine];
+                state.jumps.Ez[gofs] = bc * in.Ez[pos_mine];
+            }
+        }
+    }
+}
+
+static void
+compute_field_jumps_H(solver_state& state, const field& in)
+{
+    for (auto& ed : state.eds)
+    {
+        size_t num_all_fluxes = num_elems_all_orientations(ed) * ed.num_faces_per_elem * ed.num_fluxes;
+ 
+        #pragma omp parallel for
+        for (size_t i = 0; i < num_all_fluxes; i++)
+        {
+            auto lofs = i;
+            auto gofs = ed.flux_base + i;
+            if (ed.fluxdofs_neigh[lofs] != NOT_PRESENT)
+            {
+                auto pos_mine  = ed.fluxdofs_mine[lofs];
+                auto pos_neigh = ed.fluxdofs_neigh[lofs];
+                state.jumps.Hx[gofs] = in.Hx[pos_mine] - in.Hx[pos_neigh];
+                state.jumps.Hy[gofs] = in.Hy[pos_mine] - in.Hy[pos_neigh];
+                state.jumps.Hz[gofs] = in.Hz[pos_mine] - in.Hz[pos_neigh];
+            }
+            else
+            {
+                auto bc = state.matparams.bc_coeffs[gofs];
+                auto pos_mine  = ed.fluxdofs_mine[lofs];
+                state.jumps.Hx[gofs] = (2.0 - bc) * in.Hx[pos_mine];
+                state.jumps.Hy[gofs] = (2.0 - bc) * in.Hy[pos_mine];
+                state.jumps.Hz[gofs] = (2.0 - bc) * in.Hz[pos_mine];
+            }
+        }
+    }
+}
+
 static void
 compute_fluxes_planar(solver_state& state)
 {
@@ -192,6 +287,82 @@ compute_fluxes_planar(solver_state& state)
     }
 }
 
+/* In leapfrog we use centered fluxes for the moment */
+static void
+compute_fluxes_planar_E(solver_state& state)
+{
+    for (auto& ed : state.eds)
+    {
+        size_t num_all_faces = num_elems_all_orientations(ed) * ed.num_faces_per_elem;
+
+        #pragma omp parallel for
+        for (size_t iF = 0; iF < num_all_faces; iF++)
+        {
+            auto n = ed.normals.row(iF);
+            auto nx = n[0];
+            auto ny = n[1];
+            auto nz = n[2];
+            auto face_det = ed.face_determinants(iF);
+
+            for (size_t i = 0; i < ed.num_fluxes; i++)
+            {
+                auto lofs = iF*ed.num_fluxes + i;
+                auto gofs = ed.flux_base + lofs;
+
+                /* Sources are imposed on jumps */
+                auto jHx = state.jumps.Hx[gofs] + state.bndsrcs.Hx[gofs];
+                auto jHy = state.jumps.Hy[gofs] + state.bndsrcs.Hy[gofs];
+                auto jHz = state.jumps.Hz[gofs] + state.bndsrcs.Hz[gofs];
+
+                auto aE = face_det * 0.5;
+
+                /* Compute fluxes */
+                state.fluxes.Ex[gofs] = aE*(nz*jHy - ny*jHz);
+                state.fluxes.Ey[gofs] = aE*(nx*jHz - nz*jHx);
+                state.fluxes.Ez[gofs] = aE*(ny*jHx - nx*jHy);
+            }
+        }
+    }
+}
+
+/* In leapfrog we use centered fluxes for the moment */
+static void
+compute_fluxes_planar_H(solver_state& state)
+{
+    for (auto& ed : state.eds)
+    {
+        size_t num_all_faces = num_elems_all_orientations(ed) * ed.num_faces_per_elem;
+
+        #pragma omp parallel for
+        for (size_t iF = 0; iF < num_all_faces; iF++)
+        {
+            auto n = ed.normals.row(iF);
+            auto nx = n[0];
+            auto ny = n[1];
+            auto nz = n[2];
+            auto face_det = ed.face_determinants(iF);
+
+            for (size_t i = 0; i < ed.num_fluxes; i++)
+            {
+                auto lofs = iF*ed.num_fluxes + i;
+                auto gofs = ed.flux_base + lofs;
+
+                /* Sources are imposed on jumps */
+                auto jEx = state.jumps.Ex[gofs] - state.bndsrcs.Ex[gofs];
+                auto jEy = state.jumps.Ey[gofs] - state.bndsrcs.Ey[gofs];
+                auto jEz = state.jumps.Ez[gofs] - state.bndsrcs.Ez[gofs];
+
+                auto aH = face_det * 0.5;
+
+                /* Compute fluxes */
+                state.fluxes.Hx[gofs] = aH*(ny*jEz - nz*jEy);
+                state.fluxes.Hy[gofs] = aH*(nz*jEx - nx*jEz);
+                state.fluxes.Hz[gofs] = aH*(nx*jEy - ny*jEx);
+            }
+        }
+    }
+}
+
 static void
 compute_fluxes(solver_state& state, const field& in, field& out)
 {
@@ -209,13 +380,41 @@ compute_fluxes(solver_state& state, const field& in, field& out)
     }
 }
 
+static void
+compute_fluxes_E(solver_state& state, const field& in, field& out)
+{
+    compute_field_jumps_E(state, in);
+    compute_fluxes_planar_E(state);
+    
+    for (auto& ed : state.eds)
+    {
+        compute_flux_lifting(ed, state.fluxes.Hx, out.Hx);
+        compute_flux_lifting(ed, state.fluxes.Hy, out.Hy);
+        compute_flux_lifting(ed, state.fluxes.Hz, out.Hz);
+    }
+}
+
+static void
+compute_fluxes_H(solver_state& state, const field& in, field& out)
+{
+    compute_field_jumps_H(state, in);
+    compute_fluxes_planar_H(state);
+    
+    for (auto& ed : state.eds)
+    {
+        compute_flux_lifting(ed, state.fluxes.Ex, out.Ex);
+        compute_flux_lifting(ed, state.fluxes.Ey, out.Ey);
+        compute_flux_lifting(ed, state.fluxes.Ez, out.Ez);
+    }
+}
+
 static void
 compute_euler_update(solver_state& state, const field& y, const field& k, double dt, field& out)
 {
     #pragma omp parallel for
     for (size_t i = 0; i < out.num_dofs; i++)
     {
-        auto CR = 1 - dt*state.matparams.sigma_over_epsilon[i];
+        auto CR = 1.0 - dt*state.matparams.sigma_over_epsilon[i];
         auto CC = dt*state.matparams.inv_epsilon[i];
         out.Ex[i] = y.Ex[i]*CR + k.Ex[i]*CC;
         out.Ey[i] = y.Ey[i]*CR + k.Ey[i]*CC;
@@ -238,7 +437,7 @@ compute_rk4_weighted_sum(solver_state& state, const field& in, double dt, field&
     #pragma omp parallel for
     for (size_t i = 0; i < out.num_dofs; i++)
     {
-        auto CR = 1 - dt*state.matparams.sigma_over_epsilon[i];
+        auto CR = 1.0 - dt*state.matparams.sigma_over_epsilon[i];
         auto CC = dt*state.matparams.inv_epsilon[i];
         out.Ex[i] = in.Ex[i]*CR + CC*(state.k1.Ex[i] + 2*state.k2.Ex[i] + 2*state.k3.Ex[i] + state.k4.Ex[i])/6;
         out.Ey[i] = in.Ey[i]*CR + CC*(state.k1.Ey[i] + 2*state.k2.Ey[i] + 2*state.k3.Ey[i] + state.k4.Ey[i])/6;
@@ -262,26 +461,65 @@ apply_operator(solver_state& state, const field& curr, field& next)
     compute_fluxes(state, curr, next);
 }
 
+static void
+leapfrog(solver_state& state)
+{
+    auto dt = state.delta_t;
+    compute_curls_H(state, state.emf_curr, state.tmp);
+    compute_fluxes_H(state, state.emf_curr, state.tmp);
+    #pragma omp parallel for
+    for (size_t i = 0; i < state.emf_next.num_dofs; i++)
+    {
+        auto CRM = 1.0 - 0.5*dt*state.matparams.sigma_over_epsilon[i];
+        auto CRP = 1.0 + 0.5*dt*state.matparams.sigma_over_epsilon[i];
+        auto CR = CRM/CRP;
+        auto CC = dt*state.matparams.inv_epsilon[i]/CRP;
+        state.emf_next.Ex[i] = state.emf_curr.Ex[i]*CR + state.tmp.Ex[i]*CC; 
+        state.emf_next.Ey[i] = state.emf_curr.Ey[i]*CR + state.tmp.Ey[i]*CC; 
+        state.emf_next.Ez[i] = state.emf_curr.Ez[i]*CR + state.tmp.Ez[i]*CC; 
+    }
+
+    compute_curls_E(state, state.emf_next, state.tmp);
+    compute_fluxes_E(state, state.emf_next, state.tmp);
+    #pragma omp parallel for
+    for (size_t i = 0; i < state.emf_next.num_dofs; i++)
+    {
+        auto CC = dt*state.matparams.inv_mu[i];
+        state.emf_next.Hx[i] = state.emf_curr.Hx[i] + state.tmp.Hx[i]*CC; 
+        state.emf_next.Hy[i] = state.emf_curr.Hy[i] + state.tmp.Hy[i]*CC; 
+        state.emf_next.Hz[i] = state.emf_curr.Hz[i] + state.tmp.Hz[i]*CC; 
+    }
+}
+
 void
-timestep(solver_state& state)
+timestep(solver_state& state, time_integrator_type ti)
 {
-    /*
-    apply_operator(state, state.emf_curr, state.tmp);
-    compute_euler_update(state, state.emf_curr, state.tmp, state.delta_t, state.emf_next);
-    */
+    if (ti == time_integrator_type::EULER)
+    {
+        apply_operator(state, state.emf_curr, state.tmp);
+        compute_euler_update(state, state.emf_curr, state.tmp, state.delta_t, state.emf_next);
+    }
 
-    apply_operator(state, state.emf_curr, state.k1);
+    if (ti == time_integrator_type::RK4)
+    {    
+        apply_operator(state, state.emf_curr, state.k1);
 
-    compute_euler_update(state, state.emf_curr, state.k1, 0.5*state.delta_t, state.tmp);
-    apply_operator(state, state.tmp, state.k2);
+        compute_euler_update(state, state.emf_curr, state.k1, 0.5*state.delta_t, state.tmp);
+        apply_operator(state, state.tmp, state.k2);
 
-    compute_euler_update(state, state.emf_curr, state.k2, 0.5*state.delta_t, state.tmp);
-    apply_operator(state, state.tmp, state.k3);
+        compute_euler_update(state, state.emf_curr, state.k2, 0.5*state.delta_t, state.tmp);
+        apply_operator(state, state.tmp, state.k3);
         
-    compute_euler_update(state, state.emf_curr, state.k3, state.delta_t, state.tmp);
-    apply_operator(state, state.tmp, state.k4);
+        compute_euler_update(state, state.emf_curr, state.k3, state.delta_t, state.tmp);
+        apply_operator(state, state.tmp, state.k4);
+
+        compute_rk4_weighted_sum(state, state.emf_curr, state.delta_t, state.emf_next);
+    }
 
-    compute_rk4_weighted_sum(state, state.emf_curr, state.delta_t, state.emf_next);
+    if (ti == time_integrator_type::LEAPFROG)
+    {
+        leapfrog(state);
+    }
 
     state.curr_time += state.delta_t;
     state.curr_timestep += 1;

src/maxwell_gpu.cpp

@@ -204,28 +204,32 @@ compute_rk4_weighted_sum(solver_state_gpu& state, const field_gpu& in,
 }
 
 void
-timestep(solver_state_gpu& state)
+timestep(solver_state_gpu& state, time_integrator_type ti)
 {
     //timecounter tc;
     //tc.tic();
 
-    /*
-    apply_operator(state, state.emf_curr, state.tmp);
-    compute_euler_update(state, state.emf_curr, state.tmp, state.delta_t, state.emf_next);
-    */
+    if (ti == time_integrator_type::EULER)
+    {
+        apply_operator(state, state.emf_curr, state.tmp);
+        compute_euler_update(state, state.emf_curr, state.tmp, state.delta_t, state.emf_next);
+    }
 
-    apply_operator(state, state.emf_curr, state.k1);
+    if (ti == time_integrator_type::RK4)
+    {
+        apply_operator(state, state.emf_curr, state.k1);
 
-    compute_euler_update(state, state.emf_curr, state.k1, 0.5*state.delta_t, state.tmp);
-    apply_operator(state, state.tmp, state.k2);
+        compute_euler_update(state, state.emf_curr, state.k1, 0.5*state.delta_t, state.tmp);
+        apply_operator(state, state.tmp, state.k2);
 
-    compute_euler_update(state, state.emf_curr, state.k2, 0.5*state.delta_t, state.tmp);
-    apply_operator(state, state.tmp, state.k3);
+        compute_euler_update(state, state.emf_curr, state.k2, 0.5*state.delta_t, state.tmp);
+        apply_operator(state, state.tmp, state.k3);
         
-    compute_euler_update(state, state.emf_curr, state.k3, state.delta_t, state.tmp);
-    apply_operator(state, state.tmp, state.k4);
+        compute_euler_update(state, state.emf_curr, state.k3, state.delta_t, state.tmp);
+        apply_operator(state, state.tmp, state.k4);
 
-    compute_rk4_weighted_sum(state, state.emf_curr, state.delta_t, state.emf_next);
+        compute_rk4_weighted_sum(state, state.emf_curr, state.delta_t, state.emf_next);
+    }
 
     state.curr_time += state.delta_t;
     state.curr_timestep += 1;

src/maxwell_solver.cpp

@@ -80,11 +80,13 @@ void test_it(const model& mod, State& state, maxwell::parameter_loader& mpl)
     omp_set_num_threads(4);
 #endif
 
+    auto ti = mpl.sim_timeIntegrator();
+
     prepare_sources(mod, state, mpl);
     for(size_t i = 0; i < num_timesteps; i++)
     {
         mpl.call_timestep_callback(i);
-        timestep(state);
+        timestep(state, ti);
         do_sources(mod, state, mpl);
         
 

src/param_loader.cpp

@@ -22,6 +22,7 @@ parameter_loader_base::init(void)
     //lua["sim"]["dt"];
     //lua["sim"]["timesteps"];
     //lua["sim"]["gmsh_model"];
+    lua["sim"]["time_integrator"] = "rk4";
 
     lua["materials"] = lua.create_table();
     lua["bndconds"] = lua.create_table();
@@ -210,3 +211,26 @@ parameter_loader_base::source_was_cleared(void)
 {
     m_source_changed_state = false;
 }
+
+time_integrator_type
+parameter_loader_base::sim_timeIntegrator(void) const
+{
+    std::string ti = lua["sim"]["time_integrator"];
+
+    if (ti == "rk4")
+        return time_integrator_type::RK4;
+
+    if (ti == "leapfrog")
+        return time_integrator_type::LEAPFROG;
+
+    if (ti == "euler")
+    {
+        std::cout << "[CONFIG] warning: Euler time integrator is only ";
+        std::cout << "for testing purposes" << std::endl;
+        return time_integrator_type::EULER;
+    }
+
+    std::cout << "[CONFIG] warning: invalid time integrator '";
+    std::cout << ti << "'" << std::endl;
+    return time_integrator_type::RK4;
+}
\ No newline at end of file