From f54d0d7832b946ff1f61a3953b37326cc7111d21 Mon Sep 17 00:00:00 2001
From: Florian <florian.gacon@inria.fr>
Date: Fri, 18 Sep 2020 09:35:06 +0200
Subject: [PATCH] Add an ExplicitEuler Scheme
---
src/bvpy/ivp.py | 33 +++++++++++++++++++++++----------
1 file changed, 23 insertions(+), 10 deletions(-)
diff --git a/src/bvpy/ivp.py b/src/bvpy/ivp.py
index a19fba6..6eddfbf 100644
--- a/src/bvpy/ivp.py
+++ b/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)
--
GitLab