diff --git a/include/maxwell_interface.h b/include/maxwell_interface.h
index 7e55f0fe5d711bff5fbf1330d00e26401eb75e84..35defba4a86db21a279fecfe47accb61b30f1c0a 100644
--- a/include/maxwell_interface.h
+++ b/include/maxwell_interface.h
@@ -87,6 +87,7 @@ struct field
 
     size_t  num_dofs;
 
+    void    zero(void);
     void    resize(size_t);
 };
 
@@ -107,6 +108,7 @@ struct pinned_field
     pinned_field(const pinned_field&) = delete;
     pinned_field& operator=(const pinned_field&) = delete;
     
+    void    zero(void);
     void    resize(size_t);
 };
 #endif /* ENABLE_GPU_SOLVER */
@@ -290,6 +292,9 @@ void init_matparams(const model&, solver_state&, const parameter_loader&);
 //void apply_operator(solver_state&, const field&, field&);
 void export_to_silo(const model&, const solver_state&, const std::string&);
 void timestep(solver_state&);
+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&);
 
 #ifdef ENABLE_GPU_SOLVER
 void init_from_model(const model&, solver_state_gpu&);
@@ -297,6 +302,9 @@ void init_matparams(const model&, solver_state_gpu&, const parameter_loader&);
 //void apply_operator(solver_state_gpu&, const field_gpu&, field_gpu&);
 void export_to_silo(const model&, const solver_state_gpu&, const std::string&);
 void timestep(solver_state_gpu&);
+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&);
 #endif /* ENABLE_GPU_SOLVER */
 
 void init_boundary_conditions(const model&, const parameter_loader&, vecxd&);
diff --git a/include/param_loader.h b/include/param_loader.h
index 20630eeff1c02a95b997a4db6431e9b59158d4f7..d33fe330dae39bba2c85dcabecb420d8ee18601a 100644
--- a/include/param_loader.h
+++ b/include/param_loader.h
@@ -9,6 +9,11 @@
 
 class parameter_loader_base
 {
+    bool        m_bnd_sources_active;
+    bool        m_ifc_sources_active;
+    bool        m_vol_sources_active;
+    bool        m_source_changed_state;
+
 protected:
     sol::state lua;
 
@@ -30,6 +35,19 @@ public:
     bool            sim_usegpu(void) const;
     size_t          postpro_siloOutputRate(void) const;
     size_t          postpro_cyclePrintRate(void) const;
+
+    void            enable_boundary_sources(bool);
+    void            enable_interface_sources(bool);
+    void            enable_volume_sources(bool);
+
+    bool            boundary_sources_enabled(void) const;
+    bool            interface_sources_enabled(void) const;
+    bool            volume_sources_enabled(void) const;
+
+    bool            source_has_changed_state(void) const;
+    void            source_was_cleared(void);
+
+    void            call_timestep_callback(size_t);
 };
 
 
diff --git a/src/maxwell_cpu.cpp b/src/maxwell_cpu.cpp
index ee56a2a75bd5afb404ab81578ca903b7e7e14752..8da503207ae6f71b8e08987848996a92fe575b0e 100644
--- a/src/maxwell_cpu.cpp
+++ b/src/maxwell_cpu.cpp
@@ -2,6 +2,7 @@
 
 #include "gmsh_io.h"
 #include "maxwell_interface.h"
+#include "maxwell_common.h"
 #include "kernels_cpu.h"
 #include "silo_output.hpp"
 #include "timecounter.h"
@@ -343,6 +344,36 @@ timestep(solver_state& state)
     state.curr_timestep += 1;
 }
 
+void
+do_sources(const model& mod, maxwell::solver_state& state,
+    maxwell::parameter_loader& mpl)
+{
+    if ( mpl.source_has_changed_state() )
+    {
+        state.bndsrcs.zero();
+        mpl.source_was_cleared();
+    }
+
+    if ( mpl.boundary_sources_enabled() )
+        maxwell::eval_boundary_sources(mod, mpl, state, state.bndsrcs);
+    
+    if ( mpl.interface_sources_enabled() )
+        maxwell::eval_interface_sources(mod, mpl, state, state.bndsrcs);
+}
+
+void
+prepare_sources(const model& mod, maxwell::solver_state& state,
+    maxwell::parameter_loader& mpl)
+{
+    do_sources(mod, state, mpl);
+}
+
+void
+swap(maxwell::solver_state& state)
+{
+    std::swap(state.emf_curr, state.emf_next);
+}
+
 void
 export_to_silo(const model& mod, const solver_state& state, const std::string& filename)
 {
diff --git a/src/maxwell_gpu.cpp b/src/maxwell_gpu.cpp
index 98d3984d09316a65a458a7edbc07e11722b694a1..589dbce355cb03b396ea6b00be046a811c485939 100644
--- a/src/maxwell_gpu.cpp
+++ b/src/maxwell_gpu.cpp
@@ -2,6 +2,7 @@
 
 #include "gmsh_io.h"
 #include "maxwell_interface.h"
+#include "maxwell_common.h"
 #include "kernels_cpu.h"
 #include "kernels_gpu.h"
 #include "silo_output.hpp"
@@ -233,6 +234,46 @@ timestep(solver_state_gpu& state)
     //std::cout << "Timestep " << state.curr_timestep << ", " << dofs_per_sec << " DOFs/s" << std::endl;
 }
 
+void
+prepare_sources(const model& mod, maxwell::solver_state_gpu& state,
+    maxwell::parameter_loader& mpl)
+{
+    if ( mpl.boundary_sources_enabled() )
+        maxwell::eval_boundary_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);
+}
+
+void
+do_sources(const model& mod, maxwell::solver_state_gpu& state,
+    maxwell::parameter_loader& mpl)
+{
+    if ( mpl.source_has_changed_state() )
+    {
+        state.bndsrcs_cpu.zero();
+        mpl.source_was_cleared();
+    }
+
+    auto be = mpl.boundary_sources_enabled();
+    auto ie = mpl.interface_sources_enabled();
+    auto ve = mpl.volume_sources_enabled();
+
+    if ( be )
+        maxwell::eval_boundary_sources(mod, mpl, state, state.bndsrcs_cpu);
+
+    if ( be or ie or ve )
+        state.bndsrcs_buf.copyin(state.bndsrcs_cpu, state.memcpy_stream);
+}
+
+void
+swap(maxwell::solver_state_gpu& state)
+{
+    cudaDeviceSynchronize();
+    std::swap(state.bndsrcs, state.bndsrcs_buf);
+    std::swap(state.emf_curr, state.emf_next);
+}
 
 void
 export_to_silo(const model& mod, const solver_state_gpu& state, const std::string& filename)
diff --git a/src/maxwell_interface.cpp b/src/maxwell_interface.cpp
index b1e311e4f108ff02440657e48d56b57fdbc8d49e..4a72a54ad01c5fdc3eb170af0af422752302ead2 100644
--- a/src/maxwell_interface.cpp
+++ b/src/maxwell_interface.cpp
@@ -2,6 +2,17 @@
 
 namespace maxwell {
 
+void
+field::zero()
+{
+    Ex = vecxd::Zero(num_dofs);
+    Ey = vecxd::Zero(num_dofs);
+    Ez = vecxd::Zero(num_dofs);
+    Hx = vecxd::Zero(num_dofs);
+    Hy = vecxd::Zero(num_dofs);
+    Hz = vecxd::Zero(num_dofs);
+}
+
 void
 field::resize(size_t p_num_dofs)
 {
@@ -22,6 +33,17 @@ pinned_field::pinned_field()
       num_dofs(0)
 {}
 
+void
+pinned_field::zero()
+{
+    memset(Ex, '\0', num_dofs*sizeof(double));
+    memset(Ey, '\0', num_dofs*sizeof(double));
+    memset(Ez, '\0', num_dofs*sizeof(double));
+    memset(Hx, '\0', num_dofs*sizeof(double));
+    memset(Hy, '\0', num_dofs*sizeof(double));
+    memset(Hz, '\0', num_dofs*sizeof(double));
+}
+
 void
 pinned_field::resize(size_t sz)
 {
diff --git a/src/maxwell_solver.cpp b/src/maxwell_solver.cpp
index 9dda72de647cbd9a2b273d5dfc4cbe212613687f..2aab09f72faffb75dc80cbf508ba1bfeb52d7e29 100644
--- a/src/maxwell_solver.cpp
+++ b/src/maxwell_solver.cpp
@@ -7,7 +7,6 @@
 #include "gmsh.h"
 #include "silo_output.hpp"
 #include "maxwell_interface.h"
-#include "maxwell_common.h"
 #include "param_loader.h"
 
 #if 0
@@ -60,57 +59,9 @@ void initialize_solver(const model& mod, State& state, const maxwell::parameter_
     init_matparams(mod, state, mpl);
 }
 
-void
-do_sources(const model& mod, maxwell::solver_state& state,
-    const maxwell::parameter_loader& mpl)
-{
-    maxwell::eval_boundary_sources(mod, mpl, state, state.bndsrcs);
-    maxwell::eval_interface_sources(mod, mpl, state, state.bndsrcs);
-}
-
-void
-prepare_sources(const model& mod, maxwell::solver_state& state,
-    const maxwell::parameter_loader& mpl)
-{
-    do_sources(mod, state, mpl);
-}
-
-void
-swap(maxwell::solver_state& state)
-{
-    std::swap(state.emf_curr, state.emf_next);
-}
-
-#ifdef ENABLE_GPU_SOLVER
-void
-prepare_sources(const model& mod, maxwell::solver_state_gpu& state,
-    const maxwell::parameter_loader& mpl)
-{
-    maxwell::eval_boundary_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);
-}
-
-void
-do_sources(const model& mod, maxwell::solver_state_gpu& state,
-    const maxwell::parameter_loader& mpl)
-{
-    maxwell::eval_boundary_sources(mod, mpl, state, state.bndsrcs_cpu);
-    state.bndsrcs_buf.copyin(state.bndsrcs_cpu, state.memcpy_stream);
-}
-
-void
-swap(maxwell::solver_state_gpu& state)
-{
-    cudaDeviceSynchronize();
-    std::swap(state.bndsrcs, state.bndsrcs_buf);
-    std::swap(state.emf_curr, state.emf_next);
-}
-#endif
 
 template<typename State>
-void test_it(const model& mod, State& state, const maxwell::parameter_loader& mpl)
+void test_it(const model& mod, State& state, maxwell::parameter_loader& mpl)
 {
     initialize_solver(mod, state, mpl);
 
@@ -126,6 +77,7 @@ void test_it(const model& mod, State& state, const maxwell::parameter_loader& mp
     prepare_sources(mod, state, mpl);
     for(size_t i = 0; i < num_timesteps; i++)
     {
+        mpl.call_timestep_callback(i);
         timestep(state);
         do_sources(mod, state, mpl);
         std::stringstream ss;
diff --git a/src/param_loader.cpp b/src/param_loader.cpp
index b6576fb87a66c53147a901dc9aee4e5ebee0793d..034b0b728580a8e956b89a111dda61f783384ca4 100644
--- a/src/param_loader.cpp
+++ b/src/param_loader.cpp
@@ -2,6 +2,8 @@
 #include "param_loader.h"
 
 parameter_loader_base::parameter_loader_base()
+    : m_bnd_sources_active(true), m_ifc_sources_active(true),
+      m_vol_sources_active(true)
 {
     lua.open_libraries(sol::lib::base, sol::lib::math);
     init();
@@ -28,6 +30,18 @@ parameter_loader_base::init(void)
     lua["postpro"] = lua.create_table();
     lua["postpro"]["silo_output_rate"] = 0;
     lua["postpro"]["cycle_print_rate"] = 10;
+
+    lua.set_function("enable_boundary_sources",
+                     &parameter_loader_base::enable_boundary_sources,
+                     this);
+
+    lua.set_function("enable_interface_sources",
+                     &parameter_loader_base::enable_interface_sources,
+                     this);
+
+    lua.set_function("enable_volume_sources",
+                     &parameter_loader_base::enable_volume_sources,
+                     this);
 }
 
 bool
@@ -131,3 +145,68 @@ parameter_loader_base::postpro_siloOutputRate(void) const
     return lua["postpro"]["silo_output_rate"];
 }
 
+void
+parameter_loader_base::enable_boundary_sources(bool enable)
+{
+    m_bnd_sources_active = enable;
+    m_source_changed_state = true;
+}
+
+void
+parameter_loader_base::enable_interface_sources(bool enable)
+{
+    m_ifc_sources_active = enable;
+    m_source_changed_state = true;
+}
+
+void
+parameter_loader_base::enable_volume_sources(bool enable)
+{
+    m_vol_sources_active = enable;
+    m_source_changed_state = true;
+}
+
+bool
+parameter_loader_base::boundary_sources_enabled(void) const
+{
+    return m_bnd_sources_active;
+}
+
+bool
+parameter_loader_base::interface_sources_enabled(void) const
+{
+    return m_ifc_sources_active;
+}
+
+bool
+parameter_loader_base::volume_sources_enabled(void) const
+{
+    return m_vol_sources_active;
+}
+
+void
+parameter_loader_base::call_timestep_callback(size_t timestep)
+{
+    if ( lua["on_timestep"].valid() )
+    {
+        try {
+            lua["on_timestep"](timestep);
+        }
+        catch (...)
+        {
+            std::cout << "An error was thrown calling on_timestep()" << std::endl;
+        }
+    }
+}
+
+bool
+parameter_loader_base::source_has_changed_state(void) const
+{
+    return m_source_changed_state;
+}
+
+void
+parameter_loader_base::source_was_cleared(void)
+{
+    m_source_changed_state = false;
+}