From b028cfff563166abca33b1e0fdb2b596b77b11b6 Mon Sep 17 00:00:00 2001
From: Eric Dalissier <Eric.Dalissier@imag.fr>
Date: Tue, 6 Nov 2012 14:03:57 +0000
Subject: [PATCH] forward euler test add
---
.../hysop/operator/differentialOperator_d.py | 34 +++++-
HySoP/hysop/operator/stretching_d.py | 3 +-
.../integrator/runge_kutta4.py | 32 ++---
.../test_particular_solvers/test_euler.py | 114 ++++++++++++++++++
4 files changed, 167 insertions(+), 16 deletions(-)
create mode 100644 HySoP/hysop/test/test_particular_solvers/test_euler.py
diff --git a/HySoP/hysop/operator/differentialOperator_d.py b/HySoP/hysop/operator/differentialOperator_d.py
index 739b6e4ab..1c3fd582d 100755
--- a/HySoP/hysop/operator/differentialOperator_d.py
+++ b/HySoP/hysop/operator/differentialOperator_d.py
@@ -100,7 +100,39 @@ class DifferentialOperator_d(DiscreteOperator):
self.result[2][...]= temp1 + temp2 +temp3
elif self.choice == 'curl':
- pass
+ temp1 = (1.0 * self.field1[2][:,np.roll(ind,+2),:] -
+ 8.0 * self.field1[2][:,np.roll(ind,+1),:] +\
+ 8.0 * self.field1[2][:,np.roll(ind,-1),:] -\
+ 1.0 * self.field1[2][:,np.roll(ind,-2),:] ) / (12. * mesh_size[1])
+
+ temp2 = (1.0 * self.field1[1][...,np.roll(ind,+2)] -\
+ 8.0 * self.field1[1][...,np.roll(ind,+1)] +\
+ 8.0 * self.field1[1][...,np.roll(ind,-1)] -\
+ 1.0 * self.field1[1][...,np.roll(ind,-2)] ) / (12. * mesh_size[2])
+ self.result[0]= temp1 - temp2
+
+ temp1 = (1.0 * self.field1[2][np.roll(ind,+2),...] -
+ 8.0 * self.field1[2][np.roll(ind,+1),...] +\
+ 8.0 * self.field1[2][np.roll(ind,-1),...] -\
+ 1.0 * self.field1[2][np.roll(ind,-2),...] ) / (12. * mesh_size[0])
+
+ temp2 = (1.0 * self.field1[0][...,np.roll(ind,+2)] -\
+ 8.0 * self.field1[0][...,np.roll(ind,+1)] +\
+ 8.0 * self.field1[0][...,np.roll(ind,-1)] -\
+ 1.0 * self.field1[0][...,np.roll(ind,-2)] ) / (12. * mesh_size[2])
+ self.result[1]= temp1 - temp2
+
+ temp1 = (1.0 * self.field1[1][...,np.roll(ind,+2)] -
+ 8.0 * self.field1[1][...,np.roll(ind,+1)] +\
+ 8.0 * self.field1[1][...,np.roll(ind,-1)] -\
+ 1.0 * self.field1[1][...,np.roll(ind,-2)] ) / (12. * mesh_size[0])
+
+ temp2 = (1.0 * self.field1[0][:,np.roll(ind,+2),:] -\
+ 8.0 * self.field1[0][:,np.roll(ind,+1),:] +\
+ 8.0 * self.field1[0][:,np.roll(ind,-1),:] -\
+ 1.0 * self.field1[0][:,np.roll(ind,-2),:] ) / (12. * mesh_size[1])
+ self.result[2]= temp1 - temp2
+
else:
raise ValueError("Unknown differiential operator: " + str(self.choice))
diff --git a/HySoP/hysop/operator/stretching_d.py b/HySoP/hysop/operator/stretching_d.py
index 51c0e221f..4b02ba8f4 100755
--- a/HySoP/hysop/operator/stretching_d.py
+++ b/HySoP/hysop/operator/stretching_d.py
@@ -66,7 +66,8 @@ class Stretching_d(DiscreteOperator):
self.stretchTerm.discreteOperator.apply()
# Time integration of dw/dt = div(wu)
- methodInt=self.method()
+ methodInt=self.method(self.stretchTerm.discreteOperator)
+ print dir(methodInt)
for i in range(self.vorticity.topology.dim):
self.res_vorticity[i][...] = methodInt.integrate(self.vorticity, self.stretchTermResult, t, dt, direction=i)
diff --git a/HySoP/hysop/particular_solvers/integrator/runge_kutta4.py b/HySoP/hysop/particular_solvers/integrator/runge_kutta4.py
index fe8f05794..b528ca612 100755
--- a/HySoP/hysop/particular_solvers/integrator/runge_kutta4.py
+++ b/HySoP/hysop/particular_solvers/integrator/runge_kutta4.py
@@ -15,7 +15,7 @@ class RK4(ODESolver):
y(t_n+1)= y(t_n) + dt*y'[y(t_n)+dt/2*y'(y(t_n))].
"""
- def __init__(self, f=None, conditions=lambda x: x, dim=3):
+ def __init__(self, f=None, dim=3):
"""
Constructor.
@@ -24,11 +24,7 @@ class RK4(ODESolver):
@param conditions : function to apply boundary conditions.
"""
ODESolver.__init__(self, f)
- # Boundary conditions function.
- self.boundaryConditions = conditions
- self.gpu_kernel = """
- """
- self.nb_flop = 5
+
def integrate(self, f, yp , t, dt):
"""
@@ -36,18 +32,26 @@ class RK4(ODESolver):
yp= y + dt*y'[y+dt/2*y'(y)].
@param y : position at time t.
- @param fy : y'(y) value at time t.
+ @param f : operator to calcul y'(y) value at time t.
@param yp : result position at time t + dt
@param t : current time.
@param dt : time step.
"""
- K1 = yp + dt * f(t, yp)
- K2 = yp + dt * f(t + 0.5 * dt, (yp + K1) * 0.5)
- K3 = yp + dt * f(t + 0.5 * dt, (yp + K2) * 0.5)
- return (yp + (1/6.0) * dt * ( f(t, yp) +
- 2. * f(t + 0.5 * dt, (yp + K1) * 0.5) +
- 2. * f(t + 0.5 * dt, (yp + K2) * 0.5) +
- f(t + dt, K3 )) )
+ dt2 = dt / 2.0
+ K1 = f.field1() * dt
+ yt = f(field1= yp + K1*0.5, f.field2() )
+ yt.apply()
+ K2 = yt.field1() * dt
+ yt = f(field1= yp + K2*0.5, f.field2() )
+ yt.apply()
+ K3 = yt.field1() * dt
+ yt = f(field1= yp + K3, f.field2() )
+ yt.apply()
+ K4 = yt.field1() * dt
+ unew = yp + (1/6.0) * (K1 + K2 *2. + K3 * 2. + K4)
+ return unew
+
+
def __str__(self):
"""ToString method"""
diff --git a/HySoP/hysop/test/test_particular_solvers/test_euler.py b/HySoP/hysop/test/test_particular_solvers/test_euler.py
new file mode 100644
index 000000000..12b0dac2d
--- /dev/null
+++ b/HySoP/hysop/test/test_particular_solvers/test_euler.py
@@ -0,0 +1,114 @@
+# -*- coding: utf-8 -*-
+import unittest
+from parmepy.fields.discrete import *
+from parmepy.fields.continuous import *
+from parmepy.fields.analytical import *
+from parmepy.domain.topology import *
+from parmepy.domain.box import *
+from parmepy.constants import *
+from parmepy.particular_solvers.integrator.integrator import ODESolver
+from parmepy.particular_solvers.integrator.euler import Euler
+import numpy as np
+import numpy.testing as npt
+import math
+
+
+class test_Euler(unittest.TestCase):
+ """
+ Euler test
+ solve u'(t) = f( u(t),t )
+ u(t) = a t + U0
+ f(u,t) = a + (u - (at + U0))^4
+ """
+ def setUp(self):
+ self.e = 0.0002 # Accepted error between result and analytical result
+ self.dim = 3
+ self.boxLength = [1., 1., 1.]
+ self.boxMin = [ 0., 0., 0.]
+ self.nbPts = [6, 6, 6]
+ self.box = Box(dimension=self.dim,
+ length=self.boxLength,
+ origin=self.boxMin)
+ self.dt = 0.1
+ self.FinalTime = 3
+ self.U0 = [1.0 , 2.0, 3.0 ]
+
+ def testEulerInt(self):
+ # Continuous fields and operator declaration
+ self.velo = AnalyticalField(domain=self.box, formula=self.vitesse, name='Velocity', vector=True)
+ # Topology definition / Fields and operator discretization
+ self.result = [np.ones((self.nbPts), dtype=PARMES_REAL, order=ORDER) for d in xrange(self.dim)]
+ self.topo3D = CartesianTopology(domain=self.box, resolution=self.nbPts, dim=self.dim)
+ self.velo.discretize(self.topo3D)
+ self.velo.initialize()
+ self.resol = self.topo3D.mesh.resolution
+ #test for t=0 and initialization of the fonction f
+ t = 0.
+ fctInter = [np.ones((self.nbPts), dtype=PARMES_REAL, order=ORDER) for d in xrange(self.dim)]
+ method = Euler()
+
+ #Solution calculated with Euler
+ fctInter = self.fctTest(t,self.velo.discreteField[self.velo._fieldId].data[0]\
+ ,self.velo.discreteField[self.velo._fieldId].data[1]\
+ ,self.velo.discreteField[self.velo._fieldId].data[2])
+ for d in xrange(self.dim) :
+ self.result[d][...] = method.integrate(self.velo.discreteField[self.velo._fieldId], fctInter, t, self.dt, d )
+
+ t=t+self.dt
+
+ #Analytical solution
+ self.anal=self.analyticalSolution(t,self.velo.discreteField[self.velo._fieldId].data[0]\
+ ,self.velo.discreteField[self.velo._fieldId].data[1]\
+ ,self.velo.discreteField[self.velo._fieldId].data[2])
+ npt.assert_array_less(abs(self.anal[0]-self.result[0]), self.e)
+ npt.assert_array_less(abs(self.anal[1]-self.result[1]), self.e)
+ npt.assert_array_less(abs(self.anal[2]-self.result[2]), self.e)
+
+ # time loop
+ while t < self.FinalTime :
+ #print "T=" , t
+ fctInter = self.fctTest(t,self.result[0],self.result[1],self.result[2])
+ for d in xrange(self.dim) :
+ self.result[d][...] = method.integrate(self.result, fctInter, t, self.dt, d )
+
+ t = t + self.dt
+
+ self.anal=self.analyticalSolution(t,self.velo.discreteField[self.velo._fieldId].data[0]\
+ ,self.velo.discreteField[self.velo._fieldId].data[1]\
+ ,self.velo.discreteField[self.velo._fieldId].data[2])
+ # Comparison with analytical solution
+ npt.assert_array_less(abs(self.anal[0]-self.result[0]), self.e)
+ npt.assert_array_less(abs(self.anal[1]-self.result[1]), self.e)
+ npt.assert_array_less(abs(self.anal[2]-self.result[2]), self.e)
+
+
+ def vitesse(self, x, y, z):
+ vx = 1.0
+ vy = 2.0
+ vz = 3.0
+ return vx, vy, vz
+
+ def fctTest(self, t, x, y, z):
+ wx = 0.2 +( x - self.analyticalSolution(t, 1., 2., 3.)[0])**4
+ wy = 0.5 +( y - self.analyticalSolution(t, 1., 2., 3.)[1])**4
+ wz = 1. +( z - self.analyticalSolution(t, 1., 2., 3.)[2])**4
+ return [wx,wy,wz]
+
+ def analyticalSolution(self, t, x, y, z):
+ sx = 0.2 *t + x
+ sy = 0.5 *t + y
+ sz = 1. *t + z
+ return [sx, sy, sz]
+
+ def runTest(self):
+ self.setUp()
+ self.testOperatorDiv()
+
+
+def suite():
+ suite = unittest.TestSuite()
+ suite.addTest(unittest.makeSuite(test_Euler))
+ return suite
+
+if __name__ == "__main__":
+ unittest.TextTestRunner(verbosity=2).run(suite())
--
GitLab