diff --git a/makefile b/makefile
index c2befaa55..fb1d88f5e 100644
--- a/makefile
+++ b/makefile
@@ -43,7 +43,7 @@ make solvers
 	 Builds each solver executable.
 make {solver}
 	 Builds a solver executable,
-	 solver can be acoustics/advection/bns/cns/elliptic/fokkerPlanck/gradient/ins.
+	 solver can be acoustics/advection/bns/cns/elliptic/fokkerPlanck/gradient/ins/maxwell.
 make clean
 	 Cleans all solver executables, libraries, and object files.
 make clean-{solver}
@@ -64,7 +64,7 @@ Can use "make verbose=true" for verbose output.
 endef
 
 ifeq (,$(filter solvers \
-				acoustics advection bns cns elliptic fokkerPlanck gradient ins \
+				acoustics advection bns cns elliptic fokkerPlanck gradient ins maxwell \
 				lib clean clean-kernels \
 				realclean info help test,$(MAKECMDGOALS)))
 ifneq (,$(MAKECMDGOALS))
@@ -88,12 +88,12 @@ LIBP_CORE_LIBS=timeStepper linearSolver parAlmond mesh ogs linAlg core
 SOLVER_DIR   =${LIBP_DIR}/solvers
 
 .PHONY: all solvers libp_libs \
-			acoustics advection bns lbs cns elliptic fokkerPlanck gradient ins \
+			acoustics advection bns lbs cns elliptic fokkerPlanck gradient ins maxwell\
 			clean clean-libs realclean help info
 
 all: solvers
 
-solvers: acoustics advection bns lbs cns elliptic fokkerPlanck gradient ins
+solvers: acoustics advection bns lbs cns elliptic fokkerPlanck gradient ins maxwell
 
 libp_libs:
 ifneq (,${verbose})
@@ -174,10 +174,19 @@ else
 	@${MAKE} -C ${SOLVER_DIR}/$(@F) --no-print-directory
 endif
 
+maxwell: libp_libs
+ifneq (,${verbose})
+	${MAKE} -C ${SOLVER_DIR}/$(@F) verbose=${verbose}
+else
+	@printf "%b" "$(SOL_COLOR)Building $(@F) solver$(NO_COLOR)\n";
+	@${MAKE} -C ${SOLVER_DIR}/$(@F) --no-print-directory
+endif
+
+
 #cleanup
 clean: clean-acoustics clean-advection clean-bns clean-lbs clean-cns \
 	   clean-elliptic clean-fokkerPlanck clean-gradient clean-ins \
-	   clean-libs
+	   clean-maxwell clean-libs
 
 clean-acoustics:
 	${MAKE} -C ${SOLVER_DIR}/acoustics clean
@@ -206,6 +215,9 @@ clean-gradient:
 clean-ins:
 	${MAKE} -C ${SOLVER_DIR}/ins clean
 
+clean-maxwell:
+	${MAKE} -C ${SOLVER_DIR}/maxwell clean
+
 clean-libs:
 	${MAKE} -C ${LIBP_LIBS_DIR} clean
 
diff --git a/solvers/maxwell/data/maxwellGaussian2D.h b/solvers/maxwell/data/maxwellGaussian2D.h
new file mode 100644
index 000000000..965387690
--- /dev/null
+++ b/solvers/maxwell/data/maxwellGaussian2D.h
@@ -0,0 +1,45 @@
+/*
+
+The MIT License (MIT)
+
+Copyright (c) 2017-2022 Tim Warburton, Noel Chalmers, Jesse Chan, Ali Karakus
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+
+*/
+
+
+// Boundary conditions
+/* PEC 1 */
+#define maxwellDirichletConditions2D(bc, t, x, y, nx, ny, HxM, HyM, EzM, HxB, HyB, EzB) \
+  {									\
+    if(bc==1){								\
+      *(HxB) = HxM;							\
+      *(HyB) = HyM;							\
+      *(EzB) = -EzM;							\
+    }									\
+  }
+
+// Initial conditions
+#define maxwellInitialConditions2D(t, x, y, Hx, Hy, Ez) \
+  {							\
+    *(Hx) = 0.0;					\
+    *(Hy) = 0.0;					\
+    *(Ez) = exp(-3*(x*x+y*y));				\
+  }
diff --git a/solvers/maxwell/data/maxwellGaussian3D.h b/solvers/maxwell/data/maxwellGaussian3D.h
new file mode 100644
index 000000000..c834ee9c8
--- /dev/null
+++ b/solvers/maxwell/data/maxwellGaussian3D.h
@@ -0,0 +1,52 @@
+/*
+
+The MIT License (MIT)
+
+Copyright (c) 2017-2022 Tim Warburton, Noel Chalmers, Jesse Chan, Ali Karakus
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+
+*/
+
+
+
+// Boundary conditions
+/* PEC=1 */
+#define maxwellDirichletConditions3D(bc, t, x, y, z, nx, ny, nz, HxM, HyM, HzM, ExM, EyM, EzM, HxB, HyB, HzB, ExB, EyB, EzB) \
+  {									\
+    if(bc==1){								\
+      *(HxB) = HxM;							\
+      *(HyB) = HyM;							\
+      *(HzB) = HzM;							\
+      *(ExB) = -ExM;							\
+      *(EyB) = -EyM;							\
+      *(EzB) = -EzM;							\
+    }									\
+  }
+
+// Initial conditions
+#define maxwellInitialConditions3D(t, x, y, z, Hx, Hy, Hz, Ex, Ey, Ez) \
+  {							   \
+    *(Hx) = 0.0;					   \
+    *(Hy) = 0.0;					   \
+    *(Hz) = 0.0;					   \
+    *(Ex) = 0.0;					   \
+    *(Ey) = 0.0;					   \
+    *(Ez) = exp(-3*(x*x+y*y+z*z));			   \
+  }
diff --git a/solvers/maxwell/data/maxwellMode2D.h b/solvers/maxwell/data/maxwellMode2D.h
new file mode 100644
index 000000000..a45cd4910
--- /dev/null
+++ b/solvers/maxwell/data/maxwellMode2D.h
@@ -0,0 +1,48 @@
+/*
+
+The MIT License (MIT)
+
+Copyright (c) 2017-2022 Tim Warburton, Noel Chalmers, Jesse Chan, Ali Karakus
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+
+*/
+
+
+
+// Boundary conditions
+/* PEC=1 */
+#define maxwellDirichletConditions2D(bc, t, x, y, nx, ny, HxM, HyM,  EzM, HxB, HyB, EzB) \
+  {									\
+    if(bc==1){								\
+      *(HxB) = HxM;							\
+      *(HyB) = HyM;							\
+      *(EzB) = -EzM;							\
+    }									\
+  }
+
+// Initial conditions
+#define maxwellInitialConditions2D(t, x, y, Hx, Hy, Ez)			\
+  {									\
+    dfloat n = 1., m = 2.;						\
+    dfloat omega = M_PI*sqrt(n*n+m*m);					\
+    *(Hx) = -M_PI*n*sin(m*M_PI*x)*cos(n*M_PI*y)*sin(omega*t)/omega;	\
+    *(Hy) =  M_PI*m*cos(m*M_PI*x)*sin(n*M_PI*y)*sin(omega*t)/omega;	\
+    *(Ez) = sin(m*M_PI*x)*sin(n*M_PI*y)*cos(omega*t);			\
+  }
diff --git a/solvers/maxwell/data/maxwellMode3D.h b/solvers/maxwell/data/maxwellMode3D.h
new file mode 100644
index 000000000..12737c5d6
--- /dev/null
+++ b/solvers/maxwell/data/maxwellMode3D.h
@@ -0,0 +1,54 @@
+/*
+
+The MIT License (MIT)
+
+Copyright (c) 2017-2022 Tim Warburton, Noel Chalmers, Jesse Chan, Ali Karakus
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+
+*/
+
+
+
+// Boundary conditions
+/* PEC=1 */
+#define maxwellDirichletConditions3D(bc, t, x, y, z, nx, ny, nz, HxM, HyM, HzM, ExM, EyM, EzM, HxB, HyB, HzB, ExB, EyB, EzB) \
+  {									\
+    if(bc==1){								\
+      *(HxB) = HxM;							\
+      *(HyB) = HyM;							\
+      *(HzB) = HzM;							\
+      *(ExB) = -ExM;							\
+      *(EyB) = -EyM;							\
+      *(EzB) = -EzM;							\
+    }									\
+  }
+
+// Initial conditions
+#define maxwellInitialConditions3D(t, x, y, z, Hx, Hy, Hz, Ex, Ey, Ez)	\
+  {									\
+    dfloat n = 1., m = 2.;						\
+    const dfloat omega = M_PI*sqrt(n*n+m*m);				\
+    *(Hx) = -M_PI*n*sin(m*M_PI*x)*cos(n*M_PI*y)*sin(omega*t)/omega;	\
+    *(Hy) =  M_PI*m*cos(m*M_PI*x)*sin(n*M_PI*y)*sin(omega*t)/omega;	\
+    *(Hz) = 0;								\
+    *(Ex) = 0;								\
+    *(Ey) = 0;								\
+    *(Ez) = sin(m*M_PI*x)*sin(n*M_PI*y)*cos(omega*t);			\
+  }
diff --git a/solvers/maxwell/makefile b/solvers/maxwell/makefile
new file mode 100644
index 000000000..1f99d2f91
--- /dev/null
+++ b/solvers/maxwell/makefile
@@ -0,0 +1,175 @@
+#####################################################################################
+#
+#The MIT License (MIT)
+#
+#Copyright (c) 2017-2022 Tim Warburton, Noel Chalmers, Jesse Chan, Ali Karakus
+#
+#Permission is hereby granted, free of charge, to any person obtaining a copy
+#of this software and associated documentation files (the "Software"), to deal
+#in the Software without restriction, including without limitation the rights
+#to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+#copies of the Software, and to permit persons to whom the Software is
+#furnished to do so, subject to the following conditions:
+#
+#The above copyright notice and this permission notice shall be included in all
+#copies or substantial portions of the Software.
+#
+#THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+#IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+#FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+#AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+#LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+#OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+#SOFTWARE.
+#
+#####################################################################################
+
+define MAXWELL_HELP_MSG
+
+Maxwell solver makefile targets:
+
+   make maxwellMain (default)
+   make lib
+   make clean
+   make clean-libs
+   make clean-kernels
+   make realclean
+   make info
+   make help
+   make test
+
+Usage:
+
+make maxwellMain
+   Build maxwellMain executable.
+make lib
+   Build libmaxwell.a solver library.
+make clean
+   Clean the maxwellMain executable, library, and object files.
+make clean-libs
+   In addition to "make clean", also clean needed libraries.
+make clean-kernels
+   In addition to "make clean-libs", also cleans the cached OCCA kernels.
+make realclean
+   In addition to "make clean-kernels", also clean 3rd party libraries.
+make info
+   List directories and compiler flags in use.
+make help
+   Display this help message.
+make test
+   Run tests.
+Can use "make verbose=true" for verbose output.
+
+endef
+
+ifeq (,$(filter maxwellMain lib clean clean-libs clean-kernels \
+                realclean info help test,$(MAKECMDGOALS)))
+ifneq (,$(MAKECMDGOALS))
+$(error ${MAXWELL_HELP_MSG})
+endif
+endif
+
+ifndef LIBP_MAKETOP_LOADED
+ifeq (,$(wildcard ../../make.top))
+$(error cannot locate ${PWD}/../../make.top)
+else
+include ../../make.top
+endif
+endif
+
+#libraries
+MAXWELL_LIBP_LIBS=timeStepper mesh parAdogs ogs linAlg core
+
+#includes
+INCLUDES=${LIBP_INCLUDES} \
+				 -I.
+
+#defines
+DEFINES =${LIBP_DEFINES} \
+         -DLIBP_DIR='"${LIBP_DIR}"'
+
+#.cpp compilation flags
+MAXWELL_CXXFLAGS=${LIBP_CXXFLAGS} ${DEFINES} ${INCLUDES}
+
+#link libraries
+LIBS=-L${LIBP_LIBS_DIR} $(addprefix -l,$(MAXWELL_LIBP_LIBS)) \
+     ${LIBP_LIBS}
+
+#link flags
+LFLAGS=${MAXWELL_CXXFLAGS} ${LIBS}
+
+#object dependancies
+DEPS=$(wildcard *.hpp) \
+     $(wildcard $(LIBP_INCLUDE_DIR)/*.h) \
+     $(wildcard $(LIBP_INCLUDE_DIR)/*.hpp)
+
+SRC =$(wildcard src/*.cpp)
+
+OBJS=$(SRC:.cpp=.o)
+
+.PHONY: all lib libp_libs clean clean-libs \
+		clean-kernels realclean help info
+
+all: maxwellMain
+
+lib: libmaxwell.a
+
+libp_libs:
+ifneq (,${verbose})
+	${MAKE} -C ${LIBP_LIBS_DIR} $(MAXWELL_LIBP_LIBS) verbose=${verbose}
+else
+	@${MAKE} -C ${LIBP_LIBS_DIR} $(MAXWELL_LIBP_LIBS) --no-print-directory
+endif
+
+maxwellMain:$(OBJS) maxwellMain.o libp_libs
+ifneq (,${verbose})
+	$(LIBP_LD) -o maxwellMain maxwellMain.o $(OBJS) $(MESH_OBJS) $(LFLAGS)
+else
+	@printf "%b" "$(EXE_COLOR)Linking $(@F)$(NO_COLOR)\n";
+	@$(LIBP_LD) -o maxwellMain maxwellMain.o $(OBJS) $(MESH_OBJS) $(LFLAGS)
+endif
+
+libmaxwell.a: $(OBJS)
+ifneq (,${verbose})
+	ar -cr libmaxwell.a $(OBJS)
+else
+	@printf "%b" "$(LIB_COLOR)Building library $(@F)$(NO_COLOR)\n";
+	@ar -cr libmaxwell.a $(OBJS)
+endif
+
+# rule for .cpp files
+%.o: %.cpp $(DEPS) | libp_libs
+ifneq (,${verbose})
+	$(LIBP_CXX) -o $*.o -c $*.cpp $(MAXWELL_CXXFLAGS)
+else
+	@printf "%b" "$(OBJ_COLOR)Compiling $(@F)$(NO_COLOR)\n";
+	@$(LIBP_CXX) -o $*.o -c $*.cpp $(MAXWELL_CXXFLAGS)
+endif
+
+#cleanup
+clean:
+	rm -f src/*.o *.o maxwellMain libmaxwell.a
+
+clean-libs: clean
+	${MAKE} -C ${LIBP_LIBS_DIR} clean
+
+clean-kernels: clean-libs
+	rm -rf ${LIBP_DIR}/.occa/
+
+realclean: clean
+	${MAKE} -C ${LIBP_LIBS_DIR} realclean
+
+help:
+	$(info $(value MAXWELL_HELP_MSG))
+	@true
+
+info:
+	$(info OCCA_DIR  = $(OCCA_DIR))
+	$(info LIBP_DIR  = $(LIBP_DIR))
+	$(info LIBP_ARCH = $(LIBP_ARCH))
+	$(info CXXFLAGS  = $(MAXWELL_CXXFLAGS))
+	$(info LIBS      = $(LIBS))
+	@true
+
+test: maxwellMain
+	@${MAKE} -C $(LIBP_TEST_DIR) --no-print-directory  test-maxwell
diff --git a/solvers/maxwell/maxwell.hpp b/solvers/maxwell/maxwell.hpp
new file mode 100644
index 000000000..ab0966137
--- /dev/null
+++ b/solvers/maxwell/maxwell.hpp
@@ -0,0 +1,90 @@
+/*
+
+The MIT License (MIT)
+
+Copyright (c) 2017-2022 Tim Warburton, Noel Chalmers, Jesse Chan, Ali Karakus
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+
+*/
+
+#ifndef MAXWELL_HPP
+#define MAXWELL_HPP 1
+
+#include "core.hpp"
+#include "platform.hpp"
+#include "mesh.hpp"
+#include "solver.hpp"
+#include "timeStepper.hpp"
+#include "linAlg.hpp"
+
+#define DMAXWELL LIBP_DIR"/solvers/maxwell/"
+
+using namespace libp;
+
+class maxwellSettings_t: public settings_t {
+public:
+  maxwellSettings_t(comm_t _comm);
+  void report();
+  void parseFromFile(platformSettings_t& platformSettings,
+                     meshSettings_t& meshSettings,
+                     const std::string filename);
+};
+
+class maxwell_t: public solver_t {
+public:
+  mesh_t mesh;
+
+  int Nfields;
+
+  timeStepper_t timeStepper;
+
+  ogs::halo_t traceHalo;
+
+  memory<dfloat> q;
+  deviceMemory<dfloat> o_q;
+
+  kernel_t volumeKernel;
+  kernel_t surfaceKernel;
+
+  kernel_t initialConditionKernel;
+  kernel_t errorKernel;
+
+  maxwell_t() = default;
+  maxwell_t(platform_t &_platform, mesh_t &_mesh,
+              maxwellSettings_t& _settings) {
+    Setup(_platform, _mesh, _settings);
+  }
+
+  //setup
+  void Setup(platform_t& _platform, mesh_t& _mesh,
+             maxwellSettings_t& _settings);
+
+  void Run();
+
+  void Report(dfloat time, int tstep);
+
+  void PlotFields(memory<dfloat> Q, const std::string fileName);
+
+  void rhsf(deviceMemory<dfloat>& o_q, deviceMemory<dfloat>& o_rhs, const dfloat time);
+
+  dfloat MaxWaveSpeed();
+};
+
+#endif
diff --git a/solvers/maxwell/maxwellMain.cpp b/solvers/maxwell/maxwellMain.cpp
new file mode 100644
index 000000000..2f5b4d599
--- /dev/null
+++ b/solvers/maxwell/maxwellMain.cpp
@@ -0,0 +1,68 @@
+/*
+
+The MIT License (MIT)
+
+Copyright (c) 2017-2022 Tim Warburton, Noel Chalmers, Jesse Chan, Ali Karakus
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+
+*/
+
+#include "maxwell.hpp"
+
+int main(int argc, char **argv){
+
+  // start up MPI
+  Comm::Init(argc, argv);
+
+  LIBP_ABORT("Usage: ./maxwellMain setupfile", argc!=2);
+
+  { /*Scope so everything is destructed before MPI_Finalize */
+    comm_t comm(Comm::World().Dup());
+
+    //create default settings
+    platformSettings_t platformSettings(comm);
+    meshSettings_t meshSettings(comm);
+    maxwellSettings_t maxwellSettings(comm);
+
+    //load settings from file
+    maxwellSettings.parseFromFile(platformSettings, meshSettings,
+                              argv[1]);
+
+    // set up platform
+    platform_t platform(platformSettings);
+
+    platformSettings.report();
+    meshSettings.report();
+    maxwellSettings.report();
+
+    // set up mesh
+    mesh_t mesh(platform, meshSettings, comm);
+
+    // set up maxwell solver
+    maxwell_t maxwell(platform, mesh, maxwellSettings);
+
+    // run
+    maxwell.Run();
+  }
+
+  // close down MPI
+  Comm::Finalize();
+  return LIBP_SUCCESS;
+}
diff --git a/solvers/maxwell/okl/maxwellErrorTet3D.okl b/solvers/maxwell/okl/maxwellErrorTet3D.okl
new file mode 100644
index 000000000..119696532
--- /dev/null
+++ b/solvers/maxwell/okl/maxwellErrorTet3D.okl
@@ -0,0 +1,151 @@
+/*
+
+The MIT License (MIT)
+
+Copyright (c) 2017-2022 Tim Warburton, Noel Chalmers, Jesse Chan, Ali Karakus
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+
+*/
+
+
+// isotropic maxwell
+@kernel void maxwellErrorTet3D(const dlong Nelements,
+			       const dfloat time,
+				@restrict const  dfloat *  vgeo,
+				@restrict const  dfloat *  x,
+				@restrict const  dfloat *  y,
+				@restrict const  dfloat *  z,
+				@restrict const  dfloat *  MM,
+				@restrict const  dfloat *  q,
+				@restrict dfloat * normerrq){
+
+  for(dlong e=0;e<Nelements;++e;@outer(0)){
+
+    @shared dfloat s_errHx[p_Np];
+    @shared dfloat s_errHy[p_Np];
+    @shared dfloat s_errHz[p_Np];
+    @shared dfloat s_errEx[p_Np];
+    @shared dfloat s_errEy[p_Np];
+    @shared dfloat s_errEz[p_Np];
+    @shared dfloat s_errL2[p_Np];
+
+    for(int n=0;n<p_Np;++n;@inner(0)){
+      dlong id = e*p_Np + n;
+      dfloat xn = x[id], yn = y[id], zn = z[id];
+
+      dlong qbase = e*p_Np*p_Nfields + n;
+      const dfloat Hx = q[qbase+0*p_Np];
+      const dfloat Hy = q[qbase+1*p_Np];
+      const dfloat Hz = q[qbase+2*p_Np];
+      const dfloat Ex = q[qbase+3*p_Np];
+      const dfloat Ey = q[qbase+4*p_Np];
+      const dfloat Ez = q[qbase+5*p_Np];
+      
+      dfloat exactHx = 0, exactHy = 0, exactHz = 0;
+      dfloat exactEx = 0, exactEy = 0, exactEz = 0;
+      
+      maxwellInitialConditions3D(time, xn, yn, zn, &exactHx, &exactHy, &exactHz, &exactEx, &exactEy, &exactEz);
+
+      // compute nodal error
+      s_errHx[n] = Hx - exactHx;
+      s_errHy[n] = Hy - exactHy;
+      s_errHz[n] = Hz - exactHz;
+      s_errEx[n] = Ex - exactEx;
+      s_errEy[n] = Ey - exactEy;
+      s_errEz[n] = Ez - exactEz;      
+    }
+
+    for(int n=0;n<p_Np;++n;@inner(0)){
+      
+      const dfloat J = vgeo[e*p_Nvgeo + p_JID];
+
+      dfloat MerrHx = 0, MerrHy = 0, MerrHz = 0;
+      dfloat MerrEx = 0, MerrEy = 0, MerrEz = 0;
+
+      for(int m=0;m<p_Np;++m){
+	dfloat MMnm = MM[n+m*p_Np];
+	MerrHx += MMnm*s_errHx[m];
+	MerrHy += MMnm*s_errHy[m];
+	MerrHz += MMnm*s_errHz[m];
+	MerrEx += MMnm*s_errEx[m];
+	MerrEy += MMnm*s_errEy[m];
+	MerrEz += MMnm*s_errEz[m];
+      }
+
+      dfloat err = 0;
+      err += s_errHx[n]*MerrHx;
+      err += s_errHy[n]*MerrHy;
+      err += s_errHz[n]*MerrHz;
+      err += s_errEx[n]*MerrEx;
+      err += s_errEy[n]*MerrEy;
+      err += s_errEz[n]*MerrEz;
+
+      err *= J;
+
+      s_errL2[n] = err;
+    }
+
+#if p_Np>512
+    for(int n=0;n<p_Np;++n;@inner(0)) s_errL2[n] += (n+512<p_Np) ? s_errL2[n+512]:0.;
+#endif
+
+#if p_Np>256
+    for(int n=0;n<p_Np;++n;@inner(0)) s_errL2[n] += (n+256<p_Np) ? s_errL2[n+256]:0.;
+#endif
+
+#if p_Np>128
+    for(int n=0;n<p_Np;++n;@inner(0)) s_errL2[n] += (n+128<p_Np) ? s_errL2[n+128]:0.;
+#endif
+
+#if p_Np>64
+    for(int n=0;n<p_Np;++n;@inner(0)) s_errL2[n] += (n+64<p_Np) ? s_errL2[n+64]:0.;
+#endif
+
+#if p_Np>32
+    for(int n=0;n<p_Np;++n;@inner(0)) s_errL2[n] += (n+32<p_Np) ? s_errL2[n+32]:0.;
+#endif
+
+#if p_Np>16
+    for(int n=0;n<p_Np;++n;@inner(0)) s_errL2[n] += (n+16<p_Np) ? s_errL2[n+16]:0.;
+#endif
+
+#if p_Np>8
+    for(int n=0;n<p_Np;++n;@inner(0)) s_errL2[n] += (n+8<p_Np) ? s_errL2[n+8]:0.;
+#endif
+
+#if p_Np>4
+    for(int n=0;n<p_Np;++n;@inner(0)) s_errL2[n] += (n+4<p_Np) ? s_errL2[n+4]:0.;
+#endif
+
+#if p_Np>2
+    for(int n=0;n<p_Np;++n;@inner(0)) s_errL2[n] += (n+2<p_Np) ? s_errL2[n+2]:0.;
+#endif
+
+#if p_Np>1
+    for(int n=0;n<p_Np;++n;@inner(0)) s_errL2[n] += (n+1<p_Np) ? s_errL2[n+1]:0.;
+#endif
+
+    for(int n=0;n<p_Np;++n;@inner(0)){
+      if(n==0){
+	normerrq[e] = s_errL2[n];
+      }
+    }
+  }
+}
diff --git a/solvers/maxwell/okl/maxwellErrorTri2D.okl b/solvers/maxwell/okl/maxwellErrorTri2D.okl
new file mode 100644
index 000000000..810ef7f2a
--- /dev/null
+++ b/solvers/maxwell/okl/maxwellErrorTri2D.okl
@@ -0,0 +1,134 @@
+/*
+
+The MIT License (MIT)
+
+Copyright (c) 2017-2022 Tim Warburton, Noel Chalmers, Jesse Chan, Ali Karakus
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+
+*/
+
+
+// isotropic maxwell
+@kernel void maxwellErrorTri2D(const dlong Nelements,
+			       const dfloat time,
+				@restrict const  dfloat *  vgeo,
+				@restrict const  dfloat *  x,
+				@restrict const  dfloat *  y,
+				@restrict const  dfloat *  z,
+				@restrict const  dfloat *  MM,
+				@restrict const  dfloat *  q,
+				@restrict dfloat * normerrq){
+
+  for(dlong e=0;e<Nelements;++e;@outer(0)){
+
+    @shared dfloat s_errHx[p_Np];
+    @shared dfloat s_errHy[p_Np];
+    @shared dfloat s_errEz[p_Np];
+    @shared dfloat s_errL2[p_Np];
+
+
+    for(int n=0;n<p_Np;++n;@inner(0)){
+      dlong id = e*p_Np + n;
+      dfloat xn = x[id], yn = y[id];
+
+      dlong qbase = e*p_Np*p_Nfields + n;
+      const dfloat Hx = q[qbase+0*p_Np];
+      const dfloat Hy = q[qbase+1*p_Np];
+      const dfloat Ez = q[qbase+2*p_Np];
+      
+      dfloat exactHx = 0, exactHy = 0, exactEz = 0;
+      
+      maxwellInitialConditions2D(time, xn, yn, &exactHx, &exactHy, &exactEz);
+
+      // compute nodal error
+      s_errHx[n] = Hx - exactHx;
+      s_errHy[n] = Hy - exactHy;
+      s_errEz[n] = Ez - exactEz;      
+    }
+
+    for(int n=0;n<p_Np;++n;@inner(0)){
+      
+      const dfloat J = vgeo[e*p_Nvgeo + p_JID];
+
+      dfloat MerrHx = 0, MerrHy = 0, MerrEz = 0;
+
+      for(int m=0;m<p_Np;++m){
+	dfloat MMnm = MM[n+m*p_Np];
+	MerrHx += MMnm*s_errHx[m];
+	MerrHy += MMnm*s_errHy[m];
+	MerrEz += MMnm*s_errEz[m];
+      }
+
+      dfloat err = s_errHx[n]*MerrHx;
+      err += s_errHy[n]*MerrHy;
+      err += s_errEz[n]*MerrEz;
+
+      err *= J;
+
+      s_errL2[n] = err;
+    }
+
+#if p_Np>512
+    for(int n=0;n<p_Np;++n;@inner(0)) s_errL2[n] += (n+512<p_Np) ? s_errL2[n+512]:0.;
+#endif
+
+#if p_Np>256
+    for(int n=0;n<p_Np;++n;@inner(0)) s_errL2[n] += (n+256<p_Np) ? s_errL2[n+256]:0.;
+#endif
+
+#if p_Np>128
+    for(int n=0;n<p_Np;++n;@inner(0)) s_errL2[n] += (n+128<p_Np) ? s_errL2[n+128]:0.;
+#endif
+
+#if p_Np>64
+    for(int n=0;n<p_Np;++n;@inner(0)) s_errL2[n] += (n+64<p_Np) ? s_errL2[n+64]:0.;
+#endif
+
+#if p_Np>32
+    for(int n=0;n<p_Np;++n;@inner(0)) s_errL2[n] += (n+32<p_Np) ? s_errL2[n+32]:0.;
+#endif
+
+#if p_Np>16
+    for(int n=0;n<p_Np;++n;@inner(0)) s_errL2[n] += (n+16<p_Np) ? s_errL2[n+16]:0.;
+#endif
+
+#if p_Np>8
+    for(int n=0;n<p_Np;++n;@inner(0)) s_errL2[n] += (n+8<p_Np) ? s_errL2[n+8]:0.;
+#endif
+
+#if p_Np>4
+    for(int n=0;n<p_Np;++n;@inner(0)) s_errL2[n] += (n+4<p_Np) ? s_errL2[n+4]:0.;
+#endif
+
+#if p_Np>2
+    for(int n=0;n<p_Np;++n;@inner(0)) s_errL2[n] += (n+2<p_Np) ? s_errL2[n+2]:0.;
+#endif
+
+#if p_Np>1
+    for(int n=0;n<p_Np;++n;@inner(0)) s_errL2[n] += (n+1<p_Np) ? s_errL2[n+1]:0.;
+#endif
+
+    for(int n=0;n<p_Np;++n;@inner(0)){
+      if(n==0){
+	normerrq[e] = s_errL2[n];
+      }
+    }
+  }
+}
diff --git a/solvers/maxwell/okl/maxwellInitialCondition2D.okl b/solvers/maxwell/okl/maxwellInitialCondition2D.okl
new file mode 100644
index 000000000..dd6af5510
--- /dev/null
+++ b/solvers/maxwell/okl/maxwellInitialCondition2D.okl
@@ -0,0 +1,50 @@
+/*
+
+The MIT License (MIT)
+
+Copyright (c) 2017-2022 Tim Warburton, Noel Chalmers, Jesse Chan, Ali Karakus
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+
+*/
+
+@kernel void maxwellInitialCondition2D(const dlong Nelements,
+                                         const dfloat time,
+                                         @restrict const  dfloat *  x,
+                                         @restrict const  dfloat *  y,
+                                         @restrict const  dfloat *  z,
+                                         @restrict        dfloat *  q){
+
+  for(dlong e=0;e<Nelements;++e;@outer(0)){
+    for(int n=0;n<p_Np;++n;@inner(0)){
+      const dlong id = e*p_Np + n;
+
+      dfloat Hx = 0.0;
+      dfloat Hy = 0.0;
+      dfloat Ez = 0.0;
+
+      maxwellInitialConditions2D(time, x[id], y[id], &Hx, &Hy, &Ez);
+
+      const dlong qbase = e*p_Np*p_Nfields + n;
+      q[qbase+0*p_Np] = Hx;
+      q[qbase+1*p_Np] = Hy;
+      q[qbase+2*p_Np] = Ez;
+    }
+  }
+}
diff --git a/solvers/maxwell/okl/maxwellInitialCondition3D.okl b/solvers/maxwell/okl/maxwellInitialCondition3D.okl
new file mode 100644
index 000000000..08d03e283
--- /dev/null
+++ b/solvers/maxwell/okl/maxwellInitialCondition3D.okl
@@ -0,0 +1,53 @@
+/*
+
+The MIT License (MIT)
+
+Copyright (c) 2017-2022 Tim Warburton, Noel Chalmers, Jesse Chan, Ali Karakus
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+
+*/
+
+@kernel void maxwellInitialCondition3D(const dlong Nelements,
+                                         const dfloat time,
+                                         @restrict const  dfloat *  x,
+                                         @restrict const  dfloat *  y,
+                                         @restrict const  dfloat *  z,
+                                         @restrict        dfloat *  q){
+
+  for(dlong e=0;e<Nelements;++e;@outer(0)){
+    for(int n=0;n<p_Np;++n;@inner(0)){
+      const dlong id = e*p_Np + n;
+
+      dfloat Hx = 0.0, Hy = 0.0, Hz = 0.0;
+      dfloat Ex = 0.0, Ey = 0.0, Ez = 0.0;
+
+      maxwellInitialConditions3D(time, x[id], y[id], z[id], &Hx, &Hy, &Hz, &Ex, &Ey, &Ez);
+
+      const dlong qbase = e*p_Np*p_Nfields + n;
+      q[qbase+0*p_Np] = Hx;
+      q[qbase+1*p_Np] = Hy;
+      q[qbase+2*p_Np] = Hz;
+      q[qbase+3*p_Np] = Ex;
+      q[qbase+4*p_Np] = Ey;
+      q[qbase+5*p_Np] = Ez;
+
+    }
+  }
+}
diff --git a/solvers/maxwell/okl/maxwellSurfaceTet3D.okl b/solvers/maxwell/okl/maxwellSurfaceTet3D.okl
new file mode 100644
index 000000000..b9a5dd07f
--- /dev/null
+++ b/solvers/maxwell/okl/maxwellSurfaceTet3D.okl
@@ -0,0 +1,230 @@
+/*
+
+The MIT License (MIT)
+
+Copyright (c) 2017-2022 Tim Warburton, Noel Chalmers, Jesse Chan, Ali Karakus
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+
+*/
+
+
+// Roe averaged Riemann solver
+void upwind(const dfloat nx,
+            const dfloat ny,
+            const dfloat nz,
+            const dfloat HxM,
+            const dfloat HyM,
+            const dfloat HzM,
+            const dfloat ExM,
+            const dfloat EyM,
+            const dfloat EzM,
+            const dfloat HxP,
+            const dfloat HyP,
+            const dfloat HzP,
+            const dfloat ExP,
+            const dfloat EyP,
+            const dfloat EzP,
+	    dfloat *fluxHx,
+            dfloat *fluxHy,
+            dfloat *fluxHz,
+	    dfloat *fluxEx,
+            dfloat *fluxEy,
+            dfloat *fluxEz){
+
+  /*
+    ORIGINAL H&W nodal DG book code
+
+    % form field differences at faces                                                                                                                                                     
+    dHx(:)  = Hx(vmapP)-Hx(vmapM);  dEx(:)  = Ex(vmapP)-Ex(vmapM);
+    dHy(:)  = Hy(vmapP)-Hy(vmapM);  dEy(:)  = Ey(vmapP)-Ey(vmapM);
+    dHz(:)  = Hz(vmapP)-Hz(vmapM);  dEz(:)  = Ez(vmapP)-Ez(vmapM);
+    
+    % make boundary conditions all reflective (Ez+ = -Ez-)                                                                                                                                
+    dHx(mapB) = 0;  dEx(mapB) = -2*Ex(vmapB);
+    dHy(mapB) = 0;  dEy(mapB) = -2*Ey(vmapB);
+    dHz(mapB) = 0;  dEz(mapB) = -2*Ez(vmapB);
+    
+    alpha=1; % => full upwinding                                                                                                                                                          
+    
+    ndotdH = nx.*dHx + ny.*dHy + nz.*dHz;
+    ndotdE = nx.*dEx + ny.*dEy + nz.*dEz;
+
+    fluxHx = -ny.*dEz + nz.*dEy + alpha*(dHx - ndotdH.*nx);
+    fluxHy = -nz.*dEx + nx.*dEz + alpha*(dHy - ndotdH.*ny);
+    fluxHz = -nx.*dEy + ny.*dEx + alpha*(dHz - ndotdH.*nz);
+    
+    fluxEx =  ny.*dHz - nz.*dHy + alpha*(dEx - ndotdE.*nx);
+    fluxEy =  nz.*dHx - nx.*dHz + alpha*(dEy - ndotdE.*ny);
+    fluxEz =  nx.*dHy - ny.*dHx + alpha*(dEz - ndotdE.*nz);
+  */
+
+  //  form field differences at faces
+  dfloat dHx = HxP - HxM, dEx = ExP - ExM;
+  dfloat dHy = HyP - HyM, dEy = EyP - EyM;
+  dfloat dHz = HzP - HzM, dEz = EzP - EzM;
+    
+  dfloat alpha=1; // => full upwinding
+  
+  dfloat ndotdH = nx*dHx + ny*dHy + nz*dHz;
+  dfloat ndotdE = nx*dEx + ny*dEy + nz*dEz;
+
+  *fluxHx = -ny*dEz + nz*dEy + alpha*(dHx - ndotdH*nx);
+  *fluxHy = -nz*dEx + nx*dEz + alpha*(dHy - ndotdH*ny);
+  *fluxHz = -nx*dEy + ny*dEx + alpha*(dHz - ndotdH*nz);
+    
+  *fluxEx =  ny*dHz - nz*dHy + alpha*(dEx - ndotdE*nx);
+  *fluxEy =  nz*dHx - nx*dHz + alpha*(dEy - ndotdE*ny);
+  *fluxEz =  nx*dHy - ny*dHx + alpha*(dEz - ndotdE*nz);
+}
+
+// batch process elements
+@kernel void maxwellSurfaceTet3D(const dlong Nelements,
+                                  @restrict const  dlong  *  elementIds,
+                                  @restrict const  dfloat *  sgeo,
+                                  @restrict const  dfloat *  LIFT,
+                                  @restrict const  dlong  *  vmapM,
+                                  @restrict const  dlong  *  vmapP,
+                                  @restrict const  int    *  EToB,
+                                  const dfloat time,
+                                  @restrict const  dfloat *  x,
+                                  @restrict const  dfloat *  y,
+                                  @restrict const  dfloat *  z,
+                                  @restrict const  dfloat *  q,
+                                  @restrict dfloat *  rhsq){
+
+  // for all elements
+  for(dlong eo=0;eo<Nelements;eo+=p_NblockS;@outer(0)){
+
+    // @shared storage for flux terms
+    @shared dfloat s_fluxHx[p_NblockS][p_NfacesNfp];
+    @shared dfloat s_fluxHy[p_NblockS][p_NfacesNfp];
+    @shared dfloat s_fluxHz[p_NblockS][p_NfacesNfp];
+    @shared dfloat s_fluxEx[p_NblockS][p_NfacesNfp];
+    @shared dfloat s_fluxEy[p_NblockS][p_NfacesNfp];
+    @shared dfloat s_fluxEz[p_NblockS][p_NfacesNfp];
+
+    @exclusive dlong r_e, element;
+
+    // for all face nodes of all elements
+    for(int es=0;es<p_NblockS;++es;@inner(1)){
+      for(int n=0;n<p_maxNodes;++n;@inner(0)){ // maxNodes = max(Nfp*Nfaces,Np)
+        r_e = eo + es;
+        if(r_e<Nelements){
+          element = elementIds[r_e];
+
+          if(n<p_NfacesNfp){
+            // find face that owns this node
+            const int face = n/p_Nfp;
+
+            // load surface geofactors for this face
+            const dlong sid    = p_Nsgeo*(element*p_Nfaces+face);
+            const dfloat nx   = sgeo[sid+p_NXID];
+            const dfloat ny   = sgeo[sid+p_NYID];
+            const dfloat nz   = sgeo[sid+p_NZID];
+            const dfloat sJ   = sgeo[sid+p_SJID];
+            const dfloat invJ = sgeo[sid+p_IJID];
+
+            // indices of negative and positive traces of face node
+            const dlong id  = element*p_Nfp*p_Nfaces + n;
+            const dlong idM = vmapM[id];
+            const dlong idP = vmapP[id];
+
+            // load traces
+            const dlong eM = element;
+            const dlong eP = idP/p_Np;
+            const int vidM = idM%p_Np;
+            const int vidP = idP%p_Np;
+
+            const dlong qbaseM = eM*p_Np*p_Nfields + vidM;
+            const dlong qbaseP = eP*p_Np*p_Nfields + vidP;
+
+            const dfloat HxM = q[qbaseM + 0*p_Np];
+            const dfloat HyM = q[qbaseM + 1*p_Np];
+            const dfloat HzM = q[qbaseM + 2*p_Np];
+            const dfloat ExM = q[qbaseM + 3*p_Np];
+            const dfloat EyM = q[qbaseM + 4*p_Np];
+            const dfloat EzM = q[qbaseM + 5*p_Np];
+
+            dfloat HxP = q[qbaseP + 0*p_Np];
+            dfloat HyP = q[qbaseP + 1*p_Np];
+            dfloat HzP = q[qbaseP + 2*p_Np];
+            dfloat ExP = q[qbaseP + 3*p_Np];
+            dfloat EyP = q[qbaseP + 4*p_Np];
+            dfloat EzP = q[qbaseP + 5*p_Np];
+
+            // apply boundary condition
+            const int bc = EToB[face+p_Nfaces*element];
+            if(bc>0){
+              maxwellDirichletConditions3D(bc, time, x[idM], y[idM], z[idM], nx, ny, nz, HxM, HyM, HzM, ExM, EyM, EzM, &HxP, &HyP, &HzP, &ExP, &EyP, &EzP);
+            }
+
+            // evaluate "flux" terms: (sJ/J)*(A*nx+B*ny)*(q^* - q^-)
+            const dfloat sc = p_half*invJ*sJ;
+	    
+            dfloat fluxHx, fluxHy, fluxHz;
+	    dfloat fluxEx, fluxEy, fluxEz;
+
+            upwind(nx, ny, nz, HxM, HyM, HzM, ExM, EyM, EzM, HxP, HyP, HzP, ExP, EyP, EzP, &fluxHx, &fluxHy, &fluxHz, &fluxEx, &fluxEy, &fluxEz);
+
+            s_fluxHx[es][n] = sc*(fluxHx);
+            s_fluxHy[es][n] = sc*(fluxHy);
+            s_fluxHz[es][n] = sc*(fluxHz);
+            s_fluxEx[es][n] = sc*(fluxEx);
+            s_fluxEy[es][n] = sc*(fluxEy);
+            s_fluxEz[es][n] = sc*(fluxEz);
+          }
+        }
+      }
+    }
+
+    // for each node in the element
+    for(int es=0;es<p_NblockS;++es;@inner(1)){
+      for(int n=0;n<p_maxNodes;++n;@inner(0)){
+        if(r_e<Nelements){
+          if(n<p_Np){
+            // load rhs data from volume fluxes
+            dfloat LfluxHx = 0.f, LfluxHy = 0.f, LfluxHz = 0.f;
+	    dfloat LfluxEx = 0.f, LfluxEy = 0.f, LfluxEz = 0.f;
+
+            // rhs += LIFT*((sJ/J)*(A*nx+B*ny)*(q^* - q^-))
+            #pragma unroll p_NfacesNfp
+              for(int m=0;m<p_NfacesNfp;++m){
+                const dfloat L = LIFT[n+m*p_Np];
+                LfluxHx += L*s_fluxHx[es][m];
+                LfluxHy += L*s_fluxHy[es][m];
+                LfluxHz += L*s_fluxHz[es][m];
+                LfluxEx += L*s_fluxEx[es][m];
+                LfluxEy += L*s_fluxEy[es][m];
+                LfluxEz += L*s_fluxEz[es][m];
+              }
+
+            const dlong base = element*p_Np*p_Nfields+n;
+            rhsq[base+0*p_Np] += LfluxHx;
+            rhsq[base+1*p_Np] += LfluxHy;
+            rhsq[base+2*p_Np] += LfluxHz;
+            rhsq[base+3*p_Np] += LfluxEx;
+            rhsq[base+4*p_Np] += LfluxEy;
+            rhsq[base+5*p_Np] += LfluxEz;
+          }
+        }
+      }
+    }
+  }
+}
diff --git a/solvers/maxwell/okl/maxwellSurfaceTri2D.okl b/solvers/maxwell/okl/maxwellSurfaceTri2D.okl
new file mode 100644
index 000000000..efabdfc44
--- /dev/null
+++ b/solvers/maxwell/okl/maxwellSurfaceTri2D.okl
@@ -0,0 +1,195 @@
+/*
+
+The MIT License (MIT)
+
+Copyright (c) 2017-2022 Tim Warburton, Noel Chalmers, Jesse Chan, Ali Karakus
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+
+*/
+
+
+// TM mode:
+// dHx/dt = -dEz/dy           
+// dHy/dt =  dEz/dx
+// dEz/dt =  dHy/dx - dHx/dt
+
+
+// Riemann solver
+void upwind(const dfloat nx,
+            const dfloat ny,
+            const dfloat HxM,
+            const dfloat HyM,
+	    const dfloat EzM,
+	    const dfloat HxP,
+            const dfloat HyP,
+            const dfloat EzP,
+            dfloat *fluxHx,
+            dfloat *fluxHy,
+            dfloat *fluxEz){
+
+  /* 
+     ORIGINAL H&W nodal DG book codes (MaxwellRHS2D.m)
+     % Define field differences at faces                                                                                                                                                   
+     dHx = zeros(Nfp*Nfaces,K); dHx(:) = Hx(vmapM)-Hx(vmapP);
+     dHy = zeros(Nfp*Nfaces,K); dHy(:) = Hy(vmapM)-Hy(vmapP);
+     dEz = zeros(Nfp*Nfaces,K); dEz(:) = Ez(vmapM)-Ez(vmapP);
+     
+     % Impose reflective boundary conditions (Ez+ = -Ez-)                                                                                                                                  
+     dHx(mapB) = 0; dHy(mapB) = 0; dEz(mapB) = 2*Ez(vmapB);
+     
+     % evaluate upwind fluxes                                                                                                                                                              
+     alpha = 1.0;
+     ndotdH =  nx.*dHx+ny.*dHy;
+     fluxHx =  ny.*dEz + alpha*(ndotdH.*nx-dHx);
+     fluxHy = -nx.*dEz + alpha*(ndotdH.*ny-dHy);
+     fluxEz = -nx.*dHy + ny.*dHx - alpha*dEz;
+
+  */
+
+  // boundary traces (HxP, HyP, EzP) already modified for BCs
+  dfloat dHx = HxM-HxP;
+  dfloat dHy = HyM-HyP;
+  dfloat dEz = EzM-EzP;
+  
+  // evaluate upwind fluxes                                                                                                   
+  dfloat alpha = 1.0;
+  dfloat ndotdH =  nx*dHx+ny*dHy;
+  
+  *fluxHx = p_half*( ny*dEz + alpha*(ndotdH*nx-dHx));
+  *fluxHy = p_half*(-nx*dEz + alpha*(ndotdH*ny-dHy));
+  *fluxEz = p_half*(-nx*dHy + ny*dHx - alpha*dEz);
+}
+
+// batch process elements
+@kernel void maxwellSurfaceTri2D(const dlong Nelements,
+                                  @restrict const  dlong  *  elementIds,
+                                  @restrict const  dfloat *  sgeo,
+                                  @restrict const  dfloat *  LIFT,
+                                  @restrict const  dlong  *  vmapM,
+                                  @restrict const  dlong  *  vmapP,
+                                  @restrict const  int    *  EToB,
+                                  const dfloat time,
+                                  @restrict const  dfloat *  x,
+                                  @restrict const  dfloat *  y,
+                                  @restrict const  dfloat *  z,
+                                  @restrict const  dfloat *  q,
+                                  @restrict dfloat *  rhsq){
+
+  // for all elements
+  for(dlong eo=0;eo<Nelements;eo+=p_NblockS;@outer(0)){
+
+    // @shared storage for flux terms
+    @shared dfloat s_fluxHx[p_NblockS][p_NfacesNfp];
+    @shared dfloat s_fluxHy[p_NblockS][p_NfacesNfp];
+    @shared dfloat s_fluxEz[p_NblockS][p_NfacesNfp];
+
+    @exclusive dlong r_e, element;
+
+    // for all face nodes of all elements
+    for(int es=0;es<p_NblockS;++es;@inner(1)){
+      for(int n=0;n<p_maxNodes;++n;@inner(0)){ // maxNodes = max(Nfp*Nfaces,Np)
+        r_e = eo + es;
+        if(r_e<Nelements){
+          element = elementIds[r_e];
+
+          if(n<p_NfacesNfp){
+            // find face that owns this node
+            const int face = n/p_Nfp;
+
+            // load surface geofactors for this face
+            const dlong sid   = p_Nsgeo*(element*p_Nfaces+face);
+            const dfloat nx   = sgeo[sid+p_NXID];
+            const dfloat ny   = sgeo[sid+p_NYID];
+            const dfloat sJ   = sgeo[sid+p_SJID];
+            const dfloat invJ = sgeo[sid+p_IJID];
+
+            // indices of negative and positive traces of face node
+            const dlong id  = element*p_Nfp*p_Nfaces + n;
+            const dlong idM = vmapM[id];
+            const dlong idP = vmapP[id];
+
+            // load traces
+            const dlong eM = element;
+            const dlong eP = idP/p_Np;
+            const int vidM = idM%p_Np;
+            const int vidP = idP%p_Np;
+
+            const dlong qbaseM = eM*p_Np*p_Nfields + vidM;
+            const dlong qbaseP = eP*p_Np*p_Nfields + vidP;
+
+            const dfloat HxM = q[qbaseM + 0*p_Np];
+            const dfloat HyM = q[qbaseM + 1*p_Np];
+            const dfloat EzM = q[qbaseM + 2*p_Np];
+
+            dfloat HxP = q[qbaseP + 0*p_Np];
+            dfloat HyP = q[qbaseP + 1*p_Np];
+            dfloat EzP = q[qbaseP + 2*p_Np];
+
+            // apply boundary condition
+            const int bc = EToB[face+p_Nfaces*element];
+            if(bc>0){
+              maxwellDirichletConditions2D(bc, time, x[idM], y[idM], nx, ny, HxM, HyM, EzM, &HxP, &HyP, &EzP);
+              //should also add the Neumann BC here, but need uxM, uyM, vxM, abd vyM somehow
+            }
+	    
+            // evaluate "flux" terms: (sJ/J)*(A*nx+B*ny)*(q^* - q^-)
+            const dfloat sc = invJ*sJ;
+
+            dfloat fluxHx, fluxHy, fluxEz;
+	    
+            upwind(nx, ny, HxM, HyM, EzM, HxP, HyP, EzP, &fluxHx, &fluxHy, &fluxEz);
+
+            // const dfloat hinv = sgeo[sid + p_IHID];
+            // dfloat penalty = p_Nq*p_Nq*hinv*mu;
+
+            s_fluxHx[es][n] = sc*(fluxHx);
+            s_fluxHy[es][n] = sc*(fluxHy);
+            s_fluxEz[es][n] = sc*(fluxEz);
+          }
+        }
+      }
+    }
+
+    // for each node in the element
+    for(int es=0;es<p_NblockS;++es;@inner(1)){
+      for(int n=0;n<p_maxNodes;++n;@inner(0)){
+        if(r_e<Nelements){
+          if(n<p_Np){
+            // load rhs data from volume fluxes
+            dfloat LfluxHx = 0.f, LfluxHy = 0.f, LfluxEz = 0.f;
+
+#pragma unroll p_NfacesNfp
+	    for(int m=0;m<p_NfacesNfp;++m){
+	      const dfloat L = LIFT[n+m*p_Np];
+	      LfluxHx += L*s_fluxHx[es][m];
+	      LfluxHy += L*s_fluxHy[es][m];
+	      LfluxEz += L*s_fluxEz[es][m];
+	    }
+
+            const dlong base = element*p_Np*p_Nfields+n;
+            rhsq[base+0*p_Np] += LfluxHx;
+            rhsq[base+1*p_Np] += LfluxHy;
+            rhsq[base+2*p_Np] += LfluxEz;
+          }
+        }
+      }
+    }
+  }
+}
diff --git a/solvers/maxwell/okl/maxwellVolumeTet3D.okl b/solvers/maxwell/okl/maxwellVolumeTet3D.okl
new file mode 100644
index 000000000..cbc5ed24c
--- /dev/null
+++ b/solvers/maxwell/okl/maxwellVolumeTet3D.okl
@@ -0,0 +1,192 @@
+/*
+
+  The MIT License (MIT)
+
+  Copyright (c) 2017-2022 Tim Warburton, Noel Chalmers, Jesse Chan, Ali Karakus
+
+  Permission is hereby granted, free of charge, to any person obtaining a copy
+  of this software and associated documentation files (the "Software"), to deal
+  in the Software without restriction, including without limitation the rights
+  to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+  copies of the Software, and to permit persons to whom the Software is
+  furnished to do so, subject to the following conditions:
+
+  The above copyright notice and this permission notice shall be included in all
+  copies or substantial portions of the Software.
+
+  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+  IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+  AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+  OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+  SOFTWARE.
+
+*/
+
+
+// isotropic maxwell
+
+// thread loop over elements
+@kernel void maxwellVolumeTet3D(const dlong Nelements,
+				@restrict const  dfloat *  vgeo,
+				@restrict const  dfloat *  D,
+				@restrict const  dfloat *  q,
+				@restrict dfloat *  rhsq){
+
+#define p_Nvol 1
+
+#if p_N==1
+#define p_NblockV 32
+#elif p_N==2
+#define p_NblockV 24
+#elif p_N==3
+#define p_NblockV 12
+#elif p_N==4
+#define p_NblockV 7
+#elif p_N==5
+#define p_NblockV 4
+#elif p_N==6
+#define p_NblockV 3
+#elif p_N==7
+#define p_NblockV 2
+#else
+#define p_NblockV 1
+#endif
+
+
+  for(dlong eo=0;eo<Nelements;eo+=(p_Nvol*p_NblockV);@outer(0)){
+
+    @shared dfloat s_Hx[p_Nvol][p_NblockV][p_Np];
+    @shared dfloat s_Hy[p_Nvol][p_NblockV][p_Np];
+    @shared dfloat s_Hz[p_Nvol][p_NblockV][p_Np];
+    @shared dfloat s_Ex[p_Nvol][p_NblockV][p_Np];
+    @shared dfloat s_Ey[p_Nvol][p_NblockV][p_Np];
+    @shared dfloat s_Ez[p_Nvol][p_NblockV][p_Np];
+
+    for(int et=0;et<p_NblockV;++et;@inner(1)){
+      for(int n=0;n<p_Np;++n;@inner(0)){
+
+#pragma unroll p_Nvol
+	for(int es=0;es<p_Nvol;++es){
+
+	  const dlong e = es*p_NblockV + et + eo;
+
+	  if(e<Nelements){
+
+	    const dlong  qbase = e*p_Np*p_Nfields + n;
+	    s_Hx[es][et][n] = q[qbase+0*p_Np];
+	    s_Hy[es][et][n] = q[qbase+1*p_Np];
+	    s_Hz[es][et][n] = q[qbase+2*p_Np];
+	    s_Ex[es][et][n] = q[qbase+3*p_Np];
+	    s_Ey[es][et][n] = q[qbase+4*p_Np];
+	    s_Ez[es][et][n] = q[qbase+5*p_Np];
+	  }
+	}
+      }
+    }
+
+    for(int et=0;et<p_NblockV;++et;@inner(1)){
+      for(int n=0;n<p_Np;++n;@inner(0)){
+
+        dfloat r_dHxdr[p_Nvol] = {0.f}, r_dHxds[p_Nvol] = {0.f}, r_dHxdt[p_Nvol] = {0.f};
+        dfloat r_dHydr[p_Nvol] = {0.f}, r_dHyds[p_Nvol] = {0.f}, r_dHydt[p_Nvol] = {0.f};
+        dfloat r_dHzdr[p_Nvol] = {0.f}, r_dHzds[p_Nvol] = {0.f}, r_dHzdt[p_Nvol] = {0.f};
+
+        dfloat r_dExdr[p_Nvol] = {0.f}, r_dExds[p_Nvol] = {0.f}, r_dExdt[p_Nvol] = {0.f};
+        dfloat r_dEydr[p_Nvol] = {0.f}, r_dEyds[p_Nvol] = {0.f}, r_dEydt[p_Nvol] = {0.f};
+        dfloat r_dEzdr[p_Nvol] = {0.f}, r_dEzds[p_Nvol] = {0.f}, r_dEzdt[p_Nvol] = {0.f};
+
+	
+#pragma unroll p_Np
+	for(int m=0;m<p_Np;++m){
+
+	  const dfloat Drnm = D[n+m*p_Np];
+	  const dfloat Dsnm = D[n+m*p_Np+1*p_Np*p_Np];
+	  const dfloat Dtnm = D[n+m*p_Np+2*p_Np*p_Np];
+
+#pragma unroll p_Nvol
+	  for(int es=0;es<p_Nvol;++es){
+	    const dfloat Hxm = s_Hx[es][et][m];
+	    const dfloat Hym = s_Hy[es][et][m];
+	    const dfloat Hzm = s_Hz[es][et][m];
+	    const dfloat Exm = s_Ex[es][et][m];
+	    const dfloat Eym = s_Ey[es][et][m];
+	    const dfloat Ezm = s_Ez[es][et][m];
+
+	    r_dHxdr[es] += Drnm*Hxm;
+	    r_dHxds[es] += Dsnm*Hxm;
+	    r_dHxdt[es] += Dtnm*Hxm;
+	    r_dHydr[es] += Drnm*Hym;
+	    r_dHyds[es] += Dsnm*Hym;
+	    r_dHydt[es] += Dtnm*Hym;
+	    r_dHzdr[es] += Drnm*Hzm;
+	    r_dHzds[es] += Dsnm*Hzm;
+	    r_dHzdt[es] += Dtnm*Hzm;
+
+	    r_dExdr[es] += Drnm*Exm;
+	    r_dExds[es] += Dsnm*Exm;
+	    r_dExdt[es] += Dtnm*Exm;
+	    r_dEydr[es] += Drnm*Eym;
+	    r_dEyds[es] += Dsnm*Eym;
+	    r_dEydt[es] += Dtnm*Eym;
+	    r_dEzdr[es] += Drnm*Ezm;
+	    r_dEzds[es] += Dsnm*Ezm;
+	    r_dEzdt[es] += Dtnm*Ezm;
+	  }
+	}
+
+#pragma unroll p_Nvol
+	for(int es=0;es<p_Nvol;++es){
+
+	  const dlong e = es*p_NblockV + et + eo;
+
+	  if(e<Nelements){
+	    // prefetch geometric factors (constant on triangle)
+	    const dfloat drdx = vgeo[e*p_Nvgeo + p_RXID];
+	    const dfloat drdy = vgeo[e*p_Nvgeo + p_RYID];
+	    const dfloat drdz = vgeo[e*p_Nvgeo + p_RZID];
+	    const dfloat dsdx = vgeo[e*p_Nvgeo + p_SXID];
+	    const dfloat dsdy = vgeo[e*p_Nvgeo + p_SYID];
+	    const dfloat dsdz = vgeo[e*p_Nvgeo + p_SZID];
+	    const dfloat dtdx = vgeo[e*p_Nvgeo + p_TXID];
+	    const dfloat dtdy = vgeo[e*p_Nvgeo + p_TYID];
+	    const dfloat dtdz = vgeo[e*p_Nvgeo + p_TZID];
+
+	    const dlong base = e*p_Np*p_Nfields + n;
+
+	    const dfloat dHxdy = drdy*r_dHxdr[es] + dsdy*r_dHxds[es] + dtdy*r_dHxdt[es];
+	    const dfloat dHxdz = drdz*r_dHxdr[es] + dsdz*r_dHxds[es] + dtdz*r_dHxdt[es];
+
+	    const dfloat dHydx = drdx*r_dHydr[es] + dsdx*r_dHyds[es] + dtdx*r_dHydt[es];
+	    const dfloat dHydz = drdz*r_dHydr[es] + dsdz*r_dHyds[es] + dtdz*r_dHydt[es];
+
+	    const dfloat dHzdx = drdx*r_dHzdr[es] + dsdx*r_dHzds[es] + dtdx*r_dHzdt[es];
+	    const dfloat dHzdy = drdy*r_dHzdr[es] + dsdy*r_dHzds[es] + dtdy*r_dHzdt[es];
+
+	    const dfloat dExdy = drdy*r_dExdr[es] + dsdy*r_dExds[es] + dtdy*r_dExdt[es];
+	    const dfloat dExdz = drdz*r_dExdr[es] + dsdz*r_dExds[es] + dtdz*r_dExdt[es];
+
+	    const dfloat dEydx = drdx*r_dEydr[es] + dsdx*r_dEyds[es] + dtdx*r_dEydt[es];
+	    const dfloat dEydz = drdz*r_dEydr[es] + dsdz*r_dEyds[es] + dtdz*r_dEydt[es];
+
+	    const dfloat dEzdx = drdx*r_dEzdr[es] + dsdx*r_dEzds[es] + dtdx*r_dEzdt[es];
+	    const dfloat dEzdy = drdy*r_dEzdr[es] + dsdy*r_dEzds[es] + dtdy*r_dEzdt[es];
+
+	    // dH/dt = -curl E
+	    rhsq[base+0*p_Np] = -(dEzdy-dEydz);
+	    rhsq[base+1*p_Np] = -(dExdz-dEzdx);
+	    rhsq[base+2*p_Np] = -(dEydx-dExdy);
+
+	    // dE/dt =  curl H
+	    rhsq[base+3*p_Np] =  (dHzdy-dHydz);
+	    rhsq[base+4*p_Np] =  (dHxdz-dHzdx);
+	    rhsq[base+5*p_Np] =  (dHydx-dHxdy);
+	  }
+	}
+      }
+    }
+  }
+}
+
+
diff --git a/solvers/maxwell/okl/maxwellVolumeTri2D.okl b/solvers/maxwell/okl/maxwellVolumeTri2D.okl
new file mode 100644
index 000000000..f18f75ae4
--- /dev/null
+++ b/solvers/maxwell/okl/maxwellVolumeTri2D.okl
@@ -0,0 +1,107 @@
+/*
+
+The MIT License (MIT)
+
+Copyright (c) 2017-2022 Tim Warburton, Noel Chalmers, Jesse Chan, Ali Karakus
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+
+*/
+
+
+// isotropic maxwell
+@kernel void maxwellVolumeTri2D(const dlong Nelements,
+				@restrict const  dfloat *  vgeo,
+				@restrict const  dfloat *  D,
+				@restrict const  dfloat *  q,
+				@restrict dfloat *  rhsq){
+
+  for(dlong e=0;e<Nelements;++e;@outer(0)){
+
+    @shared dfloat s_F[p_Nfields][p_Np];
+    @shared dfloat s_G[p_Nfields][p_Np];
+
+    for(int n=0;n<p_Np;++n;@inner(0)){
+
+      // prefetch geometric factors (constant on triangle)
+      const dfloat drdx = vgeo[e*p_Nvgeo + p_RXID];
+      const dfloat drdy = vgeo[e*p_Nvgeo + p_RYID];
+      const dfloat dsdx = vgeo[e*p_Nvgeo + p_SXID];
+      const dfloat dsdy = vgeo[e*p_Nvgeo + p_SYID];
+
+      const dlong  qbase = e*p_Np*p_Nfields + n;
+      const dfloat Hx = q[qbase+0*p_Np];
+      const dfloat Hy = q[qbase+1*p_Np];
+      const dfloat Ez = q[qbase+2*p_Np];
+
+      // dH/dt = -curl E
+      // dE/dt =  curl H
+      
+      {
+	// dHx/dt = -dEz/dy
+        const dfloat f = 0;
+        const dfloat g = -Ez;
+        s_F[0][n] = drdx*f + drdy*g;
+        s_G[0][n] = dsdx*f + dsdy*g;
+      }
+
+      {
+	// dHy/dt =  dEz/dx
+        const dfloat f = Ez;
+        const dfloat g = 0;
+        s_F[1][n] = drdx*f + drdy*g;
+        s_G[1][n] = dsdx*f + dsdy*g;
+      }
+
+      {
+	// dEz/dt = dHy/dx - dHx/dy
+        const dfloat f =  Hy;
+        const dfloat g = -Hx;
+        s_F[2][n] = drdx*f + drdy*g;
+        s_G[2][n] = dsdx*f + dsdy*g;
+      }
+    }
+
+    for(int n=0;n<p_Np;++n;@inner(0)){
+
+      dfloat rhsq0 = 0, rhsq1 = 0, rhsq2 = 0;
+
+      for(int i=0;i<p_Np;++i){
+        const dfloat Drni = D[n+i*p_Np+0*p_Np*p_Np];
+        const dfloat Dsni = D[n+i*p_Np+1*p_Np*p_Np];
+
+        rhsq0 += Drni*s_F[0][i]
+                +Dsni*s_G[0][i];
+        rhsq1 += Drni*s_F[1][i]
+                +Dsni*s_G[1][i];
+        rhsq2 += Drni*s_F[2][i]
+                +Dsni*s_G[2][i];
+      }
+
+      const dlong base = e*p_Np*p_Nfields + n;
+
+      // move to rhs
+      rhsq[base+0*p_Np] = rhsq0;
+      rhsq[base+1*p_Np] = rhsq1;
+      rhsq[base+2*p_Np] = rhsq2;
+    }
+  }
+}
+
+
diff --git a/solvers/maxwell/setups/setupTet3D.rc b/solvers/maxwell/setups/setupTet3D.rc
new file mode 100644
index 000000000..40184f092
--- /dev/null
+++ b/solvers/maxwell/setups/setupTet3D.rc
@@ -0,0 +1,70 @@
+[FORMAT]
+2.0
+
+[DATA FILE]
+#data/maxwellGaussian3D.h
+data/maxwellMode3D.h
+
+[MESH FILE]
+BOX
+
+[MESH DIMENSION]
+3
+
+[ELEMENT TYPE] # number of edges
+6
+
+[BOX NX]
+6
+
+[BOX NY]
+6
+
+[BOX NZ]
+6
+
+[BOX DIMX]
+2
+
+[BOX DIMY]
+2
+
+[BOX DIMZ]
+2
+
+
+[BOX BOUNDARY FLAG]
+1
+
+[POLYNOMIAL DEGREE]
+6
+
+[THREAD MODEL]
+CUDA
+
+[PLATFORM NUMBER]
+0
+
+[DEVICE NUMBER]
+0
+
+[TIME INTEGRATOR]
+DOPRI5
+
+[START TIME]
+0
+
+[FINAL TIME]
+10
+
+[OUTPUT INTERVAL]
+1
+
+[OUTPUT TO FILE]
+FALSE
+
+[OUTPUT FILE NAME]
+maxwell
+
+[COMPUTE ERROR]
+TRUE
\ No newline at end of file
diff --git a/solvers/maxwell/setups/setupTri2D.rc b/solvers/maxwell/setups/setupTri2D.rc
new file mode 100644
index 000000000..8f725e360
--- /dev/null
+++ b/solvers/maxwell/setups/setupTri2D.rc
@@ -0,0 +1,70 @@
+[FORMAT]
+2.0
+
+[DATA FILE]
+data/maxwellMode2D.h
+#data/maxwellGaussian2D.h
+
+[MESH FILE]
+BOX
+
+[MESH DIMENSION]
+2
+
+[ELEMENT TYPE] # number of edges
+3
+
+[BOX NX]
+10
+
+[BOX NY]
+10
+
+[BOX NZ]
+10
+
+[BOX DIMX]
+2
+
+[BOX DIMY]
+2
+
+[BOX DIMZ]
+2
+
+[BOX BOUNDARY FLAG]
+1
+
+[POLYNOMIAL DEGREE]
+7
+
+[THREAD MODEL]
+CUDA
+
+[PLATFORM NUMBER]
+0
+
+[DEVICE NUMBER]
+0
+
+[TIME INTEGRATOR]
+#LSERK4
+DOPRI5
+
+[START TIME]
+0
+
+[FINAL TIME]
+10
+
+[OUTPUT INTERVAL]
+1
+
+[OUTPUT TO FILE]
+TRUE
+
+[OUTPUT FILE NAME]
+maxwell
+
+[COMPUTE ERROR]
+TRUE
\ No newline at end of file
diff --git a/solvers/maxwell/src/maxwellPlotFields.cpp b/solvers/maxwell/src/maxwellPlotFields.cpp
new file mode 100644
index 000000000..231a77987
--- /dev/null
+++ b/solvers/maxwell/src/maxwellPlotFields.cpp
@@ -0,0 +1,174 @@
+/*
+
+The MIT License (MIT)
+
+Copyright (c) 2017-2022 Tim Warburton, Noel Chalmers, Jesse Chan, Ali Karakus
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+
+*/
+
+#include "maxwell.hpp"
+
+// interpolate data to plot nodes and save to file (one per process
+void maxwell_t::PlotFields(memory<dfloat> Q, const std::string fileName){
+
+  FILE *fp;
+
+  fp = fopen(fileName.c_str(), "w");
+
+  fprintf(fp, "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"BigEndian\">\n");
+  fprintf(fp, "  <UnstructuredGrid>\n");
+  fprintf(fp, "    <Piece NumberOfPoints=\"%d\" NumberOfCells=\"%d\">\n",
+          mesh.Nelements*mesh.plotNp,
+          mesh.Nelements*mesh.plotNelements);
+
+  // write out nodes
+  fprintf(fp, "      <Points>\n");
+  fprintf(fp, "        <DataArray type=\"Float32\" NumberOfComponents=\"3\" Format=\"ascii\">\n");
+
+  //scratch space for interpolation
+  size_t Nscratch = std::max(mesh.Np, mesh.plotNp);
+  memory<dfloat> scratch(2*Nscratch);
+
+  memory<dfloat> Ix(mesh.plotNp);
+  memory<dfloat> Iy(mesh.plotNp);
+  memory<dfloat> Iz(mesh.plotNp);
+
+  // compute plot node coordinates on the fly
+  for(dlong e=0;e<mesh.Nelements;++e){
+    mesh.PlotInterp(mesh.x + e*mesh.Np, Ix, scratch);
+    mesh.PlotInterp(mesh.y + e*mesh.Np, Iy, scratch);
+    if(mesh.dim==3)
+      mesh.PlotInterp(mesh.z + e*mesh.Np, Iz, scratch);
+
+    if (mesh.dim==2) {
+      for(int n=0;n<mesh.plotNp;++n){
+        fprintf(fp, "       ");
+        fprintf(fp, "%g %g %g\n", Ix[n],Iy[n],0.0);
+      }
+    } else {
+      for(int n=0;n<mesh.plotNp;++n){
+        fprintf(fp, "       ");
+        fprintf(fp, "%g %g %g\n", Ix[n],Iy[n],Iz[n]);
+      }
+    }
+  }
+  fprintf(fp, "        </DataArray>\n");
+  fprintf(fp, "      </Points>\n");
+
+  memory<dfloat> IHx(mesh.plotNp);
+  memory<dfloat> IHy(mesh.plotNp);
+  memory<dfloat> IHz(mesh.plotNp);
+
+  // write out magnetic field
+  fprintf(fp, "      <PointData Scalars=\"scalars\">\n");
+  fprintf(fp, "        <DataArray type=\"Float32\" Name=\"MagneticField\" NumberOfComponents=\"%d\" Format=\"ascii\">\n", mesh.dim);
+  for(dlong e=0;e<mesh.Nelements;++e){
+    mesh.PlotInterp(Q + 0*mesh.Np + e*mesh.Np*Nfields, IHx, scratch);
+    mesh.PlotInterp(Q + 1*mesh.Np + e*mesh.Np*Nfields, IHy, scratch);
+    if(mesh.dim==3)
+      mesh.PlotInterp(Q + 2*mesh.Np + e*mesh.Np*Nfields, IHz, scratch);
+
+    for(int n=0;n<mesh.plotNp;++n){
+      fprintf(fp, "       ");
+      fprintf(fp, "       ");
+      if (mesh.dim==2)
+        fprintf(fp, "%f %f\n", IHx[n], IHy[n]);
+      else
+        fprintf(fp, "%f %f %f\n", IHx[n], IHy[n], IHz[n]);
+    }
+  }
+  fprintf(fp, "       </DataArray>\n");
+
+  memory<dfloat> IEx(mesh.plotNp);
+  memory<dfloat> IEy(mesh.plotNp);
+  memory<dfloat> IEz(mesh.plotNp);
+  
+  // write out electric field
+  if(mesh.dim==2){
+    fprintf(fp, "        <DataArray type=\"Float32\" Name=\"ElectricField\" NumberOfComponents=\"%d\" Format=\"ascii\">\n", 1);
+  }else{
+    fprintf(fp, "        <DataArray type=\"Float32\" Name=\"ElectricField\" NumberOfComponents=\"%d\" Format=\"ascii\">\n", mesh.dim);
+  }
+
+  for(dlong e=0;e<mesh.Nelements;++e){
+    if(mesh.dim==2){
+      mesh.PlotInterp(Q + 2*mesh.Np + e*mesh.Np*Nfields, IEz, scratch);
+    }
+    else{
+      mesh.PlotInterp(Q + 3*mesh.Np + e*mesh.Np*Nfields, IEx, scratch);
+      mesh.PlotInterp(Q + 4*mesh.Np + e*mesh.Np*Nfields, IEy, scratch);
+      mesh.PlotInterp(Q + 5*mesh.Np + e*mesh.Np*Nfields, IEz, scratch);
+    }
+
+    for(int n=0;n<mesh.plotNp;++n){
+      fprintf(fp, "       ");
+      fprintf(fp, "       ");
+      if (mesh.dim==2)
+        fprintf(fp, "%f\n", IEz[n]);
+      else
+        fprintf(fp, "%f %f %f\n", IEx[n], IEy[n], IEz[n]);
+    }
+  }
+  fprintf(fp, "       </DataArray>\n");
+  fprintf(fp, "     </PointData>\n");
+
+  
+  fprintf(fp, "    <Cells>\n");
+  fprintf(fp, "      <DataArray type=\"Int32\" Name=\"connectivity\" Format=\"ascii\">\n");
+
+  for(dlong e=0;e<mesh.Nelements;++e){
+    for(int n=0;n<mesh.plotNelements;++n){
+      fprintf(fp, "       ");
+      for(int m=0;m<mesh.plotNverts;++m){
+        fprintf(fp, "%d ", e*mesh.plotNp + mesh.plotEToV[n*mesh.plotNverts+m]);
+      }
+      fprintf(fp, "\n");
+    }
+  }
+  fprintf(fp, "        </DataArray>\n");
+
+  fprintf(fp, "        <DataArray type=\"Int32\" Name=\"offsets\" Format=\"ascii\">\n");
+  dlong cnt = 0;
+  for(dlong e=0;e<mesh.Nelements;++e){
+    for(int n=0;n<mesh.plotNelements;++n){
+      cnt += mesh.plotNverts;
+      fprintf(fp, "       ");
+      fprintf(fp, "%d\n", cnt);
+    }
+  }
+  fprintf(fp, "       </DataArray>\n");
+
+  fprintf(fp, "       <DataArray type=\"Int32\" Name=\"types\" Format=\"ascii\">\n");
+  for(dlong e=0;e<mesh.Nelements;++e){
+    for(int n=0;n<mesh.plotNelements;++n){
+      if(mesh.dim==2)
+        fprintf(fp, "5\n");
+      else
+        fprintf(fp, "10\n");
+    }
+  }
+  fprintf(fp, "        </DataArray>\n");
+  fprintf(fp, "      </Cells>\n");
+  fprintf(fp, "    </Piece>\n");
+  fprintf(fp, "  </UnstructuredGrid>\n");
+  fprintf(fp, "</VTKFile>\n");
+  fclose(fp);
+}
diff --git a/solvers/maxwell/src/maxwellReport.cpp b/solvers/maxwell/src/maxwellReport.cpp
new file mode 100644
index 000000000..13d2c21a5
--- /dev/null
+++ b/solvers/maxwell/src/maxwellReport.cpp
@@ -0,0 +1,69 @@
+/*
+
+The MIT License (MIT)
+
+Copyright (c) 2017-2022 Tim Warburton, Noel Chalmers, Jesse Chan, Ali Karakus
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+
+*/
+
+#include "maxwell.hpp"
+
+void maxwell_t::Report(dfloat time, int tstep){
+
+  static int frame=0;
+  dfloat L2error = 0;
+
+  if (settings.compareSetting("COMPUTE ERROR","TRUE")) {
+    deviceMemory<dfloat> o_err = platform.reserve<dfloat>(mesh.Nelements);
+    // only set up for Tri and Tet elements
+    errorKernel(mesh.Nelements, time, mesh.o_vgeo, mesh.o_x, mesh.o_y, mesh.o_z, mesh.o_MM, o_q, o_err);
+    L2error = sqrt(platform.linAlg().sum(mesh.Nelements, o_err, mesh.comm));
+  }
+
+  //compute q.M*q
+  dlong Nentries = mesh.Nelements*mesh.Np*Nfields;
+  deviceMemory<dfloat> o_Mq = platform.reserve<dfloat>(Nentries);
+  mesh.MassMatrixApply(o_q, o_Mq);
+
+  dfloat norm2 = sqrt(platform.linAlg().innerProd(Nentries, o_q, o_Mq, mesh.comm));
+
+  if(mesh.rank==0) {
+    if (settings.compareSetting("COMPUTE ERROR","TRUE")) {
+      printf("%5.2f (%d), %5.2f, %5.2e (time, timestep, norm, L2error)\n", time, tstep, norm2, L2error);
+    } else {
+      printf("%5.2f (%d), %5.2f (time, timestep, norm)\n", time, tstep, norm2);
+    }
+  }
+
+  if (settings.compareSetting("OUTPUT TO FILE","TRUE")) {
+
+    // copy data back to host
+    o_q.copyTo(q);
+
+    // output field files
+    std::string name;
+    settings.getSetting("OUTPUT FILE NAME", name);
+    char fname[BUFSIZ];
+    sprintf(fname, "%s_%04d_%04d.vtu", name.c_str(), mesh.rank, frame++);
+
+    PlotFields(q, std::string(fname));
+  }
+}
diff --git a/solvers/maxwell/src/maxwellRun.cpp b/solvers/maxwell/src/maxwellRun.cpp
new file mode 100644
index 000000000..e37ae6263
--- /dev/null
+++ b/solvers/maxwell/src/maxwellRun.cpp
@@ -0,0 +1,66 @@
+/*
+
+The MIT License (MIT)
+
+Copyright (c) 2017-2022 Tim Warburton, Noel Chalmers, Jesse Chan, Ali Karakus
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+
+*/
+
+#include "maxwell.hpp"
+
+void maxwell_t::Run(){
+
+  dfloat startTime, finalTime;
+  settings.getSetting("START TIME", startTime);
+  settings.getSetting("FINAL TIME", finalTime);
+
+  initialConditionKernel(mesh.Nelements,
+                         startTime,
+                         mesh.o_x,
+                         mesh.o_y,
+                         mesh.o_z,
+                         o_q);
+
+  dfloat cfl=1.0;
+  settings.getSetting("CFL NUMBER", cfl);
+
+  // set time step
+  dfloat hmin = mesh.MinCharacteristicLength();
+  dfloat vmax = MaxWaveSpeed();
+
+  dfloat dt = cfl*hmin/(vmax*(mesh.N+1.)*(mesh.N+1.));
+  timeStepper.SetTimeStep(dt);
+
+  timeStepper.Run(*this, o_q, startTime, finalTime);
+
+  // output norm of final solution
+  {
+    //compute q.M*q
+    dlong Nentries = mesh.Nelements*mesh.Np*Nfields;
+    deviceMemory<dfloat> o_Mq = platform.reserve<dfloat>(Nentries);
+    mesh.MassMatrixApply(o_q, o_Mq);
+
+    dfloat norm2 = sqrt(platform.linAlg().innerProd(Nentries, o_q, o_Mq, mesh.comm));
+
+    if(mesh.rank==0)
+      printf("Solution norm = %17.15lg\n", norm2);
+  }
+}
diff --git a/solvers/maxwell/src/maxwellSettings.cpp b/solvers/maxwell/src/maxwellSettings.cpp
new file mode 100644
index 000000000..d9938894b
--- /dev/null
+++ b/solvers/maxwell/src/maxwellSettings.cpp
@@ -0,0 +1,109 @@
+/*
+
+The MIT License (MIT)
+
+Copyright (c) 2017-2022 Tim Warburton, Noel Chalmers, Jesse Chan, Ali Karakus
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+
+*/
+
+#include "maxwell.hpp"
+
+//settings for maxwell solver
+maxwellSettings_t::maxwellSettings_t(comm_t _comm):
+  settings_t(_comm) {
+
+  newSetting("DATA FILE",
+             "data/maxwellGaussian2D.h",
+             "Boundary and Initial conditions header");
+
+  newSetting("TIME INTEGRATOR",
+             "DOPRI5",
+             "Time integration method",
+             {"AB3", "DOPRI5", "LSERK4"});
+
+  newSetting("CFL NUMBER",
+             "1.0",
+             "Multiplier for timestep stability bound");
+
+  newSetting("START TIME",
+             "0",
+             "Start time for time integration");
+
+  newSetting("FINAL TIME",
+             "10",
+             "End time for time integration");
+
+  newSetting("OUTPUT INTERVAL",
+             ".1",
+             "Time between printing output data");
+
+  newSetting("OUTPUT TO FILE",
+             "FALSE",
+             "Flag for writing fields to VTU files",
+             {"TRUE", "FALSE"});
+
+  newSetting("OUTPUT FILE NAME",
+             "maxwell");
+
+  newSetting("COMPUTE ERROR",
+             "FALSE",
+             "Flag for computing L2 error (uses InitialCondition from data file)",
+             {"TRUE", "FALSE"});
+
+}
+
+void maxwellSettings_t::report() {
+
+  if (comm.rank()==0) {
+    std::cout << "Maxwell Settings:\n\n";
+    reportSetting("DATA FILE");
+    reportSetting("TIME INTEGRATOR");
+    reportSetting("START TIME");
+    reportSetting("FINAL TIME");
+    reportSetting("OUTPUT INTERVAL");
+    reportSetting("OUTPUT TO FILE");
+    reportSetting("OUTPUT FILE NAME");
+    reportSetting("COMPUTE ERROR");
+  }
+}
+
+void maxwellSettings_t::parseFromFile(platformSettings_t& platformSettings,
+                                  meshSettings_t& meshSettings,
+                                  const std::string filename) {
+  //read all settings from file
+  settings_t s(comm);
+  s.readSettingsFromFile(filename);
+
+  for(auto it = s.settings.begin(); it != s.settings.end(); ++it) {
+    setting_t& set = it->second;
+    const std::string name = set.getName();
+    const std::string val = set.getVal<std::string>();
+    if (platformSettings.hasSetting(name))
+      platformSettings.changeSetting(name, val);
+    else if (meshSettings.hasSetting(name))
+      meshSettings.changeSetting(name, val);
+    else if (hasSetting(name)) //self
+      changeSetting(name, val);
+    else  {
+      LIBP_FORCE_ABORT("Unknown setting: [" << name << "] requested");
+    }
+  }
+}
diff --git a/solvers/maxwell/src/maxwellSetup.cpp b/solvers/maxwell/src/maxwellSetup.cpp
new file mode 100644
index 000000000..d830d874c
--- /dev/null
+++ b/solvers/maxwell/src/maxwellSetup.cpp
@@ -0,0 +1,139 @@
+/*
+
+The MIT License (MIT)
+
+Copyright (c) 2017-2022 Tim Warburton, Noel Chalmers, Jesse Chan, Ali Karakus
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+
+*/
+
+#include "maxwell.hpp"
+
+void maxwell_t::Setup(platform_t& _platform, mesh_t& _mesh,
+                        maxwellSettings_t& _settings){
+
+  platform = _platform;
+  mesh = _mesh;
+  comm = _mesh.comm;
+  settings = _settings;
+
+  Nfields = (mesh.dim==3) ? 6:3;
+
+  dlong Nlocal = mesh.Nelements*mesh.Np*Nfields;
+  dlong Nhalo  = mesh.totalHaloPairs*mesh.Np*Nfields;
+
+  //Trigger JIT kernel builds
+  ogs::InitializeKernels(platform, ogs::Dfloat, ogs::Add);
+
+  //setup linear algebra module
+  platform.linAlg().InitKernels({"innerProd", "sum", "norm2"});
+
+  /*setup trace halo exchange */
+  traceHalo = mesh.HaloTraceSetup(Nfields);
+
+  //setup timeStepper
+  if (settings.compareSetting("TIME INTEGRATOR","AB3")){
+    timeStepper.Setup<TimeStepper::ab3>(mesh.Nelements,
+                                        mesh.totalHaloPairs,
+                                        mesh.Np, Nfields, platform, comm);
+  } else if (settings.compareSetting("TIME INTEGRATOR","LSERK4")){
+    timeStepper.Setup<TimeStepper::lserk4>(mesh.Nelements,
+                                           mesh.totalHaloPairs,
+                                           mesh.Np, Nfields, platform, comm);
+  } else if (settings.compareSetting("TIME INTEGRATOR","DOPRI5")){
+    timeStepper.Setup<TimeStepper::dopri5>(mesh.Nelements,
+                                           mesh.totalHaloPairs,
+                                           mesh.Np, Nfields, platform, comm);
+  }
+
+  // set penalty parameter
+  dfloat Lambda2 = 0.5;
+
+  // compute samples of q at interpolation nodes
+  q.malloc(Nlocal+Nhalo);
+  o_q = platform.malloc<dfloat>(Nlocal+Nhalo);
+
+  mesh.MassMatrixKernelSetup(Nfields); // mass matrix operator
+
+  // OCCA build stuff
+  properties_t kernelInfo = mesh.props; //copy base occa properties
+
+  //add boundary data to kernel info
+  std::string dataFileName;
+  settings.getSetting("DATA FILE", dataFileName);
+  kernelInfo["includes"] += dataFileName;
+
+  kernelInfo["defines/" "p_Nfields"]= Nfields;
+
+  const dfloat p_half = 1./2.;
+  kernelInfo["defines/" "p_half"]= p_half;
+
+  int maxNodes = std::max(mesh.Np, (mesh.Nfp*mesh.Nfaces));
+  kernelInfo["defines/" "p_maxNodes"]= maxNodes;
+
+  int blockMax = 256;
+  if (platform.device.mode() == "CUDA") blockMax = 512;
+
+  int NblockV = std::max(1, blockMax/mesh.Np);
+  kernelInfo["defines/" "p_NblockV"]= NblockV;
+
+  int NblockS = std::max(1, blockMax/maxNodes);
+  kernelInfo["defines/" "p_NblockS"]= NblockS;
+
+  kernelInfo["defines/" "p_Lambda2"]= Lambda2;
+
+  // set kernel name suffix
+  std::string suffix = mesh.elementSuffix();
+  std::string oklFilePrefix = DMAXWELL "/okl/";
+  std::string oklFileSuffix = ".okl";
+
+  std::string fileName, kernelName;
+
+  // kernels from volume file
+  fileName   = oklFilePrefix + "maxwellVolume" + suffix + oklFileSuffix;
+  kernelName = "maxwellVolume" + suffix;
+
+  volumeKernel =  platform.buildKernel(fileName, kernelName,
+                                         kernelInfo);
+  // kernels from surface file
+  fileName   = oklFilePrefix + "maxwellSurface" + suffix + oklFileSuffix;
+  kernelName = "maxwellSurface" + suffix;
+
+  surfaceKernel = platform.buildKernel(fileName, kernelName,
+                                         kernelInfo);
+
+  // kernels from surface file
+  fileName   = oklFilePrefix + "maxwellError" + suffix + oklFileSuffix;
+  kernelName = "maxwellError" + suffix;
+
+  errorKernel = platform.buildKernel(fileName, kernelName,
+				     kernelInfo);
+  
+  if (mesh.dim==2) {
+    fileName   = oklFilePrefix + "maxwellInitialCondition2D" + oklFileSuffix;
+    kernelName = "maxwellInitialCondition2D";
+  } else {
+    fileName   = oklFilePrefix + "maxwellInitialCondition3D" + oklFileSuffix;
+    kernelName = "maxwellInitialCondition3D";
+  }
+
+  initialConditionKernel = platform.buildKernel(fileName, kernelName,
+                                                  kernelInfo);
+}
diff --git a/solvers/maxwell/src/maxwellStep.cpp b/solvers/maxwell/src/maxwellStep.cpp
new file mode 100644
index 000000000..5c7262797
--- /dev/null
+++ b/solvers/maxwell/src/maxwellStep.cpp
@@ -0,0 +1,78 @@
+/*
+
+The MIT License (MIT)
+
+Copyright (c) 2017-2022 Tim Warburton, Noel Chalmers, Jesse Chan, Ali Karakus
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+
+*/
+
+#include "maxwell.hpp"
+
+dfloat maxwell_t::MaxWaveSpeed(){
+  //wavespeed is constant 1 everywhere
+  const dfloat vmax = 1.0;
+  return vmax;
+}
+
+//evaluate ODE rhs = f(q,t)
+void maxwell_t::rhsf(deviceMemory<dfloat>& o_Q, deviceMemory<dfloat>& o_RHS, const dfloat T){
+
+  // extract q halo on DEVICE
+  traceHalo.ExchangeStart(o_Q, 1);
+
+  volumeKernel(mesh.Nelements,
+               mesh.o_vgeo,
+               mesh.o_D,
+               o_Q,
+               o_RHS);
+
+  if (mesh.NinternalElements)
+    surfaceKernel(mesh.NinternalElements,
+                  mesh.o_internalElementIds,
+                  mesh.o_sgeo,
+                  mesh.o_LIFT,
+                  mesh.o_vmapM,
+                  mesh.o_vmapP,
+                  mesh.o_EToB,
+                  T,
+                  mesh.o_x,
+                  mesh.o_y,
+                  mesh.o_z,
+                  o_Q,
+                  o_RHS);
+
+  traceHalo.ExchangeFinish(o_Q, 1);
+
+  if (mesh.NhaloElements)
+    surfaceKernel(mesh.NhaloElements,
+                  mesh.o_haloElementIds,
+                  mesh.o_sgeo,
+                  mesh.o_LIFT,
+                  mesh.o_vmapM,
+                  mesh.o_vmapP,
+                  mesh.o_EToB,
+                  T,
+                  mesh.o_x,
+                  mesh.o_y,
+                  mesh.o_z,
+                  o_Q,
+                  o_RHS);
+}
diff --git a/test/makefile b/test/makefile
index 2e2b3da30..6b3dc5d6d 100644
--- a/test/makefile
+++ b/test/makefile
@@ -45,7 +45,7 @@ Can use "make verbose=true" for verbose output.
 
 endef
 
-ifeq (,$(filter info help test test-mesh test-gradient test-advection test-acoustics \
+ifeq (,$(filter info help test test-mesh test-gradient test-advection test-acoustics test-maxwell \
 				test-elliptic test-fpe test-cns test-bns test-lbs test-ins test-initial-guess test-core,$(MAKECMDGOALS)))
 ifneq (,$(MAKECMDGOALS))
 $(error ${TEST_HELP_MSG})
@@ -68,7 +68,7 @@ MESH_DIR     =${LIBP_LIBS_DIR}/mesh
 CORE_DIR     =${LIBP_DIR}/core
 TEST_DIR     =${LIBP_DIR}/test
 
-.PHONY: all help info test test-mesh test-gradient test-advection test-acoustics \
+.PHONY: all help info test test-mesh test-gradient test-advection test-acoustics test-maxwell \
 				test-elliptic test-fpe test-cns test-bns test-ins test-initial-guess test-core
 
 
@@ -119,6 +119,9 @@ test-bns:
 test-lbs:
 	@./testLbs.py
 
+test-maxwell:
+	@./testMaxwell.py
+
 test-ins:
 	@./testIns.py
 
diff --git a/test/test.py b/test/test.py
index f182435c8..0ba3df9e9 100755
--- a/test/test.py
+++ b/test/test.py
@@ -47,6 +47,7 @@
 bnsDir           = solverDir + "/bns"
 lbsDir           = solverDir + "/lbs"
 insDir           = solverDir + "/ins"
+maxwellDir       = solverDir + "/maxwell"
 
 gradientBin  = gradientDir      + "/gradientMain"
 advectionBin = advectionDir     + "/advectionMain"
@@ -57,6 +58,7 @@
 bnsBin       = bnsDir           + "/bnsMain"
 lbsBin       = lbsDir           + "/lbsMain"
 insBin       = insDir           + "/insMain"
+maxwellBin   = maxwellDir       + "/maxwellMain"
 
 inputRC = testDir + "/setup.rc"
 
@@ -164,6 +166,7 @@ def test(name, cmd, settings, referenceNorm, ranks=1):
   import testMesh
   import testGradient
   import testAdvection
+  import testMaxwell
   import testAcoustics
   import testElliptic
   import testFokkerPlanck
@@ -182,6 +185,7 @@ def test(name, cmd, settings, referenceNorm, ranks=1):
   failCount+=testParAdogs.main()
   failCount+=testGradient.main()
   failCount+=testAdvection.main()
+  failCount+=testMaxwell.main()
   failCount+=testAcoustics.main()
   failCount+=testElliptic.main()
   failCount+=testFokkerPlanck.main()
@@ -194,4 +198,5 @@ def test(name, cmd, settings, referenceNorm, ranks=1):
   failCount+=testLinearSolver.main()
   failCount+=testParAlmond.main()
 
+
   sys.exit(failCount)
diff --git a/test/testMaxwell.py b/test/testMaxwell.py
new file mode 100755
index 000000000..8668fee48
--- /dev/null
+++ b/test/testMaxwell.py
@@ -0,0 +1,87 @@
+#!/usr/bin/env python3
+
+#####################################################################################
+#
+#The MIT License (MIT)
+#
+#Copyright (c) 2017-2022 Tim Warburton, Noel Chalmers, Jesse Chan, Ali Karakus
+#
+#Permission is hereby granted, free of charge, to any person obtaining a copy
+#of this software and associated documentation files (the "Software"), to deal
+#in the Software without restriction, including without limitation the rights
+#to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+#copies of the Software, and to permit persons to whom the Software is
+#furnished to do so, subject to the following conditions:
+#
+#The above copyright notice and this permission notice shall be included in all
+#copies or substantial portions of the Software.
+#
+#THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+#IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+#FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+#AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+#LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+#OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+#SOFTWARE.
+#
+#####################################################################################
+
+from test import *
+
+data2D = maxwellDir + "/data/maxwellGaussian2D.h"
+data3D = maxwellDir + "/data/maxwellGaussian3D.h"
+
+def maxwellSettings(rcformat="2.0", data_file=data2D,
+                     mesh="BOX", dim=2, element=4, nx=10, ny=10, nz=10, boundary_flag=-1,
+                     degree=4, thread_model=device, platform_number=0, device_number=0,
+                      time_integrator="DOPRI5", cfl=1.0, start_time=0.0, final_time=1.0,
+                      output_to_file="FALSE"):
+  return [setting_t("FORMAT", rcformat),
+          setting_t("DATA FILE", data_file),
+          setting_t("MESH FILE", mesh),
+          setting_t("MESH DIMENSION", dim),
+          setting_t("ELEMENT TYPE", element),
+          setting_t("BOX GLOBAL NX", nx),
+          setting_t("BOX GLOBAL NY", ny),
+          setting_t("BOX GLOBAL NZ", nz),
+          setting_t("BOX BOUNDARY FLAG", boundary_flag),
+          setting_t("POLYNOMIAL DEGREE", degree),
+          setting_t("THREAD MODEL", thread_model),
+          setting_t("PLATFORM NUMBER", platform_number),
+          setting_t("DEVICE NUMBER", device_number),
+          setting_t("TIME INTEGRATOR", time_integrator),
+          setting_t("CFL NUMBER", cfl),
+          setting_t("START TIME", start_time),
+          setting_t("FINAL TIME", final_time),
+          setting_t("OUTPUT TO FILE", output_to_file)]
+
+def main():
+  failCount=0;
+
+  failCount += test(name="testMaxwellTri",
+                    cmd=maxwellBin,
+                    settings=maxwellSettings(element=3,data_file=data2D,dim=2),
+                    referenceNorm=0.723919093728707)
+
+  failCount += test(name="testMaxwellTet",
+                    cmd=maxwellBin,
+                    settings=maxwellSettings(element=6,data_file=data3D,dim=3,
+                                               degree=2),
+                    referenceNorm=0.494395985627954)
+
+  failCount += test(name="testMaxwellTri_MPI", ranks=4,
+                    cmd=maxwellBin,
+                    settings=maxwellSettings(element=3,data_file=data2D,dim=2,output_to_file="TRUE"),
+                    referenceNorm=0.723919093728707)
+
+  #clean up
+  for file_name in os.listdir(testDir):
+    if file_name.endswith('.vtu'):
+      os.remove(testDir + "/" + file_name)
+
+  return failCount
+
+if __name__ == "__main__":
+  failCount=0;
+  failCount+=main()
+  sys.exit(failCount)