From f950b51e8fbe200d5a99a464ce67e6a7c8ad7617 Mon Sep 17 00:00:00 2001
From: Chloe Mimeau <Chloe.Mimeau@imag.fr>
Date: Wed, 3 Jul 2013 07:06:09 +0000
Subject: [PATCH] updates
---
Examples/Attic/NavierStokes3d_vortRing.py | 108 ++++++++++---------
HySoP/hysop/numerics/differentialOperator.py | 45 +++-----
HySoP/hysop/numerics/fct2op.py | 2 +-
HySoP/hysop/operator/discrete/multiphase.py | 88 +++++++++------
HySoP/hysop/operator/discrete/poisson_fft.py | 2 +-
HySoP/hysop/operator/discrete/stretching.py | 2 +-
6 files changed, 131 insertions(+), 116 deletions(-)
diff --git a/Examples/Attic/NavierStokes3d_vortRing.py b/Examples/Attic/NavierStokes3d_vortRing.py
index 52e4a1cef..c609bf47d 100755
--- a/Examples/Attic/NavierStokes3d_vortRing.py
+++ b/Examples/Attic/NavierStokes3d_vortRing.py
@@ -19,6 +19,7 @@ from parmepy.problem.simulation import Simulation
from parmepy.operator.monitors.energy_enstrophy import Energy_enstrophy
from parmepy.operator.monitors.reprojection_criterion import Reprojection_criterion
from parmepy.operator.monitors.printer import Printer
+from parmepy.dataloader import DataLoader
#from numpy import linalg as LA
@@ -32,48 +33,51 @@ print "Mpi process number ", rank
print " ========= Start test for Navier-Stokes 3D (Taylor Green benchmark)========="
## Opening/reading input data file
-inputData = open("input.dat", "r")
-
-Ox = np.float64(inputData.readline().rstrip('\n\r'))
-Oy = np.float64(inputData.readline().rstrip('\n\r'))
-Oz = np.float64(inputData.readline().rstrip('\n\r'))
-Lx = np.float64(inputData.readline().rstrip('\n\r'))
-Ly = np.float64(inputData.readline().rstrip('\n\r'))
-Lz = np.float64(inputData.readline().rstrip('\n\r'))
-resX = np.uint32(inputData.readline().rstrip('\n\r'))
-resY = np.uint32(inputData.readline().rstrip('\n\r'))
-resZ = np.uint32(inputData.readline().rstrip('\n\r'))
-dt = np.float64(inputData.readline().rstrip('\n\r'))
-finalTime = np.float64(inputData.readline().rstrip('\n\r'))
-restart = np.uint32(inputData.readline().rstrip('\n\r'))
-dumpfreq = np.uint32(inputData.readline().rstrip('\n\r'))
-dumpfile = inputData.readline().rstrip('\n\r')
-outModuloPrinter = np.uint32(inputData.readline().rstrip('\n\r'))
-outputPrinter = inputData.readline().rstrip('\n\r')
-outModuloForces = np.uint32(inputData.readline().rstrip('\n\r'))
-outputForces = inputData.readline().rstrip('\n\r')
-outModuloEnergies = np.uint32(inputData.readline().rstrip('\n\r'))
-outputEnergies = inputData.readline().rstrip('\n\r')
-visco = np.float64(inputData.readline().rstrip('\n\r'))
-vortProj = np.uint32(inputData.readline().rstrip('\n\r'))
-methodAdv = inputData.readline().rstrip('\n\r')
-methodStretch = inputData.readline().rstrip('\n\r')
-multires = np.uint32(inputData.readline().rstrip('\n\r'))
-resFilterX = np.uint32(inputData.readline().rstrip('\n\r'))
-resFilterY = np.uint32(inputData.readline().rstrip('\n\r'))
-resFilterZ = np.uint32(inputData.readline().rstrip('\n\r'))
-
-inputData.close()
+data = DataLoader('inputData_TG.dat')
+
+#inputData = open("input.dat", "r")
+
+#Ox = np.float64(inputData.readline().rstrip('\n\r'))
+#Oy = np.float64(inputData.readline().rstrip('\n\r'))
+#Oz = np.float64(inputData.readline().rstrip('\n\r'))
+#Lx = np.float64(inputData.readline().rstrip('\n\r'))
+#Ly = np.float64(inputData.readline().rstrip('\n\r'))
+#Lz = np.float64(inputData.readline().rstrip('\n\r'))
+#resX = np.uint32(inputData.readline().rstrip('\n\r'))
+#resY = np.uint32(inputData.readline().rstrip('\n\r'))
+#resZ = np.uint32(inputData.readline().rstrip('\n\r'))
+#dt = np.float64(inputData.readline().rstrip('\n\r'))
+#finalTime = np.float64(inputData.readline().rstrip('\n\r'))
+#restart = np.uint32(inputData.readline().rstrip('\n\r'))
+#dumpfreq = np.uint32(inputData.readline().rstrip('\n\r'))
+#dumpfile = inputData.readline().rstrip('\n\r')
+#outModuloPrinter = np.uint32(inputData.readline().rstrip('\n\r'))
+#outputPrinter = inputData.readline().rstrip('\n\r')
+#outModuloForces = np.uint32(inputData.readline().rstrip('\n\r'))
+#outputForces = inputData.readline().rstrip('\n\r')
+#outModuloEnergies = np.uint32(inputData.readline().rstrip('\n\r'))
+#outputEnergies = inputData.readline().rstrip('\n\r')
+#visco = np.float64(inputData.readline().rstrip('\n\r'))
+#vortProj = np.uint32(inputData.readline().rstrip('\n\r'))
+#methodAdv = inputData.readline().rstrip('\n\r')
+#methodStretch = inputData.readline().rstrip('\n\r')
+#multires = np.uint32(inputData.readline().rstrip('\n\r'))
+#resFilterX = np.uint32(inputData.readline().rstrip('\n\r'))
+#resFilterY = np.uint32(inputData.readline().rstrip('\n\r'))
+#resFilterZ = np.uint32(inputData.readline().rstrip('\n\r'))
+
+#inputData.close()
# Parameters
dim = 3
-nbElem = [resX, resY, resZ]
-nbElemFilter = [resFilterX, resFilterY, resFilterZ]
+nbElem = [data.resX, data.resY, data.resZ]
+nbElemFilter = [data.resFilterX, data.resFilterY, data.resFilterZ]
t0 = time.time()
## Domain
-box = pp.Box(dim, length=[2.0 * pi, 2.0 * pi, 2.0 * pi], origin=[Ox, Oy, Oz])
+box = pp.Box(dim, length=[2.0 * pi, 2.0 * pi, 2.0 * pi],
+ origin=[data.Ox, data.Oy, data.Oz])
## Fields declaration
def computeVel(x, y, z):
@@ -145,28 +149,28 @@ vorti = pp.Field(domain=box, formula=computeVort,
advec = Advection(velo, vorti,
resolutions={velo: nbElem,
vorti: nbElem},
- method=methodAdv
+ method=data.methodAdv
)
stretch = Stretching(velo, vorti,
resolutions={velo: nbElem,
vorti: nbElem},
- method=methodStretch,
- propertyOp='divConservation'
+ method=data.methodStretch,
+ propertyOp=data.propStretch,
)
diffusion = Diffusion(vorti,
resolution=nbElem,
method='fftw',
- viscosity=visco
+ viscosity=data.visco
)
poisson = Poisson(velo, vorti,
resolutions={velo: nbElem,
vorti: nbElem},
method='fftw',
- projection=vortProj,
- multires=multires,
+ projection=data.vortProj,
+ multires=data.multires,
filterSize=box.length / \
(np.asarray(nbElemFilter,
dtype=PARMES_REAL) - 1)
@@ -176,37 +180,39 @@ poisson = Poisson(velo, vorti,
energy = Energy_enstrophy(velo, vorti,
resolutions={velo: nbElem,
vorti: nbElem},
- viscosity=visco,
- frequency=outModuloEnergies,
- prefix=outputEnergies)
+ viscosity=data.visco,
+ frequency=data.outModuloEnergies,
+ prefix=data.outputEnergies)
reproj = Reprojection_criterion(velo, vorti,
resolutions={velo: nbElem,
vorti: nbElem},
method='FD_order4',
frequency=1,
- prefix='./res/Reproj.dat')
+ prefix='./Reproj.dat')
printer = Printer(fields=[vorti, velo],
- frequency=outModuloPrinter,
- prefix=outputPrinter)
+ frequency=data.outModuloPrinter,
+ prefix=data.outputPrinter,
+ ext='.vtk')
distrAdvStr = Redistribute([vorti, velo], advec, stretch)
distrStrPoiss = Redistribute([vorti, velo], stretch, poisson)
## Definition of simulation parameters
-simu = Simulation(tinit=0.0, tend=finalTime, timeStep=dt, iterMax=1000000)
+simu = Simulation(tinit=0.0, tend=data.finalTime, timeStep=data.timeStep, iterMax=1000000)
## Define the problem to solve
#pb = NSProblem([diffusion, poisson],
pb = NSProblem(operators=[advec, distrAdvStr, stretch, distrStrPoiss, diffusion, poisson],
- simulation=simu, monitors=[printer, energy, reproj], dumpFreq=int(dumpfreq), name=dumpfile)
+ simulation=simu, monitors=[printer, energy, reproj],
+ dumpFreq=int(data.dumpfreq), name=data.dumpfile)
## Setting solver to Problem
pb.setUp()
## Restarting Problem from dumped field values
-if restart :
- pb.restart(dumpfile)
+if data.restart :
+ pb.restart(data.dumpfile)
t1 = time.time()
diff --git a/HySoP/hysop/numerics/differentialOperator.py b/HySoP/hysop/numerics/differentialOperator.py
index 6dce5f092..7e377277e 100755
--- a/HySoP/hysop/numerics/differentialOperator.py
+++ b/HySoP/hysop/numerics/differentialOperator.py
@@ -3,9 +3,11 @@
@file differentialOperator.py
Discrete stretching representation
"""
+import numpy as np
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_S import SynchronizeS
class DifferentialOperator(NumMethod):
@@ -31,7 +33,7 @@ class DifferentialOperator(NumMethod):
self.choice = choice
self.method = method
self.compute_time = 0.
- self.topology = self.field1.topology
+ self.topology = self.field2.topology
self.meshSize = self.topology.mesh.space_step
# def setUp(self):
@@ -71,8 +73,8 @@ class DifferentialOperator(NumMethod):
# time.sleep(2)
# print 'Aligne 3', field1[linenb,linenb,:]
# time.sleep(2)
-# OpSynchronize = Synchronize(self.topology)
-# OpSynchronize.apply(self.field2, self.field1)
+ OpSynchronize = Synchronize(self.topology)
+ OpSynchronize.apply(self.field2, self.field1)
synchroOp.apply()
if self.method.find('FD_order2') >= 0:
@@ -185,7 +187,8 @@ class DifferentialOperator(NumMethod):
self.field2[0][ind0a:ind0b, ind1a:ind1b, ind2+2]
) / (12. * self.meshSize[2])
- tmp1 = np.array([temp1 + temp2 + temp3])
+# tmp1 = np.array([temp1 + temp2 + temp3])
+ tmp1 = temp1 + temp2 + temp3
# Y components of temp and result
temp1 = self.field2[1][...] * 0.
@@ -224,7 +227,8 @@ class DifferentialOperator(NumMethod):
self.field2[1][ind0a:ind0b, ind1a:ind1b, ind2+2]
) / (12. * self.meshSize[2])
- tmp2 = np.array([temp1 + temp2 + temp3])
+# tmp2 = np.array([temp1 + temp2 + temp3])
+ tmp2 = temp1 + temp2 + temp3
# Z components of temp and result
temp1 = self.field2[2][...] * 0.
@@ -263,10 +267,12 @@ class DifferentialOperator(NumMethod):
self.field2[2][ind0a:ind0b, ind1a:ind1b, ind2+2]
) / (12. * self.meshSize[2])
- tmp3 = np.array([temp1 + temp2 + temp3])
+# tmp3 = np.array([temp1 + temp2 + temp3])
+ tmp3 = temp1 + temp2 + temp3
result = np.concatenate((
np.array(tmp1), np.array(tmp2), np.array(tmp3)))
- return result
+# return result
+ return tmp1, tmp2, tmp3
elif self.choice.find('gradV') >= 0:
@@ -292,8 +298,6 @@ 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.
@@ -304,8 +308,6 @@ 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.
@@ -318,7 +320,6 @@ 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(
@@ -331,8 +332,6 @@ 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.
@@ -343,8 +342,6 @@ 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.
@@ -357,7 +354,6 @@ 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(
@@ -370,8 +366,6 @@ 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.
@@ -382,8 +376,6 @@ 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.
@@ -396,7 +388,6 @@ class DifferentialOperator(NumMethod):
maxstr3 = np.max(abs(temp7)+abs(temp8)+abs(temp9))
maxadv3 = np.max(abs(temp9))
- max9 = np.max(abs(temp9))
# grad(V) result
result = np.concatenate((
@@ -407,21 +398,11 @@ class DifferentialOperator(NumMethod):
# 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
diff --git a/HySoP/hysop/numerics/fct2op.py b/HySoP/hysop/numerics/fct2op.py
index 9cbf48054..ee3ea47d0 100644
--- a/HySoP/hysop/numerics/fct2op.py
+++ b/HySoP/hysop/numerics/fct2op.py
@@ -82,7 +82,7 @@ class Fct2Op(object):
if self.choice == 'divWU': # field2 = velocity field
result = DifferentialOperator(
y, self.field2, self.method, self.choice).__call__()
- return result
+ return result[0], result[1], result[2]
if self.choice == 'hessianU': # field2 = velocity field
gradU, maxgersh = DifferentialOperator(
diff --git a/HySoP/hysop/operator/discrete/multiphase.py b/HySoP/hysop/operator/discrete/multiphase.py
index 4ea3323e9..8e90160a6 100644
--- a/HySoP/hysop/operator/discrete/multiphase.py
+++ b/HySoP/hysop/operator/discrete/multiphase.py
@@ -5,30 +5,32 @@ Discrete MultiPhase Rot Grad P
"""
from discrete import DiscreteOperator
from parmepy.numerics.differentialOperator import DifferentialOperator
+from parmepy.numerics.fct2op import Fct2Op
from parmepy.constants import np, debug
from parmepy.tools.timers import timed_function
-class Pressure_d(DiscreteOperator):
+class Baroclinic_d(DiscreteOperator):
"""
"""
@debug
- def __init__(self, velocity, vorticity, density, method=None):
+ def __init__(self, velocity, vorticity, density, viscosity, method=None):
"""
Constructor.
Create the old velocity field to compute the -Dp/rho
in N.S equation
@param operator.
"""
- DiscreteOperator.__init__(self, [velocity, vorticity, density], method,
- name="Pressure")
+ DiscreteOperator.__init__(self, [velocity, vorticity, density, viscosity], method,
+ name="Baroclinic_d")
- self.velocity_old = np.copy(velocity)
- self.velocity_new = velocity
+ self.velocity = velocity
+ self.velocity_old = np.zeros_like(self.velocity.data)
self.vorticity = vorticity
self.density = density
- self.input = [self.velocity_new, self.vorticity, self.density]
+ self.viscosity = viscosity
+ self.input = [self.velocity, self.vorticity, self.density, self.viscosity]
self.output = [self.vorticity]
@debug
@@ -37,34 +39,60 @@ class Pressure_d(DiscreteOperator):
if simulation is None:
raise ValueError("Missing simulation value for computation.")
+ t = simulation.time
dt = simulation.timeStep
- result = (self.velocity_new - self.velocity_old) / dt
- gradU = DifferentialOperator(self.velocity_new, self.velocity_new,
- choice='gradV',
- topology=self.topology).apply()
-
- result = result + np.dot(self.velocity_new, gradU)
- ## result = - gradP/rho
- result = result - self.viscosity *\
- DifferentialOperator(self.velocity_new,
- self.velocity_new,
- choice='laplacianV',
- topology=self.topology).apply()
- ## result = -(gradP/rho) x (gradRho/rho)
- result = np.cross(result,
- DifferentialOperator(self.density,
+ iCompute = self.velocity.topology.mesh.iCompute
+
+ if t > dt :
+ result = np.zeros_like(self.velocity.data)
+
+ # result = du/dt
+ for d in xrange(self.velocity.dimension):
+ result[d][iCompute] = (self.velocity[d][iCompute] - self.velocity_old[d][iCompute]) / dt
+
+ # result = result + (u. grad)u
+ UgradU = DifferentialOperator(self.velocity,
+ self.velocity,
+ method = self.method,
+ choice='divWU').__call__()
+
+ for d in xrange(self.velocity.dimension):
+ result[d][iCompute] = result[d][iCompute] + UgradU[d][iCompute]
+
+ ## result = result -nu*\Laplacian u = - gradP/rho
+ viscousTerm = DifferentialOperator(self.velocity,
+ self.velocity,
+ method = self.method,
+ choice='laplacianV').__call__()
+ for d in xrange(self.velocity.dimension):
+ viscousTerm[d][iCompute] = self.viscosity[iCompute] * viscousTerm[d][iCompute]
+ result[d][iCompute] = result[d][iCompute] - viscousTerm[d][iCompute]
+
+ ## barotrop = -(gradP/rho) x (gradRho/rho)
+ gradRho_rho = DifferentialOperator(self.density,
self.density,
- choice='gradS',
- topology=self.topology
- ).apply() / self.density)
+ method = self.method,
+ choice='gradS').__call__()
+ for d in xrange(self.velocity.dimension):
+ gradRho_rho[d][iCompute] = gradRho_rho[d][iCompute] / self.density[iCompute]
+ barotrop = np.zeros_like(result)
+ barotrop[0][iCompute] = result[1][iCompute] * gradRho_rho[2][iCompute] - \
+ result[2][iCompute] * gradRho_rho[1][iCompute]
+ barotrop[1][iCompute] = result[2][iCompute] * gradRho_rho[0][iCompute] - \
+ result[0][iCompute] * gradRho_rho[2][iCompute]
+ barotrop[2][iCompute] = result[0][iCompute] * gradRho_rho[1][iCompute] - \
+ result[1][iCompute] * gradRho_rho[0][iCompute]
+
+ ## vorti(n+1) = vorti(n) + dt * barotrop
+ for d in xrange(self.vorticity.dimension):
+ self.vorticity[d][iCompute] = self.vorticity[d][iCompute] + dt * barotrop[d][iCompute]
- ## vorti(n+1) = vorti(n) + dt * result
- self.vorticity.data = self.vorticity.data + dt * result
- self.velocity_old = np.copy(self.velocity_new.data)
+ for d in xrange(self.velocity.dimension):
+ self.velocity_old[d][iCompute] = self.velocity.data[d][iCompute]
if __name__ == "__main__":
print __doc__
- print "- Provided class : Pressure_d"
- print Pressure_d.__doc__
+ print "- Provided class : Baroclinic_d"
+ print Baroclinic_d.__doc__
diff --git a/HySoP/hysop/operator/discrete/poisson_fft.py b/HySoP/hysop/operator/discrete/poisson_fft.py
index f0616c3ef..f92567cf5 100644
--- a/HySoP/hysop/operator/discrete/poisson_fft.py
+++ b/HySoP/hysop/operator/discrete/poisson_fft.py
@@ -67,7 +67,7 @@ class PoissonFFT(DiscreteOperator):
elif self.case is '3D':
#=== SOLENOIDAL PROJECTION OF VORTICITY FIELD ===
- if (self.projection and ite % 10 == 0):
+ if (self.projection and ite % 1 == 0):
self.vorticity.data[0], self.vorticity.data[1], \
self.vorticity.data[2] = \
fftw2py.projection_om_3d(self.vorticity.data[0],
diff --git a/HySoP/hysop/operator/discrete/stretching.py b/HySoP/hysop/operator/discrete/stretching.py
index 46ab540e3..48c6fd9b8 100755
--- a/HySoP/hysop/operator/discrete/stretching.py
+++ b/HySoP/hysop/operator/discrete/stretching.py
@@ -86,7 +86,7 @@ class Stretching_d(DiscreteOperator):
use default divConservation'
self.propertyOp = 'divConservation'
- if self.propertyOp == 'massConservation':
+ if self.propertyOp.find('massConservation') >= 0:
## Determination of Stretching time step (stability condition)
self.ndt = 1
self.dtstretch = dt
--
GitLab