From d8bd1e28f5a877a2c1c07a8aeecf105c8fc188b5 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Chlo=C3=A9=20Mimeau?= <Chloe.Mimeau@imag.fr>
Date: Wed, 3 Dec 2014 14:02:10 +0100
Subject: [PATCH] update sphere Example

---
 Examples/FlowAroundSphere.py | 109 ++++++++++++++++++++---------------
 1 file changed, 63 insertions(+), 46 deletions(-)

diff --git a/Examples/FlowAroundSphere.py b/Examples/FlowAroundSphere.py
index 6e43288b5..eee8df58d 100644
--- a/Examples/FlowAroundSphere.py
+++ b/Examples/FlowAroundSphere.py
@@ -31,6 +31,7 @@ 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, IO_params
 
@@ -50,9 +51,9 @@ dim = 3
 Nx = 129
 Ny = Nz = 65
 g = 2
-boxlength = npw.realarray([10.24, 5.12, 5.12])
-boxorigin = npw.realarray([-2.0, -2.56, -2.56])
-box = Box(dim, length=boxlength, origin=boxorigin)
+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])
@@ -81,7 +82,6 @@ def computeVort(res, x, y, z, t):
     res[2][...] = 0.
     return res
 
-
 #  ====== Time-dependant required-flowrate (Variable Parameter) ======
 def computeFlowrate(simu):
     # === Time-dependant flow rate ===
@@ -149,6 +149,7 @@ op['diffusion'] = Diffusion(viscosity=VISCOSITY, vorticity=vorti,
 
 rate = VariableParameter(formula=computeFlowrate)
 op['poisson'] = Poisson(velo, vorti, discretization=d3d, flowrate=rate)
+#op['poisson'] = Poisson(velo, vorti, discretization=d3d)
 
 # ===== Discretization of computational operators ======
 for ope in op.values():
@@ -158,8 +159,10 @@ topofft = op['poisson'].discreteFields[vorti].topology
 topoadvec = op['advection'].discreteFields[vorti].topology
 
 # =====  Penalization of the vorticity on a sphere inside the domain =====
-from hysop.operator.penalize_vorticity import PenalizeVorticity
-op['penalVort'] = PenalizeVorticity(velocity=velo, vorticity=vorti,
+#from hysop.operator.penalize_vorticity import PenalizeVorticity
+from hysop.operator.penalization_and_curl import PenalizationAndCurl
+#op['penalVort'] = PenalizeVorticity(velocity=velo, vorticity=vorti,
+op['penalVort'] = PenalizationAndCurl(velocity=velo, vorticity=vorti,
                                     discretization=topo_with_ghosts,
                                     obstacles=[sphere], coeff=1e8,
                                     method={SpaceDiscretisation: FD_C_4})
@@ -173,7 +176,7 @@ distr['fft2str'] = RedistributeIntra(source=op['poisson'],
                                      target=op['stretching'],
                                      variables=[velo, vorti])
 distr['str2fft'] = RedistributeIntra(source=op['stretching'],
-                                     target=op['diffusion'],
+                                     target=op['poisson'],
                                      variables=[velo, vorti])
 distr['fft2advec'] = RedistributeIntra(source=op['poisson'],
                                        target=op['advection'],
@@ -182,15 +185,17 @@ distr['advec2fft'] = RedistributeIntra(source=op['advection'],
                                        target=op['poisson'],
                                        variables=[velo, vorti])
 # ========= Monitoring operators =========
-
+monitors = {}
 iop = IO_params('fields', frequency=100)
 monitors['writer'] = HDF_Writer(variables={velo: topofft, vorti: topofft},
                                 io_params=iop)
 
 io_ener = IO_params('energy_enstrophy')
 monitors['energy'] = EnergyEnstrophy(velo, vorti, discretization=topofft,
-                                     io_params=io_ener)
+                                     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)
@@ -198,40 +203,51 @@ cbpos[...] = boxorigin[...]
 cbpos +=  15 * ref_step
 cblength[...] = boxlength[...]
 cblength -= 30 * ref_step
-cb = ControlBox(box, cbpos, cblength)
+cb = ControlBox(parent=box, origin=cbpos, length=cblength)
 coeffForce = 1. / (0.5 * uinf ** 2 * pi * RADIUS ** 2)
 
-io_forces=IO_params('drag_and_lift')
-monitors['forces'] = DragAndLift(velo, vorti, VISCOSITY, coeffForce,
-                                 discretization=topo_with_ghosts, cb, 
-                                 obstacles=[sphere], io_params=io_forces)
-
-io_forcesPenal=IO_params('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 = IO_params('sliceXY', frequency=200)
-thickSliceXY = ControlBox(box, origin=[-2.0, -2.56, -2.0 * step_dir], 
-                          lengths=[10.24, 5.12, 4.0 * step_dir])
-monitors['writerSliceXY'] = HDF_Writer(variables={velo: topofft, vorti: topofft},
-                                      io_params=io_sliceXY, subset=thickSliceXY, 
-                                      xmfalways=True)
-
-io_sliceXZ = IO_params('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_subBox = IO_params('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)
+io_forces=IO_params('drag_and_lift_Noca')
+#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=IO_params('drag_and_lift_Mom')
+monitors['forcesMom'] = MomentumForces(velocity=velo, 
+                                       discretization=topo_with_ghosts,
+                                       normalization=coeffForce,
+                                       obstacles=[sphere], 
+                                       io_params=io_forcesPenal)
+
+#io_forcesPenal=IO_params('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 = IO_params('sliceXY', frequency=200)
+#thickSliceXY = ControlBox(box, origin=[-2.0, -2.56, -2.0 * step_dir], 
+#                          lengths=[10.24, 5.12, 4.0 * step_dir])
+#monitors['writerSliceXY'] = HDF_Writer(variables={velo: topofft, vorti: topofft},
+#                                      io_params=io_sliceXY, subset=thickSliceXY, 
+#                                      xmfalways=True)
+
+#io_sliceXZ = IO_params('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_subBox = IO_params('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)
 
 ## ========= Monitors discretization========= 
 for monit in monitors.values():
@@ -244,7 +260,9 @@ for ope in op.values():
 for ope in distr.values():
     ope.setup()
 for monit in monitors.values():
-    monit.setUp()
+    monit.discretize()
+for monit in monitors.values():
+    monit.setup()
 
 print '[', main_rank, '] total time for setup:', MPI.Wtime() - time_setup
 
@@ -268,14 +286,13 @@ print '[', main_rank, '] total time for init :', MPI.Wtime() - time_init
 fullseq = []
 
 def run(sequence):
-    op['poisson'].apply(simu)               # Poisson
-    op['correction'].apply(simu)            # Velo correction
-    monitors['forcesPenal'].apply(simu)     # Forces Heloise
+    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
     op['stretching'].apply(simu)            # Stretching
-    monitors['forces'].apply(simu)          # Forces Noca
+#    monitors['forcesNoca'].apply(simu)          # Forces Noca
     distr['str2fft'].apply(simu)
     distr['str2fft'].wait()
     op['diffusion'].apply(simu)             # Diffusion
-- 
GitLab