From ba842990c16f81cc734883e28726942c4748bd2f Mon Sep 17 00:00:00 2001
From: Chloe Mimeau <Chloe.Mimeau@imag.fr>
Date: Tue, 28 May 2013 13:48:20 +0000
Subject: [PATCH] Update arrays operations in stretching time integration +
momentally fix origin problem in printer
---
HySoP/hysop/numerics/differentialOperator.py | 67 ++++++++++++++-----
HySoP/hysop/numerics/fct2op.py | 4 +-
HySoP/hysop/numerics/integrators/euler.py | 10 ++-
.../numerics/integrators/runge_kutta2.py | 22 ++++--
.../numerics/integrators/runge_kutta3.py | 32 ++++++---
.../numerics/integrators/runge_kutta4.py | 39 ++++++++---
HySoP/hysop/operator/discrete/stretching.py | 21 +++---
HySoP/hysop/operator/monitors/printer.py | 3 +-
8 files changed, 141 insertions(+), 57 deletions(-)
diff --git a/HySoP/hysop/numerics/differentialOperator.py b/HySoP/hysop/numerics/differentialOperator.py
index 9fc8b2fa9..9c8fca832 100755
--- a/HySoP/hysop/numerics/differentialOperator.py
+++ b/HySoP/hysop/numerics/differentialOperator.py
@@ -5,7 +5,7 @@ Discrete stretching representation
"""
from parmepy.constants import np, PARMES_REAL, ORDER
from parmepy.numerics.method import NumMethod
-#from parmepy.tools.cpu_data_transfer import Synchronize
+from parmepy.tools.cpu_data_transfer import Synchronize
class DifferentialOperator(NumMethod):
@@ -20,8 +20,8 @@ class DifferentialOperator(NumMethod):
Create a Stretching operator on a given continuous domain.
Work on a given field2 and field1 to compute the stretching term.
- @param field1 : field n°1 (vorticity).
- @param field2 : field n°2 (velocity or density).
+ @param field1 : field n1 (vorticity).
+ @param field2 : field n2 (velocity or density).
@param method : name and order of the spacial discretization method.
@param choice : differential operator to discretize
"""
@@ -31,7 +31,7 @@ class DifferentialOperator(NumMethod):
self.choice = choice
self.method = method
self.compute_time = 0.
- self.topology = self.field2.topology
+ self.topology = self.field1.topology
self.meshSize = self.topology.mesh.space_step
# def setUp(self):
@@ -271,14 +271,14 @@ class DifferentialOperator(NumMethod):
elif self.choice.find('gradV') >= 0:
# Ghosts synchronization
-# OpSynchronize = Synchronize(self.topology)
-# OpSynchronize.apply(self.field2)
+ OpSynchronize = Synchronize(self.topology)
+ OpSynchronize.apply(self.field2)
if self.method.find('FD_order2') >= 0:
raise ValueError("2nd order scheme Not yet implemented")
else:
- maxgersh = np.zeros(2, dtype=PARMES_REAL, order=ORDER)
+ maxArray = np.zeros(5, dtype=PARMES_REAL, order=ORDER)
## Fourth order scheme
# X components of temp and result
@@ -292,6 +292,8 @@ class DifferentialOperator(NumMethod):
1.0 * self.field2[0][ind0+2, ind1a:ind1b, ind2a:ind2b]
) / (12. * self.meshSize[0])
+ max1 = np.max(abs(temp1))
+
# temp2 = np.zeros(
# (self.resolution), dtype=PARMES_REAL, order=ORDER)
temp2 = self.field2[0][...] * 0.
@@ -302,6 +304,8 @@ class DifferentialOperator(NumMethod):
1.0 * self.field2[0][ind0a:ind0b, ind1+2, ind2a:ind2b]
) / (12. * self.meshSize[1])
+ max2 = np.max(abs(temp2))
+
# temp3 = np.zeros(
# (self.resolution), dtype=PARMES_REAL, order=ORDER)
temp3 = self.field2[0][...] * 0.
@@ -314,6 +318,7 @@ class DifferentialOperator(NumMethod):
maxstr1 = np.max(abs(temp1) + abs(temp2) + abs(temp3))
maxadv1 = np.max(abs(temp1))
+ max3 = np.max(abs(temp3))
# Y components of temp and result
# temp4 = np.zeros(
@@ -326,6 +331,8 @@ class DifferentialOperator(NumMethod):
1.0 * self.field2[1][ind0+2, ind1a:ind1b, ind2a:ind2b]
) / (12. * self.meshSize[0])
+ max4 = np.max(abs(temp4))
+
# temp5 = np.zeros(
# (self.resolution), dtype=PARMES_REAL, order=ORDER)
temp5 = self.field2[1][...] * 0.
@@ -336,6 +343,8 @@ class DifferentialOperator(NumMethod):
1.0 * self.field2[1][ind0a:ind0b, ind1+2, ind2a:ind2b]
) / (12. * self.meshSize[1])
+ max5 = np.max(abs(temp5))
+
# temp6 = np.zeros(
# (self.resolution), dtype=PARMES_REAL, order=ORDER)
temp6 = self.field2[1][...] * 0.
@@ -348,6 +357,7 @@ class DifferentialOperator(NumMethod):
maxstr2 = np.max(abs(temp4)+abs(temp5)+abs(temp6))
maxadv2 = np.max(abs(temp5))
+ max6 = np.max(abs(temp6))
# Z components of temp and result
# temp7 = np.zeros(
@@ -360,6 +370,8 @@ class DifferentialOperator(NumMethod):
1.0 * self.field2[2][ind0+2, ind1a:ind1b, ind2a:ind2b]
) / (12. * self.meshSize[0])
+ max7 = np.max(abs(temp7))
+
# temp8 = np.zeros(
# (self.resolution), dtype=PARMES_REAL, order=ORDER)
temp8 = self.field2[2][...] * 0.
@@ -370,6 +382,8 @@ class DifferentialOperator(NumMethod):
1.0 * self.field2[2][ind0a:ind0b, ind1+2, ind2a:ind2b]
) / (12. * self.meshSize[1])
+ max8 = np.max(abs(temp8))
+
# temp9 = np.zeros(
# (self.resolution), dtype=PARMES_REAL, order=ORDER)
temp9 = self.field2[2][...] * 0.
@@ -382,21 +396,40 @@ class DifferentialOperator(NumMethod):
maxstr3 = np.max(abs(temp7)+abs(temp8)+abs(temp9))
maxadv3 = np.max(abs(temp9))
- maxgersh[0] = max(maxstr1, maxstr2, maxstr3)
- maxgersh[1] = max(maxadv1, maxadv2, maxadv3)
+ max9 = np.max(abs(temp9))
+
+ # grad(V) result
result = np.concatenate((
np.array([temp1, temp2, temp3]),
np.array([temp4, temp5, temp6]),
np.array([temp7, temp8, temp9])
))
-
- return result, maxgersh
+ # div(V) result
+ divergence = np.array(temp1 + temp5 + temp9)
+
+ maxdiv = np.max(abs(divergence))
+ maxMax = max(max1, max2, max3,
+ max4, max5, max6,
+ max7, max8, max9)
+ maxproj = maxdiv / maxMax
+
+ # maxima of partial derivatives : needed for stab conditions
+ # advection stab
+ maxArray[0] = max(maxstr1, maxstr2, maxstr3)
+ # stretching stab
+ maxArray[1] = max(maxadv1, maxadv2, maxadv3)
+ # solenoidal reprojection criterion
+ maxArray[2] = maxproj
+ maxArray[3] = maxdiv
+ maxArray[4] = maxMax
+
+ return result, maxArray, divergence
elif (self.choice.find('gradS') >= 0):
# Ghosts synchronization
-# OpSynchronize = SynchronizeS(self.topology)
-# OpSynchronize.apply(self.field2)
+ OpSynchronize = SynchronizeS(self.topology)
+ OpSynchronize.apply(self.field2)
if self.method.find('FD_order2') >= 0:
raise ValueError("2nd order scheme Not yet implemented")
@@ -434,8 +467,8 @@ class DifferentialOperator(NumMethod):
elif self.choice.find('laplacianV') >= 0:
# Ghosts synchronization
-# OpSynchronize = Synchronize(self.topology)
-# OpSynchronize.apply(self.field2)
+ OpSynchronize = Synchronize(self.topology)
+ OpSynchronize.apply(self.field2)
if self.method.find('FD_order2') >= 0:
raise ValueError("2nd order scheme Not yet implemented")
@@ -583,8 +616,8 @@ class DifferentialOperator(NumMethod):
elif self.choice.find('curl') >= 0:
# Ghosts synchronization
-# OpSynchronize = Synchronize(self.topology)
-# OpSynchronize.apply(self.field1)
+ OpSynchronize = Synchronize(self.topology)
+ OpSynchronize.apply(self.field1)
if self.method.find('FD_order2') >= 0:
raise ValueError("2nd order scheme Not yet implemented")
diff --git a/HySoP/hysop/numerics/fct2op.py b/HySoP/hysop/numerics/fct2op.py
index be96777fc..43b783d77 100644
--- a/HySoP/hysop/numerics/fct2op.py
+++ b/HySoP/hysop/numerics/fct2op.py
@@ -71,8 +71,8 @@ class Fct2Op(object):
y[1][ind0a:ind0b,ind1a:ind1b,ind2a:ind2b] + \
self.field2[8][ind0a:ind0b,ind1a:ind1b,ind2a:ind2b]* \
y[2][ind0a:ind0b,ind1a:ind1b,ind2a:ind2b]
- result = np.concatenate((np.array([temp0]),np.array([temp1]),np.array([temp2])))
- return result
+# result = np.concatenate((np.array([temp0]),np.array([temp1]),np.array([temp2])))
+ return temp0, temp1, temp2
if self.choice == 'divWU' : # field2 = velocity field
result = DifferentialOperator(y, self.field2, self.method, self.choice).__call__()
diff --git a/HySoP/hysop/numerics/integrators/euler.py b/HySoP/hysop/numerics/integrators/euler.py
index 8c02801bd..57bbb3f91 100755
--- a/HySoP/hysop/numerics/integrators/euler.py
+++ b/HySoP/hysop/numerics/integrators/euler.py
@@ -38,8 +38,14 @@ class Euler(ODESolver):
@param t : current time.
@param dt : time step.
# """
- result = y + dt * np.asarray(f(t, y), dtype=PARMES_REAL, order=ORDER)
- return result
+# result = y + dt * np.asarray(f(t, y), dtype=PARMES_REAL, order=ORDER)
+# return result
+
+ RHS = f(t, y)
+
+ return y[0] + dt * RHS[0], \
+ y[1] + dt * RHS[1], \
+ y[2] + dt * RHS[2]
def __str__(self):
"""ToString method"""
diff --git a/HySoP/hysop/numerics/integrators/runge_kutta2.py b/HySoP/hysop/numerics/integrators/runge_kutta2.py
index 8a66b71a6..3419c2e33 100755
--- a/HySoP/hysop/numerics/integrators/runge_kutta2.py
+++ b/HySoP/hysop/numerics/integrators/runge_kutta2.py
@@ -37,12 +37,22 @@ class RK2(ODESolver):
@param y : function Y(t,space)
"""
- # K1 = dt *f[t,y(t,x,y,z)]
- K1 = dt * np.asarray(f(t, y), dtype=PARMES_REAL, order=ORDER)
- # result = y(t,x,y,z)] + dt *f[t + dt/2 , y(t,x,y,z) + k1/2.]
- result = y + dt * np.asarray(f(t + dt / 2., y + K1 / 2.),
- dtype=PARMES_REAL, order=ORDER)
- return result
+# # K1 = dt *f[t,y(t,x,y,z)]
+# K1 = dt * np.asarray(f(t, y), dtype=PARMES_REAL, order=ORDER)
+# # result = y(t,x,y,z)] + dt *f[t + dt/2 , y(t,x,y,z) + k1/2.]
+
+# result = y + dt * np.asarray(f(t + dt / 2., y + K1 / 2.),
+# dtype=PARMES_REAL, order=ORDER)
+# return result
+
+ RHS1 = f(t, y)
+ K1 = np.array([dt * RHS1[0], dt * RHS1[1], dt * RHS1[2]])
+
+ RHS2 = f(t + dt / 2., y + K1 / 2.)
+
+ return y[0] + dt * RHS2[0], \
+ y[1] + dt * RHS2[1], \
+ y[2] + dt * RHS2[2]
def __str__(self):
"""ToString method"""
diff --git a/HySoP/hysop/numerics/integrators/runge_kutta3.py b/HySoP/hysop/numerics/integrators/runge_kutta3.py
index 454fc4604..c5de5014f 100755
--- a/HySoP/hysop/numerics/integrators/runge_kutta3.py
+++ b/HySoP/hysop/numerics/integrators/runge_kutta3.py
@@ -36,16 +36,30 @@ class RK3(ODESolver):
@param dt : time step.
@param y : function Y(t,space)
"""
- # K1 = dt *f[t,y(t,x,y,z)]
- K1 = dt * np.asarray(f(t, y), dtype=PARMES_REAL, order=ORDER)
- # result = y(t,x,y,z)] + dt *f[t + dt/2 , y(t,x,y,z) + k1/2.]
- K2 = dt * np.asarray(f(t + dt / 3., y + K1[:, :, :, :] / 3.),
- dtype=PARMES_REAL, order=ORDER)
- K3 = dt * np.asarray(f(t + dt * 2. / 3., y + K2[:, :, :, :] * 2. / 3.),
- dtype=PARMES_REAL, order=ORDER)
+# # K1 = dt *f[t,y(t,x,y,z)]
+# K1 = dt * np.asarray(f(t, y), dtype=PARMES_REAL, order=ORDER)
+# # result = y(t,x,y,z)] + dt *f[t + dt/2 , y(t,x,y,z) + k1/2.]
+# K2 = dt * np.asarray(f(t + dt / 3., y + K1[:, :, :, :] / 3.),
+# dtype=PARMES_REAL, order=ORDER)
+# K3 = dt * np.asarray(f(t + dt * 2. / 3., y + K2[:, :, :, :] * 2. / 3.),
+# dtype=PARMES_REAL, order=ORDER)
+# # result = y(t,x,y,z)] + 1/4 [ k1 + 3 k3]
+# result = y + 1. / 4. * (K1 + 3. * K3)
+# return result
+
+ RHS1 = f(t, y)
+ K1 = np.array([dt * RHS1[0], dt * RHS1[1], dt * RHS1[2]])
+
+ RHS2 = f(t + dt / 3., y + K1 / 3.)
+ K2 = np.array([dt * RHS2[0], dt * RHS2[1], dt * RHS2[2]])
+
+ RHS3 = f(t + dt * 2. / 3., y + K2 * 2. / 3.)
+ K3 = np.array([dt * RHS3[0], dt * RHS3[1], dt * RHS3[2]])
+
# result = y(t,x,y,z)] + 1/4 [ k1 + 3 k3]
- result = y + 1. / 4. * (K1 + 3. * K3)
- return result
+ return y[0] + 1. / 4. * (K1[0] + 3. * K3[0]), \
+ y[1] + 1. / 4. * (K1[1] + 3. * K3[1]), \
+ y[2] + 1. / 4. * (K1[2] + 3. * K3[2])
def __str__(self):
"""ToString method"""
diff --git a/HySoP/hysop/numerics/integrators/runge_kutta4.py b/HySoP/hysop/numerics/integrators/runge_kutta4.py
index 286d9fc1b..679352b67 100755
--- a/HySoP/hysop/numerics/integrators/runge_kutta4.py
+++ b/HySoP/hysop/numerics/integrators/runge_kutta4.py
@@ -38,18 +38,35 @@ class RK4(ODESolver):
@param y : function Y(t,space)
"""
- # K1 = dt *f[t, y(t, x, y, z)]
- K1 = dt * np.asarray(f(t, y), dtype=PARMES_REAL, order=ORDER)
- K2 = dt * np.asarray(f(t + dt / 2., y + K1[:, :, :, :] / 2.),
- dtype=PARMES_REAL, order=ORDER)
- K3 = dt * np.asarray(f(t + dt / 2., y + K2[:, :, :, :] / 2.),
- dtype=PARMES_REAL, order=ORDER)
- K4 = dt * np.asarray(f(t + dt, y + K3[:, :, :, :]),
- dtype=PARMES_REAL, order=ORDER)
+# # K1 = dt *f[t, y(t, x, y, z)]
+# K1 = dt * np.asarray(f(t, y), dtype=PARMES_REAL, order=ORDER)
+# K2 = dt * np.asarray(f(t + dt / 2., y + K1[:, :, :, :] / 2.),
+# dtype=PARMES_REAL, order=ORDER)
+# K3 = dt * np.asarray(f(t + dt / 2., y + K2[:, :, :, :] / 2.),
+# dtype=PARMES_REAL, order=ORDER)
+# K4 = dt * np.asarray(f(t + dt, y + K3[:, :, :, :]),
+# dtype=PARMES_REAL, order=ORDER)
+# # result = y(t, x, y, z)] + dt *f[t + dt/2 , y(t, x, y, z) + k1/2.]
+# result = y + 1. / 6. * (K1[:, :, :, :] + 2. * K2[:, :, :, :] +
+# 2. * K3[:, :, :, :] + K4[:, :, :, :])
+# return result
+
+ RHS1 = f(t, y)
+ K1 = np.array([dt * RHS1[0], dt * RHS1[1], dt * RHS1[2]])
+
+ RHS2 = f(t + dt / 2., y + K1 / 1.)
+ K2 = np.array([dt * RHS2[0], dt * RHS2[1], dt * RHS2[2]])
+
+ RHS3 = f(t + dt * 2., y + K2 * 2.)
+ K3 = np.array([dt * RHS3[0], dt * RHS3[1], dt * RHS3[2]])
+
+ RHS4 = f(t + dt, y + K3)
+ K4 = np.array([dt * RHS4[0], dt * RHS4[1], dt * RHS4[2]])
+
# result = y(t, x, y, z)] + dt *f[t + dt/2 , y(t, x, y, z) + k1/2.]
- result = y + 1. / 6. * (K1[:, :, :, :] + 2. * K2[:, :, :, :] +
- 2. * K3[:, :, :, :] + K4[:, :, :, :])
- return result
+ return y[0] + 1. / 6. * (K1[0] + 2. * K2[0] + 2. * K3[0] + K4[0]), \
+ y[1] + 1. / 6. * (K1[1] + 2. * K2[1] + 2. * K3[1] + K4[1]), \
+ y[2] + 1. / 6. * (K1[2] + 2. * K2[2] + 2. * K3[2] + K4[2])
def __str__(self):
"""ToString method"""
diff --git a/HySoP/hysop/operator/discrete/stretching.py b/HySoP/hysop/operator/discrete/stretching.py
index 227d9bb7a..3d55cbc29 100755
--- a/HySoP/hysop/operator/discrete/stretching.py
+++ b/HySoP/hysop/operator/discrete/stretching.py
@@ -60,9 +60,7 @@ class Stretching_d(DiscreteOperator):
else:
self.num_method = RK3
self.cststretch = 2.5127
- self.cststretch = self.cststretch / 2.
- @debug
@timed_function
def apply(self, simulation):
"""
@@ -106,13 +104,13 @@ class Stretching_d(DiscreteOperator):
method='FD_order4',
choice='gradV')
- gradResult, maxgersh = self.stretchOp()
+ gradResult, maxArray, divergence = self.stretchOp()
- print "Maxgersh", maxgersh
- print "dtAdvec", 0.5 / maxgersh[1]
+ print "maxadv, maxgersh", maxArray[1], maxArray[0]
+ print "dtAdvec", 0.5 / maxArray[1]
## Determination of Stretching time step (stability condition)
- self.ndt = self.stabilityTest(dt, maxgersh[0],
+ self.ndt = self.stabilityTest(dt, maxArray[0],
self.num_method, self.cststretch)
self.dtstretch = dt / float(self.ndt)
@@ -127,9 +125,14 @@ class Stretching_d(DiscreteOperator):
t = simulation.time
for i in range(self.ndt):
methodInt = self.num_method(self.function)
- self.vorticity.data[...] = methodInt(self.function.apply,
- t, self.dtstretch,
- self.vorticity.data)
+# self.vorticity.data[...] = methodInt(self.function.apply,
+# t, self.dtstretch,
+# self.vorticity.data)
+ self.vorticity.data[0][...], \
+ self.vorticity.data[1][...], \
+ self.vorticity.data[2][...] = methodInt(self.function.apply,
+ t, self.dtstretch,
+ self.vorticity.data)
# print 'norm omega', LA.norm(self.vorticity.data)
def stabilityTest(self, dt, maxi, num_method, cststretch):
diff --git a/HySoP/hysop/operator/monitors/printer.py b/HySoP/hysop/operator/monitors/printer.py
index 2e0f1c32d..0b7ae2597 100644
--- a/HySoP/hysop/operator/monitors/printer.py
+++ b/HySoP/hysop/operator/monitors/printer.py
@@ -117,7 +117,8 @@ class Printer(Monitoring):
print "Print to file " + filename
## VTK output \todo: Need fix in 2D, getting an IOError.
try:
- evtk.imageToVTK(filename, pointData=self._build_vtk_dict())
+ orig = self.variables[0].discreteFields.values()[0].topology.mesh.global_start
+ evtk.imageToVTK(filename, origin=(orig[0], orig[1], orig[2]), pointData=self._build_vtk_dict())
except IOError:
if self.variables[0].domain.dimension == 2:
## Standard output
--
GitLab