Added upwind to leapfrog. Looks like it works well. Sources still need to be fixed.

CMakeLists.txt

@@ -108,8 +108,6 @@ if (OPT_ENABLE_GPU_SOLVER)
     ## Find HIP stuff
     find_package(HIP QUIET MODULE)
     if (HIP_FOUND)
-        execute_process(COMMAND sh -c "/opt/rocm/hip/bin/hipify-perl ${CMAKE_SOURCE_DIR}/src/kernels_cuda.cu > ${CMAKE_SOURCE_DIR}/src/hipified_kernels_cuda.cpp")
-        execute_process(COMMAND sh -c "/opt/rocm/hip/bin/hipify-perl ${CMAKE_SOURCE_DIR}/src/maxwell_kernels_cuda.cu > ${CMAKE_SOURCE_DIR}/src/hipified_maxwell_kernels_cuda.cpp")
         add_definitions(-DHAVE_HIP)
         include_directories("/opt/rocm/include")
         set(CMAKE_HIP_LINK_EXECUTABLE "${HIP_HIPCC_CMAKE_LINKER_HELPER} ${HCC_HOME} <FLAGS> <CMAKE_CXX_LINK_FLAGS> <LINK_FLAGS> <OBJECTS> -o <TARGET> <LINK_LIBRARIES>")
@@ -123,6 +121,8 @@ if (OPT_ENABLE_GPU_SOLVER)
         set(HIP_NVCC_FLAGS "-ccbin g++-8" ${HIP_NVCC_FLAGS})
         option(OPT_USE_HIP "Use HIP GPU code path" OFF)
         if (OPT_USE_HIP)
+            execute_process(COMMAND sh -c "/opt/rocm/hip/bin/hipify-perl ${CMAKE_SOURCE_DIR}/src/kernels_cuda.cu > ${CMAKE_SOURCE_DIR}/src/hipified_kernels_cuda.cpp")
+            execute_process(COMMAND sh -c "/opt/rocm/hip/bin/hipify-perl ${CMAKE_SOURCE_DIR}/src/maxwell_kernels_cuda.cu > ${CMAKE_SOURCE_DIR}/src/hipified_maxwell_kernels_cuda.cpp")
             add_definitions(-DGPU_USE_HIP)
             add_definitions(-D__HIP_PLATFORM_HCC__)
             set(LINK_LIBS ${LINK_LIBS})
@@ -163,8 +163,8 @@ if (COMPILER_IS_CLANG)
         -Wno-padded -Wno-shorten-64-to-32 -Wno-sign-conversion -Wno-old-style-cast \
         -Wno-sign-compare -Wno-c99-extensions -Wno-extra-semi-stmt -Wno-source-uses-openmp")
 else()
-    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wno-sign-compare -Wshadow \
-        -Wno-unknown-pragmas")
+    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wextra -Wno-sign-compare -Wshadow \
+        -Wno-unknown-pragmas -Wno-unused-parameter")
 endif()
 set(CMAKE_CXX_FLAGS_DEBUG "-O0 -g -fpermissive")
 set(CMAKE_CXX_FLAGS_RELEASE "-O3 -g -DNDEBUG -fpermissive")

doc/changelog.md

@@ -18,5 +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
+- Preliminary implementation of leapfrog integration on CPU.
 

doc/lua_api.md

@@ -20,7 +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.
+- `sim.time_integrator` (string): `"rk4"`, `"leapfrog"` or `"euler"`. Euler is only for test purporses, as it is numerically instable. Leapfrog works on CPU for now, use it if you have dissipative materials.
 
 ### Postprocessing general variables
 

include/param_loader.h

@@ -46,6 +46,7 @@ public:
     std::string     sim_gmshmodel(void) const;
     bool            sim_usegpu(void) const;
     time_integrator_type    sim_timeIntegrator(void) const;
+    std::string             sim_timeIntegratorName(void) const;
     size_t          postpro_siloOutputRate(void) const;
     size_t          postpro_cyclePrintRate(void) const;
 

src/maxwell_cpu.cpp

@@ -287,7 +287,6 @@ compute_fluxes_planar(solver_state& state)
     }
 }
 
-/* In leapfrog we use centered fluxes for the moment */
 static void
 compute_fluxes_planar_E(solver_state& state)
 {
@@ -310,22 +309,27 @@ compute_fluxes_planar_E(solver_state& state)
                 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 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;
+                auto ndotE = nx*jEx + ny*jEy + nz*jEz;
+
+                auto aE = face_det * state.matparams.aE[gofs];
+                auto bE = face_det * state.matparams.bE[gofs];
 
                 /* 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);
+                state.fluxes.Ex[gofs] = aE*(nz*jHy - ny*jHz) + bE*(ndotE*nx - jEx);
+                state.fluxes.Ey[gofs] = aE*(nx*jHz - nz*jHx) + bE*(ndotE*ny - jEy);
+                state.fluxes.Ez[gofs] = aE*(ny*jHx - nx*jHy) + bE*(ndotE*nz - jEz);
             }
         }
     }
 }
 
-/* In leapfrog we use centered fluxes for the moment */
 static void
 compute_fluxes_planar_H(solver_state& state)
 {
@@ -351,13 +355,19 @@ compute_fluxes_planar_H(solver_state& state)
                 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 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 ndotH = nx*jHx + ny*jHy + nz*jHz;
 
-                auto aH = face_det * 0.5;
+                auto aH = face_det * state.matparams.aH[gofs];
+                auto bH = face_det * state.matparams.bH[gofs];
 
                 /* 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);
+                state.fluxes.Hx[gofs] = aH*(ny*jEz - nz*jEy) + bH*(ndotH*nx - jHx);
+                state.fluxes.Hy[gofs] = aH*(nz*jEx - nx*jEz) + bH*(ndotH*ny - jHy);
+                state.fluxes.Hz[gofs] = aH*(nx*jEy - ny*jEx) + bH*(ndotH*nz - jHz);
             }
         }
     }

src/maxwell_solver.cpp

@@ -75,12 +75,16 @@ void test_it(const model& mod, State& state, maxwell::parameter_loader& mpl)
     auto num_timesteps = mpl.sim_timesteps();
     auto silo_output_rate = mpl.postpro_siloOutputRate();
     auto cycle_print_rate = mpl.postpro_cyclePrintRate();
+    auto ti = mpl.sim_timeIntegrator();
 
 #ifdef _OPENMP
     omp_set_num_threads(4);
 #endif
 
-    auto ti = mpl.sim_timeIntegrator();
+    std::cout << "    BEGINNING SIMULATION" << std::endl;
+    std::cout << "I will do " << num_timesteps << " of " << state.delta_t;
+    std::cout << "s each." << std::endl;
+    std::cout << "Time integrator: " << mpl.sim_timeIntegratorName() << std::endl;
 
     prepare_sources(mod, state, mpl);
     for(size_t i = 0; i < num_timesteps; i++)

src/param_loader.cpp

@@ -234,3 +234,11 @@ parameter_loader_base::sim_timeIntegrator(void) const
     std::cout << ti << "'" << std::endl;
     return time_integrator_type::RK4;
 }
+
+std::string
+parameter_loader_base::sim_timeIntegratorName(void) const
+{
+    std::string ti = lua["sim"]["time_integrator"];
+    return ti;
+}
+