Add an ExplicitEuler Scheme

src/bvpy/ivp.py

@@ -34,10 +34,10 @@ class Time():
 
 
 class FirstOrderScheme():
-    def __init__(self, vform, dt=1):
+    def __init__(self, vform):
         self._vform = None
         self.previous_solution = None
-        self.dt = dt
+        self.dt = None
 
         self.set_vform(vform)
 
@@ -58,9 +58,6 @@ class FirstOrderScheme():
     def set_vform(self, vform):
         self._vform = vform
 
-    def set_dt(self, dt):
-        self.dt = dt
-
     @property
     def dx(self):
         return self._vform.dx
@@ -70,8 +67,8 @@ class FirstOrderScheme():
         return self._vform.ds
 
 class ImplicitEuler(FirstOrderScheme):
-    def __init__(self, vform, dt=1):
-        super().__init__(vform, dt)
+    def __init__(self, vform):
+        super().__init__(vform)
         self.type = 'implicit'
 
     def apply(self, u, v, sol):
@@ -85,6 +82,22 @@ class ImplicitEuler(FirstOrderScheme):
         # self._vform.lhs, self._vform.rhs = lhs, rhs
         return lhs, rhs
 
+class ExplicitEuler(FirstOrderScheme):
+    def __init__(self, vform):
+        super().__init__(vform)
+        self.type = 'explicit'
+
+    def apply(self, u, v, sol):
+        lhs, rhs = self._vform.get_form(self.previous_solution, v, sol)
+        F = fe.Constant(self.dt)*lhs + fe.Constant(self.dt)*rhs
+
+        F += u*v*self.dx
+        F -= self.previous_solution*v*self.dx
+
+        # self._vform.lhs, self._vform.rhs = lhs, rhs
+        return fe.lhs(F), fe.rhs(F)
+
+
 class IVP(BVP):
     def __init__(self, domain=None, vform=None, scheme=None, bc=None):
         super().__init__(domain=domain, vform=vform, bc=bc)
@@ -106,7 +119,7 @@ class IVP(BVP):
         self.time.current = time
 
     def solve(self, lhs, rhs, **kwargs):
-        if isinstance(self.vform.rhs, Zero):
+        if isinstance(rhs, Zero):
             self.solver = NonLinearSolver()
             jacobian = fe.derivative(lhs, self.solution, self.trialFunction)
             self.solver.set_problem(self.solution, lhs,
@@ -152,7 +165,7 @@ if __name__ =='__main__':
 
     # implicit_euler = ImplicitEuler(poisson, dt=0.04)
 
-    ivp = IVP(rectangle, poisson, ImplicitEuler, diri)
+    ivp = IVP(rectangle, poisson, ExplicitEuler, diri)
     sol_0 = create_expression('exp(-a*pow(x, 2) - a*pow(y, 2))', degree=2, a=5)
     ivp.set_initial_solution(sol_0)
     ivp.set_time(0)
@@ -160,5 +173,5 @@ if __name__ =='__main__':
     with fe.XDMFFile('solution.xdmf') as f:
         f.write(ivp.solution, ivp.time.current)
         for i in range(10):
-            ivp.step(0.04)
+            ivp.step(1e-4)
             f.write(ivp.solution, ivp.time.current)