From d9619d4fafbf38512642b1d8095caa05c3d18fad Mon Sep 17 00:00:00 2001
From: Eric Dalissier <Eric.Dalissier@imag.fr>
Date: Thu, 6 Dec 2012 17:33:12 +0000
Subject: [PATCH] EDO solver : validated
---
.../hysop/operator/differentialOperator_d.py | 52 +++++
HySoP/hysop/operator/stretching_d.py | 5 -
.../particular_solvers/integrator/euler.py | 29 ++-
.../integrator/integrator.py | 4 +-
.../integrator/runge_kutta2.py | 30 ++-
.../integrator/runge_kutta2stretching.py | 2 +
.../integrator/runge_kutta3.py | 69 +++++++
.../integrator/runge_kutta4.py | 35 +++-
.../integrator/runge_kutta4stretching.py | 1 +
.../test_EDO_erreur.py | 179 ++++++++++++++++++
.../test/test_particular_solvers/test_RK.py | 136 +++++++------
11 files changed, 445 insertions(+), 97 deletions(-)
create mode 100755 HySoP/hysop/particular_solvers/integrator/runge_kutta3.py
create mode 100644 HySoP/hysop/test/test_particular_solvers/test_EDO_erreur.py
diff --git a/HySoP/hysop/operator/differentialOperator_d.py b/HySoP/hysop/operator/differentialOperator_d.py
index 9493292fb..3200094aa 100755
--- a/HySoP/hysop/operator/differentialOperator_d.py
+++ b/HySoP/hysop/operator/differentialOperator_d.py
@@ -265,6 +265,58 @@ class DifferentialOperator_d(DiscreteOperator):
self.maxgersh[0] = max(max1,max2,max3)
# self.maxgersh[0] = max(maxadv1,maxadv2,maxadv3)
+ elif self.choice == 'divConservationGradU':
+
+## fourth order scheme
+# First line of Grad U matrix
+ self.result[0][...] = (1.0 * self.field2[0][np.roll(ind,+2),...] -
+ 8.0 * self.field2[0][np.roll(ind,+1),...] +\
+ 8.0 * self.field2[0][np.roll(ind,-1),...] -\
+ 1.0 * self.field2[0][np.roll(ind,-2),...]) / (12. * mesh_size[0])
+
+ self.result[1][...] = (1.0 * self.field2[0][:,np.roll(ind,+2),:] -\
+ 8.0 * self.field2[0][:,np.roll(ind,+1),:] +\
+ 8.0 * self.field2[0][:,np.roll(ind,-1),:] -\
+ 1.0 * self.field2[0][:,np.roll(ind,-2),:]) / (12. * mesh_size[1])
+
+ self.result[2][...] = (1.0 * self.field2[0][...,np.roll(ind,+2)] -\
+ 8.0 * self.field2[0][...,np.roll(ind,+1)] +\
+ 8.0 * self.field2[0][...,np.roll(ind,-1)] -\
+ 1.0 * self.field2[0][...,np.roll(ind,-2)]) / (12. * mesh_size[2])
+
+# Second line of Grad U matrix
+ self.result[3][...] = (1.0 * self.field2[1][np.roll(ind,+2),...] -
+ 8.0 * self.field2[1][np.roll(ind,+1),...] +\
+ 8.0 * self.field2[1][np.roll(ind,-1),...] -\
+ 1.0 * self.field2[1][np.roll(ind,-2),...]) / (12. * mesh_size[0])
+
+ self.result[4][...] = (1.0 * self.field2[1][:,np.roll(ind,+2),:] -\
+ 8.0 * self.field2[1][:,np.roll(ind,+1),:] +\
+ 8.0 * self.field2[1][:,np.roll(ind,-1),:] -\
+ 1.0 * self.field2[1][:,np.roll(ind,-2),:]) / (12. * mesh_size[1])
+
+ self.result[5][...] = (1.0 * self.field2[1][...,np.roll(ind,+2)] -\
+ 8.0 * self.field2[1][...,np.roll(ind,+1)] +\
+ 8.0 * self.field2[1][...,np.roll(ind,-1)] -\
+ 1.0 * self.field2[1][...,np.roll(ind,-2)]) / (12. * mesh_size[2])
+
+# Third line of Grad U matrix
+ self.result[6][...] = (1.0 * self.field2[2][np.roll(ind,+2),...] -
+ 8.0 * self.field2[2][np.roll(ind,+1),...] +\
+ 8.0 * self.field2[2][np.roll(ind,-1),...] -\
+ 1.0 * self.field2[2][np.roll(ind,-2),...]) / (12. * mesh_size[0])
+
+ self.result[7][...] = (1.0 * self.field2[2][:,np.roll(ind,+2),:] -\
+ 8.0 * self.field2[2][:,np.roll(ind,+1),:] +\
+ 8.0 * self.field2[2][:,np.roll(ind,-1),:] -\
+ 1.0 * self.field2[2][:,np.roll(ind,-2),:]) / (12. * mesh_size[1])
+
+ self.result[8][...] = (1.0 * self.field2[2][...,np.roll(ind,+2)] -\
+ 8.0 * self.field2[2][...,np.roll(ind,+1)] +\
+ 8.0 * self.field2[2][...,np.roll(ind,-1)] -\
+ 1.0 * self.field2[2][...,np.roll(ind,-2)]) / (12. * mesh_size[2])
+
+
elif self.choice == 'curl':
temp1 = (1.0 * self.field1[2][:,np.roll(ind,+2),:] -
diff --git a/HySoP/hysop/operator/stretching_d.py b/HySoP/hysop/operator/stretching_d.py
index 88ed5c5f1..0271e9b8c 100755
--- a/HySoP/hysop/operator/stretching_d.py
+++ b/HySoP/hysop/operator/stretching_d.py
@@ -72,7 +72,6 @@ class Stretching_d(DiscreteOperator):
self.stretchTerm.discreteOperator.apply()
# Time integration of dw/dt = div(wu)
-# self.method=Euler
ndt = self.stabilityTest(dt, self.maxgersh, self.method)
dtstretch = dt/float(ndt)
@@ -84,20 +83,16 @@ class Stretching_d(DiscreteOperator):
for i in range(self.vorticity.topology.dim):
self.res_vorticity[i][...] = methodInt.integrate(self.vorticity, self.stretchTermResult, t, dtstretch, direction=i)
-# self.method=RK2Stretch
if RK2Stretch == self.method :
print "RK2"
methodInt=self.method(self.stretchTerm.discreteOperator)
self.res_vorticity = methodInt.integrate(self.stretchTerm, self.stretchTermResult, self.vorticity, self.velocity, self.res_vorticity, t, dtstretch, self.vorticity.topology.dim)
-# self.method=RK3Stretch
if RK3Stretch == self.method :
print "RK3"
methodInt=self.method(self.stretchTerm.discreteOperator)
self.res_vorticity = methodInt.integrate(self.stretchTerm, self.stretchTermResult, self.vorticity, self.velocity, self.res_vorticity, t, dtstretch, self.vorticity.topology.dim)
-
-# self.method=RK4Stretch
if RK4Stretch == self.method :
print "RK4"
methodInt=self.method(self.stretchTerm.discreteOperator)
diff --git a/HySoP/hysop/particular_solvers/integrator/euler.py b/HySoP/hysop/particular_solvers/integrator/euler.py
index 8c5f8468c..7fbca8814 100755
--- a/HySoP/hysop/particular_solvers/integrator/euler.py
+++ b/HySoP/hysop/particular_solvers/integrator/euler.py
@@ -4,6 +4,7 @@
Forward Euler method.
"""
from .integrator import ODESolver
+import numpy as np
#providedClass = "Euler"
@@ -27,17 +28,31 @@ class Euler(ODESolver):
"""
ODESolver.__init__(self, f)
- def integrate(self, y, f, t, dt, direction):
+# def integrate(self, y, f, t, dt, direction):
+# """
+# Integration step for forward Euler method.
+# y(t_{n+1}) = y(t_n) + dt*f(y(t_n))
+
+# @param y : position at time t.
+# @param t : current time.
+# @param dt : time step.
+# @return y at t+dt.
+# """
+# return y[direction][...] +dt * f[direction][...]
+
+ def integrate(self, f, t , dt , u, topo ):
"""
- Integration step for forward Euler method.
- y(t_{n+1}) = y(t_n) + dt*f(y(t_n))
+ Integration step for Euler Method.
+ up = f[t, u(t,x,y,z)].
- @param y : position at time t.
+ @param u : function of position at time t.
+ @param f : function of function of position at time t.
@param t : current time.
@param dt : time step.
- @return y at t+dt.
- """
- return y[direction][...] +dt * f[direction][...]
+# """
+ result = u[:,:,:,:] + dt * np.asarray(f(t,u[:,:,:,:]))
+ return result
+
def __str__(self):
"""ToString method"""
diff --git a/HySoP/hysop/particular_solvers/integrator/integrator.py b/HySoP/hysop/particular_solvers/integrator/integrator.py
index f915d1a96..f418635a0 100644
--- a/HySoP/hysop/particular_solvers/integrator/integrator.py
+++ b/HySoP/hysop/particular_solvers/integrator/integrator.py
@@ -4,7 +4,7 @@
ODESolver interface.
"""
from ...constants import *
-
+import numpy as np
class ODESolver:
"""
@@ -23,7 +23,7 @@ class ODESolver:
if f != None:
self.f = f
else:
- self.f = lambda t, y: [0. for x in y]
+ self.f = lambda t, y: np.vectorize(f)(t,y)# [0. for x in y]
def integrate(self, y, t, dt):
"""
diff --git a/HySoP/hysop/particular_solvers/integrator/runge_kutta2.py b/HySoP/hysop/particular_solvers/integrator/runge_kutta2.py
index b9056f0c3..98982cf82 100755
--- a/HySoP/hysop/particular_solvers/integrator/runge_kutta2.py
+++ b/HySoP/hysop/particular_solvers/integrator/runge_kutta2.py
@@ -4,6 +4,7 @@
RK2 method interface.
"""
from .integrator import ODESolver
+import numpy as np
class RK2(ODESolver):
@@ -24,18 +25,35 @@ class RK2(ODESolver):
"""
ODESolver.__init__(self, f)
- def integrate(self, y, fy, yp, t, dt, dir):
+# def integrate(self, y, fy, yp, t, dt, dir):
+# """
+# Integration step for RK2 Method.
+# yp= y + dt*y'[y+dt/2*y'(y)].
+
+# @param y : position at time t.
+# @param fy : y'(y) value at time t.
+# @param yp : result position at time t + dt
+# @param t : current time.
+# @param dt : time step.
+# """
+# pass
+
+ def integrate(self, f, t , dt , u, topo ):
+
"""
Integration step for RK2 Method.
- yp= y + dt*y'[y+dt/2*y'(y)].
+ up = f[t, u(t,x,y,z)].
- @param y : position at time t.
- @param fy : y'(y) value at time t.
- @param yp : result position at time t + dt
+ @param u : function of position at time t.
+ @param f : function of function of position at time t.
@param t : current time.
@param dt : time step.
"""
- pass
+#K1 = dt *f[t,u(t,x,y,z)]
+ K1 = dt * np.asarray(f(t,u))
+#result = u(t,x,y,z)] + dt *f[t + dt/2 , u(t,x,y,z) + k1/2.]
+ result = u + dt *np.asarray(f(t+dt/2.,u + K1/2.))
+ return result
def __str__(self):
"""ToString method"""
diff --git a/HySoP/hysop/particular_solvers/integrator/runge_kutta2stretching.py b/HySoP/hysop/particular_solvers/integrator/runge_kutta2stretching.py
index 8fb0bca80..25e7b282c 100755
--- a/HySoP/hysop/particular_solvers/integrator/runge_kutta2stretching.py
+++ b/HySoP/hysop/particular_solvers/integrator/runge_kutta2stretching.py
@@ -32,6 +32,7 @@ class RK2Stretch(ODESolver):
def integrate(self, f, fd ,field1 , field2, res_vorticity, t, dt, dimension):
+# def integrate(self, f,field1 , field2, res_vorticity, t, dt, dimension):
"""
Integration step for RK4 Method
yp= y + dt*y'[y+dt/2*y'(y)].
@@ -44,6 +45,7 @@ class RK2Stretch(ODESolver):
@param t : current time.
@param dt : time step.
"""
+ print "ehoh bis"
maxgersh = np.zeros(1, dtype=PARMES_REAL, order=ORDER)
resultTmp = np.asarray([np.zeros((field2.topology.resolution), dtype=PARMES_REAL, order=ORDER) for d in xrange(field2.topology.dim)])
K1=np.asarray(fd)
diff --git a/HySoP/hysop/particular_solvers/integrator/runge_kutta3.py b/HySoP/hysop/particular_solvers/integrator/runge_kutta3.py
new file mode 100755
index 000000000..6438426ef
--- /dev/null
+++ b/HySoP/hysop/particular_solvers/integrator/runge_kutta3.py
@@ -0,0 +1,69 @@
+# -*- coding: utf-8 -*-
+"""
+@package Utils
+RK3 method interface.
+"""
+from .integrator import ODESolver
+import numpy as np
+
+
+class RK3(ODESolver):
+ """
+ ODESolver implementation for solving an equation system with RK3 method.
+ y'(t) = f(t,y)
+
+ y(t_n+1)= y(t_n) + dt*y'[y(t_n)+dt/2*y'(y(t_n))].
+
+ """
+ def __init__(self, f=None, dim=3):
+ """
+ Constructor.
+
+ @param f function f(t,y) : Right hand side of the equation to solve.
+ @param dim : dimensions
+ @param conditions : function to apply boundary conditions.
+ """
+ ODESolver.__init__(self, f)
+
+# def integrate(self, y, fy, yp, t, dt, dir):
+# """
+# Integration step for RK3 Method.
+# yp= y + dt*y'[y+dt/2*y'(y)].
+
+# @param y : position at time t.
+# @param fy : y'(y) value at time t.
+# @param yp : result position at time t + dt
+# @param t : current time.
+# @param dt : time step.
+# """
+# pass
+
+ def integrate(self, f, t , dt , u, topo ):
+
+ """
+ Integration step for RK3 Method.
+ up = f[t, u(t,x,y,z)].
+
+ @param u : function of position at time t.
+ @param f : function of function of position at time t.
+ @param t : current time.
+ @param dt : time step.
+ """
+#K1 = dt *f[t,u(t,x,y,z)]
+ K1 = dt * np.asarray(f(t,u))
+#result = u(t,x,y,z)] + dt *f[t + dt/2 , u(t,x,y,z) + k1/2.]
+ K2 = dt *np.asarray(f(t+dt/3.,u + K1[:,:,:,:]/3.))
+ K3 = dt *np.asarray(f(t+dt*2./3.,u + K2[:,:,:,:]*2./3.))
+#result = u(t,x,y,z)] + 1/4 [ k1 + 3 k3]
+ result = u + 1./4.*(K1+3.*K3)
+ return result
+
+ def __str__(self):
+ """ToString method"""
+ return "RK3 Method"
+
+
+if __name__ == "__main__":
+ print __doc__
+ print "- Provided class : RK3"
+ print RK3.__doc__
diff --git a/HySoP/hysop/particular_solvers/integrator/runge_kutta4.py b/HySoP/hysop/particular_solvers/integrator/runge_kutta4.py
index 6abd3e777..1bbedadc1 100755
--- a/HySoP/hysop/particular_solvers/integrator/runge_kutta4.py
+++ b/HySoP/hysop/particular_solvers/integrator/runge_kutta4.py
@@ -4,6 +4,7 @@
RK4 method interface.
"""
from .integrator import ODESolver
+import numpy as np
class RK4(ODESolver):
@@ -24,18 +25,38 @@ class RK4(ODESolver):
"""
ODESolver.__init__(self, f)
- def integrate(self, y, fy, yp, t, dt):
+# def integrate(self, y, fy, yp, t, dt, dir):
+# """
+# Integration step for RK4 Method.
+# yp= y + dt*y'[y+dt/2*y'(y)].
+
+# @param y : position at time t.
+# @param fy : y'(y) value at time t.
+# @param yp : result position at time t + dt
+# @param t : current time.
+# @param dt : time step.
+# """
+# pass
+
+ def integrate(self, f, t , dt , u, topo ):
+
"""
- Integration step for RK2 Method.
- yp= y + dt*y'[y+dt/2*y'(y)].
+ Integration step for RK4 Method.
+ up = f[t, u(t,x,y,z)].
- @param y : position at time t.
- @param fy : y'(y) value at time t.
- @param yp : result position at time t + dt
+ @param u : function of position at time t.
+ @param f : function of function of position at time t.
@param t : current time.
@param dt : time step.
"""
- pass
+#K1 = dt *f[t,u(t,x,y,z)]
+ K1 = dt * np.asarray(f(t,u))
+ K2 = dt *np.asarray(f(t+dt/2.,u + K1[:,:,:,:]/2.))
+ K3 = dt *np.asarray(f(t+dt/2.,u + K2[:,:,:,:]/2.))
+ K4 = dt *np.asarray(f(t+dt,u + K3[:,:,:,:]))
+#result = u(t,x,y,z)] + dt *f[t + dt/2 , u(t,x,y,z) + k1/2.]
+ result = u + 1./6. * (K1[:,:,:,:] + 2.*K2[:,:,:,:] + 2.*K3[:,:,:,:]+ K4[:,:,:,:])
+ return result
def __str__(self):
"""ToString method"""
diff --git a/HySoP/hysop/particular_solvers/integrator/runge_kutta4stretching.py b/HySoP/hysop/particular_solvers/integrator/runge_kutta4stretching.py
index 6f35ba8f2..b16238fbd 100755
--- a/HySoP/hysop/particular_solvers/integrator/runge_kutta4stretching.py
+++ b/HySoP/hysop/particular_solvers/integrator/runge_kutta4stretching.py
@@ -42,6 +42,7 @@ class RK4Stretch(ODESolver):
@param t : current time.
@param dt : time step.
"""
+ print "ehoh"
maxgersh = np.zeros(1, dtype=PARMES_REAL, order=ORDER)
resultTmp = np.asarray([np.zeros((field2.topology.resolution), dtype=PARMES_REAL, order=ORDER) for d in xrange(field2.topology.dim)])
K1=np.asarray(fd)
diff --git a/HySoP/hysop/test/test_particular_solvers/test_EDO_erreur.py b/HySoP/hysop/test/test_particular_solvers/test_EDO_erreur.py
new file mode 100644
index 000000000..64487fbc4
--- /dev/null
+++ b/HySoP/hysop/test/test_particular_solvers/test_EDO_erreur.py
@@ -0,0 +1,179 @@
+# -*- coding: utf-8 -*-
+import time
+import parmepy as pp
+from parmepy.particular_solvers.integrator.euler import Euler
+from parmepy.particular_solvers.integrator.runge_kutta2 import RK2
+from parmepy.particular_solvers.integrator.runge_kutta3 import RK3
+from parmepy.particular_solvers.integrator.runge_kutta4 import RK4
+from parmepy.particular_solvers.integrator.runge_kutta2stretching import RK2Stretch
+from parmepy.particular_solvers.integrator.runge_kutta3stretching import RK3Stretch
+from parmepy.particular_solvers.integrator.runge_kutta4stretching import RK4Stretch
+from math import *
+import unittest
+import numpy as np
+import numpy.testing as npt
+import copy
+import matplotlib.pyplot as plt
+from parmepy.constants import *
+
+
+class test_EDO(unittest.TestCase):
+ """
+ DiscreteVariable test class
+ """
+
+ def analyticalSolution(self,t, x, y, z):
+ sx = (t*np.exp(t) + 1.) * np.exp(-t)
+ sy = (t*np.exp(t) + 1.) * np.exp(-t)
+ sz = (t*np.exp(t) + 1.) * np.exp(-t)
+ return sx, sy, sz
+
+ def f (self, t, u):
+ fx = -u[0][:,:,:] + t + 1.
+ fy = -u[1][:,:,:] + t + 1.
+ fz = -u[2][:,:,:] + t + 1.
+ return fx , fy ,fz
+
+ def U (self, t, x, y, z) :
+ wx = 1.
+ wy = 1.
+ wz = 1.
+ return wx, wy, wz
+
+ def testIntegratorEDO(self):
+ # Parameters
+ nb = 32
+ timeStep = 0.1
+ finalTime = 1.0
+ multi = 2.
+ maxerror=0.
+
+ t0 = time.time()
+ self.t = 0.
+ ## Domain
+ box = pp.Box(3, length=[4.*np.pi, 4.*np.pi, 4.*np.pi],
+ origin=[- 2.*np.pi,- 2.*np.pi,- 2.*np.pi])
+
+
+#####################################################################################
+
+ ## Solver creation (discretisation of objects is done in solver initialisation)
+ topo3D = pp.CartesianTopology(domain=box, resolution=[nb, nb, nb], dim=3)
+
+ Result=np.asarray([np.zeros([nb, nb, nb],
+ dtype=PARMES_REAL, order=ORDER) for d in xrange(3)])
+ Analvorty=np.asarray([np.zeros([nb, nb, nb],
+ dtype=PARMES_REAL, order=ORDER) for d in xrange(3)])
+
+ compteur = 0
+ taille = int(finalTime / timeStep) +1
+ compteur = 0
+ dt = np.zeros([taille], dtype=PARMES_REAL, order=ORDER)
+ errEuler = np.zeros([taille], dtype=PARMES_REAL, order=ORDER)
+ errRK2 = np.zeros([taille], dtype=PARMES_REAL, order=ORDER)
+ errRK3 = np.zeros([taille], dtype=PARMES_REAL, order=ORDER)
+ errRK4 = np.zeros([taille], dtype=PARMES_REAL, order=ORDER)
+ UEuler = np.asarray([np.zeros([nb, nb, nb], dtype=PARMES_REAL,
+ order=ORDER) for d in xrange(3)])
+ URK2 = np.asarray([np.zeros([nb, nb, nb], dtype=PARMES_REAL,
+ order=ORDER) for d in xrange(3)])
+ URK3 = np.asarray([np.zeros([nb, nb, nb], dtype=PARMES_REAL,
+ order=ORDER) for d in xrange(3)])
+ URK4 = np.asarray([np.zeros([nb, nb, nb], dtype=PARMES_REAL,
+ order=ORDER) for d in xrange(3)])
+
+ UEuler[:,:,:,:] = 1.
+ URK2[:,:,:,:] = 1.
+ URK3[:,:,:,:] = 1.
+ URK4[:,:,:,:] = 1.
+
+ Analvorty[:,:,:,:]=np.vectorize(self.analyticalSolution)(timeStep,
+ topo3D.mesh.coords[0], \
+ topo3D.mesh.coords[1], \
+ topo3D.mesh.coords[2])[:]
+
+ while self.t < finalTime :
+ print "t", self.t
+
+ # Euler
+ self.method= Euler
+ methodInt=self.method(self.f)
+ Result = methodInt.integrate(self.f, self.t , timeStep , UEuler, topo3D)
+ UEuler = Result
+ errEuler[compteur]= np.max (abs(Analvorty[:,:,:,:] - Result[:,:,:,:]))
+
+ # RK2
+ self.method= RK2
+ methodInt=self.method(self.f)
+ Result = methodInt.integrate(self.f, self.t , timeStep , URK2, topo3D)
+ URK2 = Result
+ errRK2[compteur]= np.max (abs(Analvorty[:,:,:,:] - Result[:,:,:,:]))
+
+ # RK3
+ self.method= RK3
+ methodInt=self.method(self.f)
+ Result = methodInt.integrate(self.f, self.t , timeStep , URK3, topo3D)
+ URK3 = Result
+ errRK3[compteur]= np.max (abs(Analvorty[:,:,:,:] - Result[:,:,:,:]))
+
+ # RK4
+ self.method= RK4
+ methodInt=self.method(self.f)
+ Result = methodInt.integrate(self.f, self.t , timeStep , URK4, topo3D)
+ URK4 = Result
+ errRK4[compteur]= np.max (abs(Analvorty[:,:,:,:] - Result[:,:,:,:]))
+ dt[compteur] = self.t
+
+ self.t = self.t + timeStep
+ Analvorty[:,:,:,:]=np.vectorize(self.analyticalSolution)(self.t + timeStep,
+ topo3D.mesh.coords[0], \
+ topo3D.mesh.coords[1], \
+ topo3D.mesh.coords[2])[:]
+ compteur = compteur + 1
+
+ npt.assert_array_less(errEuler,timeStep)
+ npt.assert_array_less(errRK2,timeStep**2)
+ npt.assert_array_less(errRK3,timeStep**3)
+ npt.assert_array_less(errRK4,timeStep**4)
+
+ ## Table of error values
+# print "erreur Euler", errEuler
+# print "erreur RK2", errRK2
+# print "erreur RK3", errRK3
+# print "erreur RK4", errRK4
+
+ ## Plot error_scheme vs time
+# plt.figure(1)
+# plt.subplot(211)
+# plt.xlabel('dt')
+# plt.xscale('log')
+# plt.ylabel('Erreur')
+# plt.yscale('log')
+# plt.plot(dt, errEuler, '+-',dt ,errRK2, '+-' ,dt, errRK3, '+-',dt,errRK4, '+-')
+# plt.legend([u"Euler", u"RK2", u"RK3", u"RK4"], loc=4)
+# plt.subplot(212)
+# plt.xlabel('dt')
+# plt.ylabel('Erreur')
+# plt.plot(dt, errEuler, '+-',dt ,errRK2, '+-' ,dt, errRK3, '+-',dt,errRK4, '+-')
+# plt.legend([u"Euler", u"RK2", u"RK3", u"RK4"], loc=4)
+
+ plt.show()
+
+ t1 = time.time()
+
+ print "\n"
+# print "Total time : ", tf - t0, "sec (CPU)"
+# print "Init time : ", t1 - t0, "sec (CPU)"
+# print "Solving time : ", tf - t1, "sec (CPU)"
+
+
+ def runTest(self):
+ self.testIntegratorEDO()
+
+def suite():
+ suite = unittest.TestSuite()
+ suite.addTest(unittest.makeSuite(test_EDO))
+ return suite
+
+if __name__ == "__main__":
+ unittest.TextTestRunner(verbosity=2).run(suite())
diff --git a/HySoP/hysop/test/test_particular_solvers/test_RK.py b/HySoP/hysop/test/test_particular_solvers/test_RK.py
index 7281ecb4e..70ed6aa89 100644
--- a/HySoP/hysop/test/test_particular_solvers/test_RK.py
+++ b/HySoP/hysop/test/test_particular_solvers/test_RK.py
@@ -2,6 +2,9 @@
import time
import parmepy as pp
from parmepy.particular_solvers.integrator.euler import Euler
+from parmepy.particular_solvers.integrator.runge_kutta2 import RK2
+from parmepy.particular_solvers.integrator.runge_kutta3 import RK3
+from parmepy.particular_solvers.integrator.runge_kutta4 import RK4
from parmepy.particular_solvers.integrator.runge_kutta2stretching import RK2Stretch
from parmepy.particular_solvers.integrator.runge_kutta3stretching import RK3Stretch
from parmepy.particular_solvers.integrator.runge_kutta4stretching import RK4Stretch
@@ -15,7 +18,7 @@ from parmepy.constants import *
-class test_Stretching(unittest.TestCase):
+class test_RK(unittest.TestCase):
"""
DiscreteVariable test class
"""
@@ -27,43 +30,37 @@ class test_Stretching(unittest.TestCase):
return vx, vy, vz
def vorticite(self, x, y, z):
- wx = self.t+1.
- wy = self.t+1.
- wz = self.t+1.
+ wx = self.t**2+1.
+ wy = self.t**2+1.
+ wz = self.t**2+1.
return wx, wy, wz
- def analyticalSolution(self, x, y, z):
- sx = (0.5*(self.t*self.t)+self.t) * np.cos(x)
- sy = (0.5*(self.t*self.t)+self.t) * np.cos(y)
- sz = (0.5*(self.t*self.t)+self.t) * np.cos(z)
+ def analyticalSolution(self,t, x, y, z):
+ sx = t**2 - (t + t**3 /3.) * np.cos(x)
+ sy = t**2 - (t + t**3 /3.) * np.cos(y)
+ sz = t**2 - (t + t**3 /3.) * np.cos(z)
return sx, sy, sz
- def correctionRK2(self,x, y ,z):
- cx = np.cos(x) * np.cos(x)* self.t*self.t*( ( 1./2. + self.t/4. ) )
- cy = np.cos(y) * np.cos(y)* self.t*self.t*( ( 1./2. + self.t/4. ) )
- cz = np.cos(z) * np.cos(z)* self.t*self.t*( ( 1./2. + self.t/4. ) )
- return cx, cy, cz
-
- def correctionRK3(self,x, y ,z):
- cx = np.cos(x) * np.cos(x)* self.t*self.t*( ( 1./2. + self.t/2. ) + np.cos(x) * self.t *(1./6. + self.t/6.) )
- cy = np.cos(y) * np.cos(y)* self.t*self.t*( ( 1./2. + self.t/2. ) + np.cos(y) * self.t *(1./6. + self.t/6.) )
- cz = np.cos(z) * np.cos(z)* self.t*self.t*( ( 1./2. + self.t/2. ) + np.cos(z) * self.t *(1./6. + self.t/6.) )
- return cx, cy, cz
-
- def correctionRK4(self,x, y ,z):
- cx = 1./6. * (-np.cos(x) * (self.t)**2 + (np.cos(x))**2 * (self.t)**2 * (3. + 7./2. * self.t) + \
- 1./2. * (self.t)**3 * (np.cos(x))**3 * (2. + 3. * self.t) + 1./4. * (self.t)**4 * (np.cos(x))**4 * (1. + 2. * self.t))
- cy = 1./6. * (-np.cos(y) * (self.t)**2 + (np.cos(y))**2 * (self.t)**2 * (3. + 7./2. * self.t) + \
- 1./2. * (self.t)**3 * (np.cos(y))**3 * (2. + 3. * self.t) + 1./4. * (self.t)**4 * (np.cos(y))**4 * (1. + 2. * self.t))
- cz = 1./6. * (-np.cos(z) * (self.t)**2 + (np.cos(z))**2 * (self.t)**2 * (3. + 7./2. * self.t) + \
- 1./2. * (self.t)**3 * (np.cos(z))**3 * (2. + 3. * self.t) + 1./4. * (self.t)**4 * (np.cos(z))**4 * (1. + 2. * self.t))
- return cx, cy, cz
-
- def testOperatorStretching(self):
+ def f (self, t, u):
+ fx = 2.*t -u[0][:,:,:]
+ fy = 2.*t -u[1][:,:,:]
+ fz = 2.*t -u[2][:,:,:]
+ return fx , fy ,fz
+
+
+ def wgradu (self, t, x, y, z) :
+ wx =( 1 + t**2 )* np.cos(x)
+ wy =( 1 + t**2 )* np.cos(y)
+ wz =( 1 + t**2 )* np.cos(z)
+ return wx, wy, wz
+
+
+ def testIntegratorRK(self):
# Parameters
- nb = 64
+ nb = 32
timeStep = 0.01
finalTime = 0.01
+ maxerror=0.
t0 = time.time()
self.t = 0.
@@ -74,70 +71,69 @@ class test_Stretching(unittest.TestCase):
#####################################################################################
## Fields
- velo = pp.AnalyticalField(domain=box, formula=self.vitesse, name='Velocity', vector=True)
- vorti = pp.AnalyticalField(domain=box, formula=self.vorticite, name='Vorticity', vector=True)
+# velo = pp.AnalyticalField(domain=box, formula=self.vitesse, name='Velocity', vector=True)
+# vorti = pp.AnalyticalField(domain=box, formula=self.vorticite, name='Vorticity', vector=True)
## Operators
- stretch = pp.Stretching(velo,vorti)
+# stretch = pp.Stretching(velo,vorti)
## Solver creation (discretisation of objects is done in solver initialisation)
topo3D = pp.CartesianTopology(domain=box, resolution=[nb, nb, nb], dim=3)
##Problem
- pb = pp.Problem(topo3D, [stretch])
+# pb = pp.Problem(topo3D, [stretch])
## Setting solver to Problem
- pb.setSolver(finalTime, timeStep, solver_type='basic')
+# pb.setSolver(finalTime, timeStep, solver_type='basic')
+
+ Result=np.asarray([np.zeros([nb, nb, nb], dtype=PARMES_REAL, order=ORDER) for d in xrange(3)])
Analvorty=np.asarray([np.zeros([nb, nb, nb], dtype=PARMES_REAL, order=ORDER) for d in xrange(3)])
- Analvorty[:,:,:,:]= np.vectorize(self.vorticite)(topo3D.mesh.coords[0], \
- topo3D.mesh.coords[1], \
- topo3D.mesh.coords[2])[:]
- Analvorty2=np.asarray([np.zeros([nb, nb, nb], dtype=PARMES_REAL, order=ORDER) for d in xrange(3)])
- Analvorty2[:,:,:,:]= np.vectorize(self.vorticite)(topo3D.mesh.coords[0], \
+ Analvorty[:,:,:,:]= np.vectorize(self.wgradu)(self.t,topo3D.mesh.coords[0], \
topo3D.mesh.coords[1], \
topo3D.mesh.coords[2])[:]
- pb.solver.ODESolver= RK3Stretch #RK2Stretch #Euler#RK4Stretch #
+ self.method= RK3 # RK3Stretch #RK2Stretch #RK4Stretch #Euler #RK4 #
+
+ print "\nODE SOLVER :: ", self.method
+
+ methodInt=self.method(self.f)
+ Result = methodInt.integrate(self.f, self.t , timeStep , self.wgradu, topo3D)
- print "\nODE SOLVER :: ", pb.solver.ODESolver
- pb.initSolver()
+ if Euler == self.method :
+ maxerror = abs((timeStep**2 - 1./3. * timeStep**3)-.0) + timeStep**3
+ if RK2 == self.method :
+ maxerror = abs((-1./3. * timeStep**3) - (0.5* timeStep**2)) + timeStep**8
+ if RK3 == self.method :
+ maxerror = abs(0. - (0.5* timeStep**2 - 1./6. * timeStep**3)) + timeStep**3
+ if RK4 == self.method :
+ maxerror = abs((timeStep**2 - 1./3. * timeStep**3) - (7./6.* timeStep**2 - 1./2. * timeStep**3 + 1./8. * timeStep**4))
+
+
+# pb.initSolver()
t1 = time.time()
## Solve problem
- pb.solve()
- self.t = pb.timer.t
- Analvorty[:,:,:,:]=Analvorty[:,:,:,:]+np.vectorize(self.analyticalSolution)(topo3D.mesh.coords[0], \
+# pb.solve()
+ self.t = self.t + timeStep
+ Analvorty[:,:,:,:]=Analvorty[:,:,:,:]+np.vectorize(self.analyticalSolution)(self.t, topo3D.mesh.coords[0], \
topo3D.mesh.coords[1], \
topo3D.mesh.coords[2])[:]
- Analvorty2[:,:,:,:]=Analvorty2[:,:,:,:]+np.vectorize(self.analyticalSolution)(topo3D.mesh.coords[0], \
- topo3D.mesh.coords[1], \
- topo3D.mesh.coords[2])[:] + np.vectorize(self.correctionRK3)(topo3D.mesh.coords[0], \
- topo3D.mesh.coords[1], \
- topo3D.mesh.coords[2])[:]
-
-
## Visual comparison between the numerical resolution and the analytical resolution
for i in range(1):
plt.figure(1)
- plt.subplot(221)
+ plt.subplot(211)
plt.axis([-2*np.pi,2*np.pi,np.min(Analvorty[i][:,0,0]),np.max(Analvorty[i][:,0,0])])
- plt.plot(topo3D.mesh.coords[0][:,0,0], Analvorty[i][:,0,0], '-' ,topo3D.mesh.coords[0][:,0,0],vorti.discreteField[i].data[i][:,0,0], '-' )
- plt.legend([u"Solution Analytique", u"Solution Numérique"])
- plt.subplot(223)
- plt.axis([-2*np.pi,2*np.pi,0,max(abs(Analvorty[i][:,0,0] - vorti.discreteField[i].data[i][:,0,0]))])
- plt.plot(topo3D.mesh.coords[0][:,0,0], abs(Analvorty[i][:,0,0] - vorti.discreteField[i].data[i][:,0,0]), '-' )
- plt.legend([u"Erreur"])
- plt.subplot(222)
- plt.axis([-2*np.pi,2*np.pi,np.min(Analvorty2[i][:,0,0]),np.max(Analvorty2[i][:,0,0])])
- plt.plot(topo3D.mesh.coords[0][:,0,0], Analvorty2[i][:,0,0], '-' ,topo3D.mesh.coords[0][:,0,0],vorti.discreteField[i].data[i][:,0,0], '-' )
+ plt.plot(topo3D.mesh.coords[0][:,0,0], Analvorty[i][:,0,0], '-' ,topo3D.mesh.coords[0][:,0,0],Result[i][:,0,0], '-' )
plt.legend([u"Solution Analytique", u"Solution Numérique"])
- plt.subplot(224)
- plt.axis([-2*np.pi,2*np.pi,0,max(abs(Analvorty2[i][:,0,0] - vorti.discreteField[i].data[i][:,0,0]))])
- plt.plot(topo3D.mesh.coords[0][:,0,0], abs(Analvorty2[i][:,0,0] - vorti.discreteField[i].data[i][:,0,0]), '-' )
+ plt.subplot(212)
+ plt.axis([-2*np.pi,2*np.pi,0,max(abs(Analvorty[i][:,0,0] - Result[i][:,0,0]))])
+ plt.plot(topo3D.mesh.coords[0][:,0,0], abs(Analvorty[i][:,0,0] - Result[i][:,0,0]), '-' )
plt.legend([u"Erreur"])
+ plt.ax = plt.gca()
+ plt.ax.ticklabel_format(style='sci', axis='y')
plt.show()
- npt.assert_array_less(abs(Analvorty[:,:,:,:] - vorti.discreteField[0].data[:,:,:]),0.002)
+ npt.assert_array_less(abs(Analvorty[:,:,:,:] - Result[:,:,:,:]),maxerror)
tf = time.time()
print "\n"
@@ -147,11 +143,11 @@ class test_Stretching(unittest.TestCase):
def runTest(self):
- self.testOperatorStretching()
+ self.testIntegratorRK()
def suite():
suite = unittest.TestSuite()
- suite.addTest(unittest.makeSuite(test_Stretching))
+ suite.addTest(unittest.makeSuite(test_RK))
return suite
if __name__ == "__main__":
--
GitLab