diff --git a/Examples/FlowAroundSphere.py b/Examples/FlowAroundSphere.py
index f1686910ee0f9d3aae79369a59a6f43dc04ac6c7..0a9ffa1ad9b601dda7e02ae0e8e5b6898aac9b83 100644
--- a/Examples/FlowAroundSphere.py
+++ b/Examples/FlowAroundSphere.py
@@ -49,10 +49,10 @@ VISCOSITY = 1. / 300.
# ======= Domain =======
dim = 3
-#Nx = 513
-#Ny = Nz = 257
-Nx = 257
-Ny = Nz = 129
+Nx = 513
+Ny = Nz = 257
+#Nx = 257
+#Ny = Nz = 129
g = 2
boxlength = npw.asrealarray([10.24, 5.12, 5.12])
boxorigin = npw.asrealarray([-2.0, -2.56, -2.56])
diff --git a/Examples/FlowAroundSphere_pressure.py b/Examples/FlowAroundSphere_pressure.py
new file mode 100644
index 0000000000000000000000000000000000000000..31b3dc80108f57a30eca29681fcecb2b535a14ab
--- /dev/null
+++ b/Examples/FlowAroundSphere_pressure.py
@@ -0,0 +1,434 @@
+#!/usr/bin/python
+
+"""
+Taylor Green 3D : see paper van Rees 2011.
+
+All parameters are set and defined in python module dataTG.
+
+"""
+
+from hysop import Box
+from hysop.f2py import fftw2py
+import numpy as np
+from hysop.fields.continuous import Field
+from hysop.fields.variable_parameter import VariableParameter
+from hysop.mpi.topology import Cartesian
+from hysop.operator.advection import Advection
+from hysop.operator.stretching import Stretching
+from hysop.operator.absorption_BC import AbsorptionBC
+from hysop.operator.poisson import Poisson
+from hysop.operator.diffusion import Diffusion
+from hysop.operator.adapt_timestep import AdaptTimeStep
+from hysop.operator.differential import DivAdvection
+from hysop.operator.redistribute_intra import RedistributeIntra
+from hysop.operator.hdf_io import HDF_Writer
+from hysop.operator.energy_enstrophy import EnergyEnstrophy
+from hysop.problem.simulation import Simulation
+from hysop.methods_keys import Scales, TimeIntegrator, Interpolation,\
+ Remesh, Support, Splitting, dtCrit, SpaceDiscretisation, \
+ GhostUpdate, Formulation
+from hysop.numerics.integrators.runge_kutta2 import RK2 as RK2
+from hysop.numerics.integrators.runge_kutta3 import RK3 as RK3
+from hysop.numerics.integrators.runge_kutta4 import RK4 as RK4
+from hysop.numerics.finite_differences import FD_C_4, FD_C_2
+from hysop.numerics.interpolation import Linear
+from hysop.numerics.remeshing import L6_4 as rmsh
+import hysop.tools.io_utils as io
+import hysop.tools.numpywrappers as npw
+from hysop.mpi import main_rank, MPI
+from hysop.tools.parameters import Discretization, IOParams
+
+print " ========= Start Navier-Stokes 3D (Taylor Green benchmark) ========="
+
+# ====== pi constant and trigonometric functions ======
+pi = np.pi
+cos = np.cos
+sin = np.sin
+
+# ====== Flow constants =======
+uinf = 1.0
+VISCOSITY = 1. / 300.
+
+# ======= Domain =======
+dim = 3
+#Nx = 513
+#Ny = Nz = 257
+Nx = 257
+Ny = Nz = 129
+g = 2
+boxlength = npw.asrealarray([10.24, 5.12, 5.12])
+boxorigin = npw.asrealarray([-2.0, -2.56, -2.56])
+box = Box(length=boxlength, origin=boxorigin)
+
+# A global discretization with ghost points
+d3dg = Discretization([Nx, Ny, Nz], [g, g, g])
+# A global discretization, without ghost points
+d3d = Discretization([Nx, Ny, Nz])
+
+# ====== Sphere inside the domain ======
+RADIUS = 0.5
+pos = [0., 0., 0.]
+from hysop.domain.subsets.sphere import Sphere, HemiSphere
+#sphere = Sphere(origin=pos, radius=RADIUS, parent=box)
+sphere = HemiSphere(origin=pos, radius=RADIUS, parent=box)
+
+
+# ======= Function to compute initial velocity =======
+def computeVel(res, x, y, z, t):
+ res[0][...] = uinf
+ res[1][...] = 0.
+ res[2][...] = 0.
+ return res
+
+
+# ======= Function to compute initial vorticity =======
+def computeVort(res, x, y, z, t):
+ res[0][...] = 0.
+ res[1][...] = 0.
+ res[2][...] = 0.
+ return res
+
+# ======= Function to compute initial pressure =======
+def computePressure(res, x, y, z, t):
+ res[0][...] = 0.
+ return res
+
+# ====== Time-dependant required-flowrate (Variable Parameter) ======
+def computeFlowrate(simu):
+ # === Time-dependant flow rate ===
+ t = simu.tk
+ Tstart = 3.0
+ flowrate = np.zeros(3)
+ flowrate[0] = uinf * box.length[1] * box.length[2]
+ if t >= Tstart and t <= Tstart + 1.0:
+ flowrate[1] = sin(pi * (t - Tstart)) * \
+ box.length[1] * box.length[2]
+ # === Constant flow rate ===
+ # flowrate = np.zeros(3)
+ # flowrate[0] = uinf * box.length[1] * box.length[2]
+ return flowrate
+
+
+# ======= Fields =======
+velo = Field(domain=box, formula=computeVel,
+ name='Velocity', is_vector=True)
+vorti = Field(domain=box, formula=computeVort,
+ name='Vorticity', is_vector=True)
+pressure = Field(domain=box, formula=computePressure,
+ name='Pressure')
+
+# ========= Simulation setup =========
+simu = Simulation(tinit=0.0, tend=80.0, timeStep=0.0125, iterMax=10000000)
+
+
+# Adaptative timestep method : dt = min(values(dtCrit))
+# where dtCrit is a list of criterions on which the computation
+# of the adaptative time step is based
+# ex : dtCrit = ['gradU', 'cfl', 'stretch'], means :
+# dt = min (dtAdv, dtCfl, dtStretch), where dtAdv is equal to LCFL / |gradU|
+# For dtAdv, the possible choices are the following:
+# 'vort' (infinite norm of vorticity) : dtAdv = LCFL / |vort|
+# 'gradU' (infinite norm of velocity gradient), dtAdv = LCFL / |gradU|
+# 'deform' (infinite norm of deformation tensor),
+# dtAdv = LCFL / (0.5(gradU + gradU^T))
+op = {}
+iop = IOParams("time_step")
+# Default topology (i.e. 3D, with ghosts)
+topo_with_ghosts = box.create_topology(d3dg)
+
+
+op['dtAdapt'] = AdaptTimeStep(velo, vorti, simulation=simu,
+ discretization=topo_with_ghosts,
+ method={TimeIntegrator: RK3,
+ SpaceDiscretisation: FD_C_4,
+ dtCrit: ['gradU', 'stretch']},
+ io_params=iop,
+ lcfl=0.125,
+ cfl=0.5)
+
+op['advection'] = Advection(velo, vorti,
+ discretization=d3d,
+ method={Scales: 'p_M6',
+ Splitting: 'classic'}
+ )
+
+op['stretching'] = Stretching(velo, vorti,
+ discretization=topo_with_ghosts)
+
+op['diffusion'] = Diffusion(viscosity=VISCOSITY, vorticity=vorti,
+ discretization=d3d)
+
+rate = VariableParameter(formula=computeFlowrate)
+op['poisson'] = Poisson(velo, vorti, discretization=d3d, flowrate=rate)
+
+op['rhsPrsPoisson'] = DivAdvection(velo, pressure,
+ discretization=topo_with_ghosts)
+
+op['poissonPressure'] = Poisson(pressure, pressure, discretization=d3d,
+ method={SpaceDiscretisation: 'fftw',
+ GhostUpdate: True,
+ Formulation: 'pressure'})
+
+# ===== Discretization of computational operators ======
+for ope in op.values():
+ ope.discretize()
+
+topofft = op['poisson'].discreteFields[vorti].topology
+topoadvec = op['advection'].discreteFields[vorti].topology
+
+# ===== Smooth vorticity absorption at the outlet =====
+op['vort_absorption'] = AbsorptionBC(velo, vorti, discretization=topofft,
+ req_flowrate=rate,
+ x_coords_absorp=[7.24, 8.24])
+# x_coords_absorp=[1.56, 2.56])
+op['vort_absorption'].discretize()
+
+# ===== Penalization of the vorticity on a sphere inside the domain =====
+from hysop.operator.penalize_vorticity import PenalizeVorticity
+op['penalVort'] = PenalizeVorticity(velocity=velo, vorticity=vorti,
+ discretization=topo_with_ghosts,
+ obstacles=[sphere], coeff=1e8,
+ method={SpaceDiscretisation: FD_C_4})
+op['penalVort'].discretize()
+
+
+# ==== Operators to map data between the different computational operators ===
+# (i.e. between topologies)
+distr = {}
+distr['fft2str'] = RedistributeIntra(source=op['poisson'],
+ target=op['stretching'],
+ variables=[velo, vorti])
+distr['str2fft'] = RedistributeIntra(source=op['stretching'],
+ target=op['poisson'],
+ variables=[velo, vorti])
+distr['fft2advec'] = RedistributeIntra(source=op['poisson'],
+ target=op['advection'],
+ variables=[velo, vorti])
+distr['advec2fft'] = RedistributeIntra(source=op['advection'],
+ target=op['poisson'],
+ variables=[velo, vorti])
+distr['str2prs'] = RedistributeIntra(source=op['rhsPrsPoisson'],
+ target=op['poissonPressure'],
+ variables=[pressure])
+# ========= Monitoring operators =========
+monitors = {}
+iop = IOParams('fields', frequency=100)
+monitors['writer'] = HDF_Writer(variables={velo: topofft, vorti: topofft},
+ io_params=iop)
+
+io_ener = IOParams('energy_enstrophy')
+monitors['energy'] = EnergyEnstrophy(velo, vorti, discretization=topofft,
+ io_params=io_ener, is_normalized=False)
+
+from hysop.domain.subsets.control_box import ControlBox
+from hysop.operator.drag_and_lift import MomentumForces, NocaForces
+ref_step = topo_with_ghosts.mesh.space_step
+cbpos = npw.zeros(dim)
+cblength = npw.zeros(dim)
+cbpos[...] = boxorigin[...]
+cbpos += 15 * ref_step
+cblength[...] = boxlength[...]
+cblength -= 30 * ref_step
+cb = ControlBox(parent=box, origin=cbpos, length=cblength)
+coeffForce = 1. / (0.5 * uinf ** 2 * pi * RADIUS ** 2)
+
+io_forces=IOParams('drag_and_lift_NocaII')
+#monitors['forcesNoca'] = NocaForces(velo, vorti,
+# discretization=topo_with_ghosts,
+# nu=VISCOSITY,
+# volume_of_control=cb,
+# normalization=coeffForce,
+# obstacles=[sphere],
+# io_params=io_forces)
+
+io_forcesPenal=IOParams('drag_and_lift_Mom')
+monitors['forcesMom'] = MomentumForces(velocity=velo,
+ discretization=topo_with_ghosts,
+ normalization=coeffForce,
+ obstacles=[sphere],
+ penalisation_coeff=[1e8],
+ io_params=io_forcesPenal)
+
+#io_forcesPenal=IOParams('drag_and_lift_penal')
+#monitors['forcesPenal'] = DragAndLiftPenal(velo, vorti, coeffForce,
+# discretization=topofft,
+# obstacles=[sphere], factor=[1e8],
+# io_params=io_forcesPenal)
+
+step_dir = ref_step[0]
+io_sliceXY = IOParams('sliceXY', frequency=5)
+thickSliceXY = ControlBox(parent=box, origin=[-2.0, -2.56, -2.0 * step_dir],
+ length=[10.24- step_dir, 5.12- step_dir, 4.0 * step_dir])
+#thickSliceXY = ControlBox(parent=box, origin=[-2.56, -2.56, -2.0 * step_dir],
+# length=[5.12 - step_dir, 5.12 - step_dir, 4.0 * step_dir])
+monitors['writerSliceXY'] = HDF_Writer(variables={velo: topofft, vorti: topofft},
+ io_params=io_sliceXY, subset=thickSliceXY,
+ xmfalways=True)
+
+#io_sliceXZ = IOParams('sliceXZ', frequency=2000)
+#thickSliceXZ = ControlBox(box, origin=[-2.0, -2.0 * step_dir, -2.56],
+# lengths=[10.24, 4.0 * step_dir, 5.12])
+#monitors['writerSliceXZ'] = HDF_Writer(variables={velo: topofft, vorti: topofft},
+# io_params=io_sliceXZ, subset=thickSliceXZ,
+# xmfalways=True)
+
+io_slicePrs = IOParams('slicePrs', frequency=5)
+thickSlicePrs = ControlBox(parent=box, origin=[-2.0, -2.56, -2.0 * step_dir],
+ length=[10.24- step_dir, 5.12- step_dir, 4.0 * step_dir])
+monitors['writerSlicePrs'] = HDF_Writer(variables={pressure: topofft},
+ io_params=io_slicePrs, subset=thickSlicePrs,
+ xmfalways=True)
+
+#io_subBox = IOParams('subBox', frequency=2000)
+#subBox = ControlBox(box, origin=[-0.7, -2.0, -2.0], lengths=[8.0, 4.0, 4.0])
+#monitors['writerSubBox'] = HDF_Writer(variables={velo: topofft, vorti: topofft},
+# io_params=io_subBox, subset=subBox,
+# xmfalways=True)
+
+# ========= Setup for all declared operators/monitors =========
+time_setup = MPI.Wtime()
+for ope in op.values():
+ ope.setup()
+for ope in distr.values():
+ ope.setup()
+for monit in monitors.values():
+ monit.discretize()
+for monit in monitors.values():
+ monit.setup()
+
+print '[', main_rank, '] total time for setup:', MPI.Wtime() - time_setup
+
+# ========= Fields initialization =========
+# - initialize velo + vort on topostr
+# - penalize vorticity
+# - redistribute topostr --> topofft
+
+time_init = MPI.Wtime()
+ind = sphere.discretize(topofft)
+def initFields():
+ velo.initialize(topo=topo_with_ghosts)
+ vorti.initialize(topo=topo_with_ghosts)
+ pressure.initialize(topo=topo_with_ghosts)
+ op['penalVort'].apply(simu)
+ distr['str2fft'].apply(simu)
+ distr['str2fft'].wait()
+
+initFields()
+print '[', main_rank, '] total time for init :', MPI.Wtime() - time_init
+
+fullseq = []
+
+def run(sequence):
+ op['vort_absorption'].apply(simu)
+ op['poisson'].apply(simu) # Poisson + correction
+ monitors['forcesMom'].apply(simu) # Forces Heloise
+ distr['fft2str'].apply(simu)
+ distr['fft2str'].wait()
+ op['penalVort'].apply(simu) # Vorticity penalization
+# distr['str2fft'].apply(simu)
+# distr['str2fft'].wait()
+# op['poisson'].apply(simu)
+# distr['fft2str'].apply(simu)
+# distr['fft2str'].wait()
+ op['stretching'].apply(simu) # Stretching
+ op['rhsPrsPoisson'].apply(simu) # RHS computation in Pressure Poisson eq
+# monitors['forcesNoca'].apply(simu) # Forces Noca
+ distr['str2fft'].apply(simu)
+ distr['str2fft'].wait()
+ distr['str2prs'].apply(simu)
+ distr['str2prs'].wait()
+ op['diffusion'].apply(simu) # Diffusion
+ op['poissonPressure'].apply(simu) # Pressure Poisson
+ distr['fft2advec'].apply(simu)
+ distr['fft2advec'].wait()
+ op['advection'].apply(simu) # Advection (scales)
+ distr['advec2fft'].apply(simu)
+ distr['advec2fft'].wait()
+ monitors['writerSliceXY'].apply(simu)
+# monitors['writerSliceXZ'].apply(simu)
+ monitors['writerSlicePrs'].apply(simu)
+# monitors['writerSubBox'].apply(simu)
+ monitors['energy'].apply(simu) # Energy/enstrophy
+ distr['fft2str'].apply(simu)
+ distr['fft2str'].wait()
+ op['dtAdapt'].apply(simu) # Update timestep
+ op['dtAdapt'].wait()
+
+# ==== Serialize some data of the problem to a "restart" file ====
+# def dump(filename):
+# """
+# Serialize some data of the problem to file
+# (only data required for a proper restart, namely fields in self.input
+# and simulation).
+# @param filename : prefix for output file. Real name = filename_rk_N,
+# N being current process number. If None use default value from problem
+# parameters (self.filename)
+# """
+# if filename is not None:
+# filename = filename
+# filedump = filename + '_rk_' + str(main_rank)
+# db = parmesPickle(filedump, mode='store')
+# db.dump(simu, 'simulation')
+# velo.dump(filename, mode='append')
+# vorti.dump(filename, mode='append')
+
+# ## ====== Load some data of the problem from a "restart" file ======
+# def restart(filename):
+# """
+# Load serialized data to restart from a previous state.
+# self.input variables and simulation are loaded.
+# @param filename : prefix for downloaded file.
+# Real name = filename_rk_N, N being current process number.
+# If None use default value from problem
+# parameters (self.filename)
+# """
+# if filename is not None:
+# filename = filename
+# filedump = filename + '_rk_' + str(main_rank)
+# db = parmesPickle(filedump, mode='load')
+# simu = db.load('simulation')[0]
+# simu.start = simu.time - simu.timeStep
+# ite = simu.currentIteration
+# simu.initialize()
+# simu.currentIteration = ite
+# print 'simu', simu
+# print ("load ...", filename)
+# velo.load(filename)
+# vorti.load(filename)
+# return simu
+
+seq = fullseq
+
+simu.initialize()
+#doDump = False
+#doRestart = False
+#dumpFreq = 5000
+#io_default={"filename":'restart'}
+#dump_filename = io.Writer(params=io_default).filename
+#===== Restart (if needed) =====
+# if doRestart:
+# simu = restart(dump_filename)
+# # Set up for monitors and redistribute
+# for ope in distr.values():
+# ope.setUp()
+# for monit in monitors.values():
+# monit.setUp()
+
+# ======= Time loop =======
+time_run = MPI.Wtime()
+while not simu.isOver:
+ if topofft.rank == 0:
+ simu.printState()
+ run(seq)
+ simu.advance()
+# # testdump = simu.currentIteration % dumpFreq is 0
+# # if doDump and testdump:
+# # dump(dump_filename)
+print '[', main_rank, '] total time for run :', MPI.Wtime() - time_run
+
+# ======= Finalize =======
+fftw2py.clean_fftw_solver(box.dimension)
+for ope in distr.values():
+ ope.finalize()
+for monit in monitors.values():
+ monit.finalize()
diff --git a/HySoP/hysop/numerics/differential_operations.py b/HySoP/hysop/numerics/differential_operations.py
index 5dab84270755d00a15dcb0978302d5d5583b942c..6466a62e002de87e20abeee00cc38fc320803442 100755
--- a/HySoP/hysop/numerics/differential_operations.py
+++ b/HySoP/hysop/numerics/differential_operations.py
@@ -509,34 +509,68 @@ class DivAdvection(DifferentialOperation):
# dimension of the fields
self._dim = topo.domain.dimension
+ @staticmethod
+ def getWorkLengths(nb_components=None, domain_dim=None, fd_method=None):
+ """
+ fd_method = FD_C_4 --> call fd_compute
+ fd_method = FD_C_4_CAA --> call fd.compute_and_add
+ """
+ if domain_dim == 1 or nb_components == 1:
+ return 1
+ elif domain_dim == 2 or nb_components == 2:
+ return 2
+ else:
+ return 3
+
def __call__(self, var1, result):
return self.fcall(var1, result)
def FDCentral(self, var1, result):
-
-
- work = npw.zeros_like(rhs)
-
- j1 = np.arange(self.velocity.nbComponents)
- rhs[...] = 0.0
-
- for i in xrange(self.velocity.nbComponents):
- j2 = np.where(j1 != i)[0]
- self._fd_scheme.compute(vdata[j2[0]], j2[1], work)
- self._fd_scheme.compute_and_mult(vdata[j2[1]], j2[0], work)
- rhs[...] += work[...]
- work[...] = 0.0
- rhs[...] *= -1.0
-
- for i in xrange(self.velocity.nbComponents):
- j2 = np.where(j1 != i)[0]
- self._fd_scheme.compute(vdata[j2[0]], j2[0], work)
- self._fd_scheme.compute_and_mult(vdata[j2[1]], j2[1], work)
- rhs[...] += work[...]
- work[...] = 0.0
-
- rhs[...] *= 2.0
-
+ self._work = []
+ memshape = result[0].shape
+ for i in xrange(self.getWorkLengths(len(var1), self._dim)):
+ self._work.append(npw.zeros(memshape))
+ assert len(self._work) == \
+ self.getWorkLengths(len(var1), self._dim)
+ result[0][...] = 0.0
+ nbc = len(var1)
+ for comp in xrange(nbc):
+ self.fd_scheme.compute(var1[comp], comp, self._work[comp])
+ result[0][...] += self._work[XDIR] * self._work[YDIR]
+ result[0][...] += self._work[XDIR] * self._work[ZDIR]
+ result[0][...] += self._work[YDIR] * self._work[ZDIR]
+ self.fd_scheme.compute(var1[XDIR], YDIR, self._work[0])
+ self.fd_scheme.compute(var1[YDIR], XDIR, self._work[1])
+ result[0][...] -= self._work[0] * self._work[1]
+ self.fd_scheme.compute(var1[XDIR], ZDIR, self._work[0])
+ self.fd_scheme.compute(var1[ZDIR], XDIR, self._work[1])
+ result[0][...] -= self._work[0] * self._work[1]
+ self.fd_scheme.compute(var1[YDIR], ZDIR, self._work[0])
+ self.fd_scheme.compute(var1[ZDIR], YDIR, self._work[1])
+ result[0][...] -= self._work[0] * self._work[1]
+
+ result[0][...] *= 2.0
+
+# # With only 1 work array
+# nbc = len(var1)
+# j1 = np.arange(nbc)
+# result[0][...] = 0.0
+
+# for i in xrange(nbc):
+# j2 = np.where(j1 != i)[0]
+# self._fd_scheme.compute(var1[j2[0]], j2[1], self._work[0])
+# self._fd_scheme.compute_and_mult(var1[j2[1]], j2[0], self._work[0])
+# result[0][...] += self._work[0][...]
+# self._work[0][...] = 0.0
+# result[0][...] *= -1.0
+
+# for i in xrange(nbc):
+# j2 = np.where(j1 != i)[0]
+# self._fd_scheme.compute(vdata[j2[0]], j2[0], self._work[0])
+# self._fd_scheme.compute_and_mult(var1[j2[1]], j2[1], self._work[0])
+# result[0][...] += self._work[0][...]
+# self._work[0][...] = 0.0
+
+# result[0][...] *= 2.0
-
return result
diff --git a/HySoP/hysop/operator/discrete/drag_and_lift.py b/HySoP/hysop/operator/discrete/drag_and_lift.py
index 0bf0d8fd1bead9d73138ca0e0b0be3ad48106e2e..2eea997d9a7676f6a53610a6dd42ed1eec51a872 100644
--- a/HySoP/hysop/operator/discrete/drag_and_lift.py
+++ b/HySoP/hysop/operator/discrete/drag_and_lift.py
@@ -233,12 +233,12 @@ class MomentumForces(Forces):
for i in xrange(len(self._indices)):
self._buff[...] = 0.0
icurrent = self._indices[i]
- print 'shape', np.shape(icurrent), 'r =', self._topology.rank, 'i =', i
+# print 'shape', np.shape(icurrent), 'r =', self._topology.rank, 'i =', i
current = self._buff[icurrent]
coeff = self._coeff[i] * dt / (1. + self._coeff[i] * dt)
current[...] = self.velocity.data[d][icurrent] * coeff
self.force[d] += npw.real_sum(current)
- print 'd=', d, 'force =', self.force[d]
+# print 'd=', d, 'force =', self.force[d]
# i = 0
# for ind in self._indices:
diff --git a/HySoP/hysop/operator/discrete/poisson_fft.py b/HySoP/hysop/operator/discrete/poisson_fft.py
index 1b1be1875639972255a3861528b969b807e51e33..ace81936fa90c282b10120c6d7f7a1b5100aaf9d 100644
--- a/HySoP/hysop/operator/discrete/poisson_fft.py
+++ b/HySoP/hysop/operator/discrete/poisson_fft.py
@@ -11,6 +11,7 @@ except ImportError:
from hysop.operator.discrete.discrete import DiscreteOperator
from hysop.operator.discrete.reprojection import Reprojection
+from hysop.methods_keys import Formulation, SpaceDiscretisation
from hysop.constants import debug
from hysop.tools.profiler import profile
@@ -23,7 +24,7 @@ class PoissonFFT(DiscreteOperator):
@debug
def __init__(self, output_field, input_field, projection=None,
- filterSize=None, correction=None, **kwds):
+ filterSize=None, correction=None, formulation=None, **kwds):
"""
Constructor.
@param[out] output_field : discretization of the solution field
@@ -56,6 +57,7 @@ class PoissonFFT(DiscreteOperator):
if self.dim == 2:
assert self.input_field.nb_components == 1
self.correction = correction
+ self.formulation = formulation
self.input = [self.input_field]
self.output = [self.output_field]
@@ -96,6 +98,8 @@ class PoissonFFT(DiscreteOperator):
else:
if multires:
self._solve = self._solve3D_multires
+ elif self.formulation is not None:
+ self._solve = self._solve_3d_scalar_fd
else:
self._solve = self._solve3D
else:
diff --git a/HySoP/hysop/operator/poisson.py b/HySoP/hysop/operator/poisson.py
index 9a237ae88736481da444184dd3fc25433865f4be..73d29b5351b5080641d2a63d42d161db72b59b80 100644
--- a/HySoP/hysop/operator/poisson.py
+++ b/HySoP/hysop/operator/poisson.py
@@ -9,7 +9,7 @@ from hysop.operator.discrete.poisson_fft import PoissonFFT
from hysop.constants import debug
from hysop.operator.velocity_correction import VelocityCorrection
from hysop.operator.reprojection import Reprojection
-from hysop.methods_keys import SpaceDiscretisation
+from hysop.methods_keys import SpaceDiscretisation, Formulation
from hysop.operator.continuous import opsetup
import hysop.default_methods as default
@@ -65,6 +65,12 @@ class Poisson(Computational):
if self.method[SpaceDiscretisation] is not 'fftw':
raise AttributeError("Method not yet implemented.")
+ # Deterlination of the Poisson equation formulation :
+ # Velo Poisson eq or Pressure Poisson eq
+ self.formulation = None
+ if self.method[Formulation] is not 'velocity':
+ self.formulation = self.method[Formulation]
+
self.input = [self.input_field]
self.output = [self.output_field]
if flowrate is not None:
@@ -77,7 +83,7 @@ class Poisson(Computational):
self._config = kwds
if projection is not None:
- self.output_field.append(self.input_field)
+ self.output.append(self.input_field)
def discretize(self):
# Poisson solver based on fftw
@@ -125,6 +131,7 @@ class Poisson(Computational):
self.discreteFields[self.input_field],
correction=cd,
rwork=rwork, iwork=iwork,
- projection=projection_discr)
+ projection=projection_discr,
+ formulation=self.formulation)
self._is_uptodate = True
diff --git a/HySoP/hysop/operator/tests/test_poisson.py b/HySoP/hysop/operator/tests/test_poisson.py
index bdbdb1858e1bb264093fbd4cd0d0cb2edc274bd9..7676e87dbfb32d91466dae6b7a4b098ff98ea1e2 100755
--- a/HySoP/hysop/operator/tests/test_poisson.py
+++ b/HySoP/hysop/operator/tests/test_poisson.py
@@ -6,6 +6,9 @@ from hysop.operator.analytic import Analytic
from hysop.operator.reprojection import Reprojection
from hysop.problem.simulation import Simulation
from hysop.tools.parameters import Discretization
+from hysop.methods_keys import Scales, TimeIntegrator, Interpolation,\
+ Remesh, Support, Splitting, dtCrit, SpaceDiscretisation, \
+ GhostUpdate, Formulation
import numpy as np
import hysop.tools.numpywrappers as npw
import math
@@ -32,6 +35,14 @@ def computeVort(res, x, y, z, t):
res[2][...] = coeff * cos(x * cc[0]) * cos(y * cc[1]) * sin(z * cc[2])
return res
+def computePressure(res, x, y, z, t):
+ res[0][...] = -3.0 * sin(x * cc[0]) * cos(y * cc[1]) * cos(z * cc[2])
+ return res
+
+def computeRefPressure(res, x, y, z, t):
+ res[0][...] = sin(x * cc[0]) * cos(y * cc[1]) * cos(z * cc[2])
+ return res
+
# ref. field
def computeRef(res, x, y, z, t):
@@ -82,6 +93,36 @@ def computeRef2D(res, x, y, t):
return res
+def test_Poisson_Pressure_3D():
+ dom = pp.Box(length=LL)
+
+ # Fields
+ ref = pp.Field(domain=dom, name='Ref')
+ pressure = pp.Field(domain=dom, formula=computePressure, name='Pressure')
+
+ # Definition of the Poisson operator
+ poisson = Poisson(pressure, pressure, discretization=d3D,
+ method={SpaceDiscretisation: 'fftw',
+ GhostUpdate: True,
+ Formulation: 'pressure'})
+
+ poisson.discretize()
+ poisson.setup()
+ topo = poisson.discreteFields[pressure].topology
+ # Analytic operator to compute the reference field
+ refOp = Analytic(variables={ref: topo}, formula=computeRefPressure)
+ simu = Simulation(nbIter=10)
+ refOp.discretize()
+ refOp.setup()
+ pressure.initialize(topo=topo)
+ poisson.apply(simu)
+ refOp.apply(simu)
+ assert np.allclose(ref.norm(topo), pressure.norm(topo))
+ refD = ref.discretize(topo)
+ prsD = pressure.discretize(topo)
+ assert np.allclose(prsD[0], refD[0])
+ poisson.finalize()
+
def test_Poisson3D():
dom = pp.Box(length=LL)
@@ -260,6 +301,7 @@ def test_Poisson3D_projection_2():
# This may be useful to run mpi tests
if __name__ == "__main__":
+ test_Poisson_Pressure_3D()
test_Poisson3D()
test_Poisson2D()
test_Poisson3D_correction()