From 05b245bb6354b63c0a0642943f4f5994779c7a30 Mon Sep 17 00:00:00 2001
From: Jean-Matthieu Etancelin <jean-matthieu.etancelin@imag.fr>
Date: Fri, 28 Feb 2014 15:22:39 +0000
Subject: [PATCH] Add interface to Scales Vector advection. Some bug fix.
---
HySoP/ParmesToSinglePrecision.patch | 8 +-
HySoP/hysop/f2py/scales2py.f90 | 44 +++++++++--
.../gpu/tests/test_opencl_environment.py | 73 +++++++++++--------
.../operator/discrete/scales_advection.py | 41 +++++++----
HySoP/hysop/operator/monitors/monitoring.py | 3 +-
HySoP/hysop/operator/monitors/printer.py | 3 +-
HySoP/hysop/problem/problem.py | 2 +-
HySoP/hysop/problem/problem_tasks.py | 2 +-
HySoP/hysop/problem/simulation.py | 3 -
.../particles/advec_Vector.f90 | 39 ++++++++++
10 files changed, 153 insertions(+), 65 deletions(-)
diff --git a/HySoP/ParmesToSinglePrecision.patch b/HySoP/ParmesToSinglePrecision.patch
index 61079e7b3..d3a9075d4 100644
--- a/HySoP/ParmesToSinglePrecision.patch
+++ b/HySoP/ParmesToSinglePrecision.patch
@@ -51,11 +51,11 @@ index 8645a9f..718eb19 100755
start = MPI_WTime()
call r2c_3d(omega_x,omega_y,omega_z, ghosts_vort)
diff --git parmepy/f2py/scales2py.f90 parmepy/f2py/scales2py.f90
-index e574369..391e127 100755
+index cca8fea..0abc5cd 100755
--- parmepy/f2py/scales2py.f90
+++ parmepy/f2py/scales2py.f90
-@@ -5,7 +5,7 @@ use cart_topology, only : cart_create,set_group_size,discretisation_create,N_pro
- use advec, only : advec_init,advec_step,advec_step_Inter_basic,advec_step_Inter_Two
+@@ -6,7 +6,7 @@ use advec, only : advec_init,advec_step,advec_step_Inter_basic,advec_step_Inter_
+ use advec_vect, only : advec_step_Vect,advec_step_Inter_basic_Vect
use interpolation_velo, only : interpol_init
use mpi
-use parmesparam
@@ -63,7 +63,7 @@ index e574369..391e127 100755
implicit none
-@@ -91,7 +91,7 @@ contains
+@@ -92,7 +92,7 @@ contains
real(pk), dimension(size(vx,1),size(vx,2),size(vx,3)), intent(inout) :: scal
!f2py real(pk) intent(in,out), depend(size(vx,1)) :: scal
diff --git a/HySoP/hysop/f2py/scales2py.f90 b/HySoP/hysop/f2py/scales2py.f90
index e57436974..cca8fea00 100755
--- a/HySoP/hysop/f2py/scales2py.f90
+++ b/HySoP/hysop/f2py/scales2py.f90
@@ -3,6 +3,7 @@ module scales2py
use cart_topology, only : cart_create,set_group_size,discretisation_create,N_proc,coord,cart_rank,discretisation_set_mesh_Velo
use advec, only : advec_init,advec_step,advec_step_Inter_basic,advec_step_Inter_Two
+use advec_vect, only : advec_step_Vect,advec_step_Inter_basic_Vect
use interpolation_velo, only : interpol_init
use mpi
use parmesparam
@@ -98,6 +99,30 @@ contains
!!print *, "inside ...", cart_rank, ":", MPI_Wtime()-t0
end subroutine solve_advection
+ !> Particular solver for the advection of a scalar
+ !! @param[in] dt current time step
+ !! @param[in] vx x component of the velocity used for advection
+ !! @param[in] vy y component of the velocity used for advection
+ !! @param[in] vz z component of the velocity used for advection
+ !! @param[in,out] cx 3d scalar field which is advected
+ !! @param[in,out] cy 3d scalar field which is advected
+ !! @param[in,out] cz 3d scalar field which is advected
+ subroutine solve_advection_vect(dt,vx,vy,vz,cx,cy,cz)
+
+ real(pk), intent(in) :: dt
+ real(pk), dimension(:,:,:), intent(in) :: vx, vy, vz
+ real(pk), dimension(size(vx,1),size(vx,2),size(vx,3)), intent(inout) :: cx,cy,cz
+ !f2py real(pk) intent(in,out), depend(size(vx,1)) :: cx
+ !f2py real(pk) intent(in,out), depend(size(vx,1)) :: cy
+ !f2py real(pk) intent(in,out), depend(size(vx,1)) :: cz
+
+ real(8) :: t0
+
+ t0 = MPI_Wtime()
+ call advec_step_Vect(dt,vx,vy,vz,cx,cy,cz)
+ !!print *, "inside ...", cart_rank, ":", MPI_Wtime()-t0
+ end subroutine solve_advection_vect
+
!> Solve advection equation - order 2 - with basic velocity interpolation
!! @param[in] dt = time step
@@ -117,22 +142,27 @@ contains
end subroutine solve_advection_inter_basic
- !> Solve advection equation - order 2 - with improved velocity interpolation
+ !> Solve advection equation - order 2 - with basic velocity interpolation
!! @param[in] dt = time step
!! @param[in] Vx = velocity along x (could be discretised on a bigger mesh then the scalar)
!! @param[in] Vy = velocity along y
!! @param[in] Vz = velocity along z
- !! @param[in,out] scal = scalar field to advect
- subroutine solve_advection_inter(dt, Vx, Vy, Vz, scal)
+ !! @param[in,out] cx 3d scalar field which is advected
+ !! @param[in,out] cy 3d scalar field which is advected
+ !! @param[in,out] cz 3d scalar field which is advected
+ subroutine solve_advection_inter_basic_vec(dt, Vx, Vy, Vz, cx, cy, cz)
! Input/Output
real(pk), intent(in) :: dt
real(pk), dimension(:,:,:), intent(in) :: Vx, Vy, Vz
- real(pk), dimension(:,:,:), intent(inout) :: scal
- !f2py intent(in,out) :: scal
+ real(pk), dimension(:,:,:), intent(inout) :: cx,cy,cz
+ !f2py intent(in,out) :: cx
+ !f2py intent(in,out) :: cy
+ !f2py intent(in,out) :: cz
- call advec_step_Inter_Two(dt, Vx, Vy, Vz, scal)
+ call advec_step_Inter_basic_Vect(dt, Vx, Vy, Vz, cx, cy, cz)
+
+ end subroutine solve_advection_inter_basic
- end subroutine solve_advection_inter
end module scales2py
diff --git a/HySoP/hysop/gpu/tests/test_opencl_environment.py b/HySoP/hysop/gpu/tests/test_opencl_environment.py
index 46931e9af..25e061d1d 100644
--- a/HySoP/hysop/gpu/tests/test_opencl_environment.py
+++ b/HySoP/hysop/gpu/tests/test_opencl_environment.py
@@ -21,13 +21,15 @@ def test_parse_src_expand_floatN():
import StringIO
cl_env = get_opencl_environment(0, 0, 'gpu', FLOAT_GPU,)
str_as_src = """
- vstore__N__((float__N__)(gscal_loc[noBC_id(i+__NN__,nb_part)],
- ), (i + gidY*WIDTH)/__N__, gscal);
- """
+ vstore__N__((float__N__)(gscal_loc[noBC_id(i+__NN__,nb_part)],
+ ), (i + gidY*WIDTH)/__N__, gscal);
+ """
parsed_str_as_src = """
- vstore4((float4)(gscal_loc[noBC_id(i+0,nb_part)],gscal_loc[noBC_id(i+1,nb_part)],gscal_loc[noBC_id(i+2,nb_part)],gscal_loc[noBC_id(i+3,nb_part)]
- ), (i + gidY*WIDTH)/4, gscal);
- """
+ vstore4((float4)(gscal_loc[noBC_id(i+0,nb_part)],""" + \
+ """gscal_loc[noBC_id(i+1,nb_part)],""" + \
+ """gscal_loc[noBC_id(i+2,nb_part)],gscal_loc[noBC_id(i+3,nb_part)]
+ ), (i + gidY*WIDTH)/4, gscal);
+ """
buf = StringIO.StringIO(str_as_src)
res = cl_env.parse_file(buf, n=4)
assert len(parsed_str_as_src) == len(res)
@@ -62,52 +64,61 @@ def test_parse_expand_remeshed_component():
import StringIO
cl_env = get_opencl_environment(0, 0, 'gpu', FLOAT_GPU)
str_as_src = """
-__kernel void advection_and_remeshing(__global const float* gvelo,
- __RCOMP_P__global const float* pscal__ID__,
- __RCOMP_P__global float* gscal__ID__,
- __local float* gvelo_loc,
- __RCOMP_P__local float* gscal_loc__ID__,
- float dt,float min_position, float dx)
+ __kernel void advection_and_remeshing(__global const float* gvelo,
+ __RCOMP_P__global const float* pscal__ID__,
+ __RCOMP_P__global float* gscal__ID__,
+ __local float* gvelo_loc,
+ __RCOMP_P__local float* gscal_loc__ID__,
+ float dt,float min_position, float dx)
{
__RCOMP_I gscal_loc__ID__[noBC_id(i)] = 0.0;
remesh(i, dx, invdx, s, p, __RCOMP_Pgscal_loc__ID__);
test(__RCOMP_Pgscal_loc__ID__, __RCOMP_Ppscal__ID__);
__RCOMP_I gscal__ID__[i + line_index] = gscal_loc__ID__[noBC_id(i)];
__RCOMP_I vstore__N__((float__N__)(gscal_loc__ID__[noBC_id(i+__NN__)],
- ), (i + line_index)/__N__, gscal__ID__);
+ ), (i + line_index)/__N__, gscal__ID__);
"""
parsed_str_as_src_2components = """
-__kernel void advection_and_remeshing(__global const float* gvelo,
- __global const float* pscal0, __global const float* pscal1,
- __global float* gscal0, __global float* gscal1,
- __local float* gvelo_loc,
- __local float* gscal_loc0, __local float* gscal_loc1,
- float dt,float min_position, float dx)
+ __kernel void advection_and_remeshing(__global const float* gvelo,
+ """ + \
+ """__global const float* pscal0, __global const float* pscal1,
+ __global float* gscal0, __global float* gscal1,
+ __local float* gvelo_loc,
+ __local float* gscal_loc0, __local float* gscal_loc1,
+ float dt,float min_position, float dx)
{
gscal_loc0[noBC_id(i)] = 0.0; gscal_loc1[noBC_id(i)] = 0.0;
remesh(i, dx, invdx, s, p, gscal_loc0, gscal_loc1);
test(gscal_loc0, gscal_loc1, pscal0, pscal1);
- gscal0[i + line_index] = gscal_loc0[noBC_id(i)]; gscal1[i + line_index] = gscal_loc1[noBC_id(i)];
- vstore4((float4)(gscal_loc0[noBC_id(i+0)],gscal_loc0[noBC_id(i+1)],gscal_loc0[noBC_id(i+2)],gscal_loc0[noBC_id(i+3)]
- ), (i + line_index)/4, gscal0); vstore4((float4)(gscal_loc1[noBC_id(i+0)],gscal_loc1[noBC_id(i+1)],gscal_loc1[noBC_id(i+2)],gscal_loc1[noBC_id(i+3)]
- ), (i + line_index)/4, gscal1);
+ gscal0[i + line_index] = gscal_loc0[noBC_id(i)]; """ + \
+ """gscal1[i + line_index] = gscal_loc1[noBC_id(i)];
+ vstore4((float4)(gscal_loc0[noBC_id(i+0)],""" + \
+ """gscal_loc0[noBC_id(i+1)],gscal_loc0[noBC_id(i+2)],""" + \
+ """gscal_loc0[noBC_id(i+3)]
+ ), (i + line_index)/4, gscal0); """ + \
+ """vstore4((float4)(gscal_loc1[noBC_id(i+0)],""" + \
+ """gscal_loc1[noBC_id(i+1)],gscal_loc1[noBC_id(i+2)],""" + \
+ """gscal_loc1[noBC_id(i+3)]
+ ), (i + line_index)/4, gscal1);
"""
parsed_str_as_src_1components = """
-__kernel void advection_and_remeshing(__global const float* gvelo,
- __global const float* pscal0,
- __global float* gscal0,
- __local float* gvelo_loc,
- __local float* gscal_loc0,
- float dt,float min_position, float dx)
+ __kernel void advection_and_remeshing(__global const float* gvelo,
+ __global const float* pscal0,
+ __global float* gscal0,
+ __local float* gvelo_loc,
+ __local float* gscal_loc0,
+ float dt,float min_position, float dx)
{
gscal_loc0[noBC_id(i)] = 0.0;
remesh(i, dx, invdx, s, p, gscal_loc0);
test(gscal_loc0, pscal0);
gscal0[i + line_index] = gscal_loc0[noBC_id(i)];
- vstore4((float4)(gscal_loc0[noBC_id(i+0)],gscal_loc0[noBC_id(i+1)],gscal_loc0[noBC_id(i+2)],gscal_loc0[noBC_id(i+3)]
- ), (i + line_index)/4, gscal0);
+ vstore4((float4)(gscal_loc0[noBC_id(i+0)],""" + \
+ """gscal_loc0[noBC_id(i+1)],gscal_loc0[noBC_id(i+2)],""" + \
+ """gscal_loc0[noBC_id(i+3)]
+ ), (i + line_index)/4, gscal0);
"""
buf = StringIO.StringIO(str_as_src)
diff --git a/HySoP/hysop/operator/discrete/scales_advection.py b/HySoP/hysop/operator/discrete/scales_advection.py
index 98ff425f3..08bca68a4 100644
--- a/HySoP/hysop/operator/discrete/scales_advection.py
+++ b/HySoP/hysop/operator/discrete/scales_advection.py
@@ -44,15 +44,26 @@ class ScalesAdvection(DiscreteOperator):
self.input = [self.velocity]
self.output = list(variables).remove(self.velocity)
- if isMultiscale:
- # 3D interpolation of the velocity before advection
- self.scales_func = scales2py.solve_advection_inter_basic
- # Other interpolation only 2D interpolation first and
- # 1D interpolations before advections in each direction
- # (slower than basic)
- # self.scales_func = scales2py.solve_advection_inter
- else:
- self.scales_func = scales2py.solve_advection
+ # Scales functions for each field (depending if vector)
+ self._scales_func = []
+
+ for adF in self.advectedFields:
+ if adF.nbComponents == 3:
+ if isMultiscale:
+ # 3D interpolation of the velocity before advection
+ self._scales_func.append(
+ scales2py.solve_advection_inter_basic_vect)
+ # Other interpolation only 2D interpolation first and
+ # 1D interpolations before advections in each direction
+ # (slower than basic): solve_advection_inter
+ else:
+ self._scales_func.append(scales2py.solve_advection_vect)
+ else:
+ if isMultiscale:
+ self._scales_func.append(
+ scales2py.solve_advection_inter_basic)
+ else:
+ self._scales_func.append(scales2py.solve_advection)
@debug
def apply(self, simulation=None):
@@ -66,13 +77,11 @@ class ScalesAdvection(DiscreteOperator):
dt = simulation.timeStep
# Call scales advection
- for adF in self.advectedFields:
- for d in xrange(adF.nbComponents):
- adF[d][...] = self.scales_func(
- dt, self.velocity.data[0],
- self.velocity.data[1],
- self.velocity.data[2],
- adF[d][...])
+ for adF, fun in zip(self.advectedFields, self._scales_func):
+ adF = fun(dt, self.velocity.data[0],
+ self.velocity.data[1],
+ self.velocity.data[2],
+ *adF)
self._apply_timer.append_time(MPI.Wtime() - ctime)
def finalize(self):
diff --git a/HySoP/hysop/operator/monitors/monitoring.py b/HySoP/hysop/operator/monitors/monitoring.py
index c96932fdd..19d728506 100644
--- a/HySoP/hysop/operator/monitors/monitoring.py
+++ b/HySoP/hysop/operator/monitors/monitoring.py
@@ -34,7 +34,8 @@ class Monitoring(Operator):
if io_params is None:
self._writer = None
else:
- self._writer = io.Writer(io_params, topo.comm)
+ if topo is not None:
+ self._writer = io.Writer(io_params, topo.comm)
def setUp(self):
"""
diff --git a/HySoP/hysop/operator/monitors/printer.py b/HySoP/hysop/operator/monitors/printer.py
index 35022c94c..28dc1b7c8 100644
--- a/HySoP/hysop/operator/monitors/printer.py
+++ b/HySoP/hysop/operator/monitors/printer.py
@@ -71,7 +71,8 @@ class Printer(Monitoring):
prefix = os.path.join(io.io.defaultPath(), prefix)
## shortcut to topo name
self.topo = self._predefinedTopo
- io.io.checkDir(prefix, 0, self.topo.comm)
+ if self.topo is not None:
+ io.io.checkDir(prefix, 0, self.topo.comm)
self.prefix = prefix
self._xmf_data_files = []
if self.formattype == VTK:
diff --git a/HySoP/hysop/problem/problem.py b/HySoP/hysop/problem/problem.py
index 5483c7dd1..6b8d64bc5 100644
--- a/HySoP/hysop/problem/problem.py
+++ b/HySoP/hysop/problem/problem.py
@@ -177,7 +177,7 @@ class Problem(object):
At end of time step, call an io step.\n
Displays timings at simulation end.
"""
- self.simulation.init()
+ self.simulation.initialize()
if main_rank == 0:
print "\n\n Start solving ..."
while not self.simulation.isOver:
diff --git a/HySoP/hysop/problem/problem_tasks.py b/HySoP/hysop/problem/problem_tasks.py
index 3ec307560..d2da7077a 100644
--- a/HySoP/hysop/problem/problem_tasks.py
+++ b/HySoP/hysop/problem/problem_tasks.py
@@ -173,7 +173,7 @@ class ProblemTasks(Problem):
At end of time step, call an io step.\n
Displays timings at simulation end.
"""
- self.simulation.init()
+ self.simulation.initialize()
self.main_comm.Barrier()
if self._main_rank == 0:
print "\n\n Start solving ..."
diff --git a/HySoP/hysop/problem/simulation.py b/HySoP/hysop/problem/simulation.py
index fb16b8076..50265dac2 100644
--- a/HySoP/hysop/problem/simulation.py
+++ b/HySoP/hysop/problem/simulation.py
@@ -58,9 +58,6 @@ class Simulation(object):
assert (self.start + self.timeStep) <= self.end,\
'start + step is bigger than end.'
- def init(self):
- self.currentIteration = 1
-
def advance(self):
"""
Proceed to next time.
diff --git a/HySoP/src/scalesInterface/particles/advec_Vector.f90 b/HySoP/src/scalesInterface/particles/advec_Vector.f90
index a6722a535..573f3488f 100644
--- a/HySoP/src/scalesInterface/particles/advec_Vector.f90
+++ b/HySoP/src/scalesInterface/particles/advec_Vector.f90
@@ -62,6 +62,45 @@ function type_part_solver()
type_part_solver = type_part_solv
end function
+!> Solve advection equation - order 2 - with basic velocity interpolation
+!! @param[in] dt = time step
+!! @param[in] Vx = velocity along x (could be discretised on a bigger mesh then the scalar)
+!! @param[in] Vy = velocity along y
+!! @param[in] Vz = velocity along z
+!! @param[in,out] VectX = X component of vector to advect
+!! @param[in,out] VectY = Y component of vector to advect
+!! @param[in,out] VectZ = Z component of vector to advect
+subroutine advec_step_Inter_basic_Vect(dt, Vx, Vy, Vz, vectX, vectY, vectZ)
+
+ use Interpolation_velo, only : Interpol_3D
+
+ ! Input/Output
+ real(WP), intent(in) :: dt
+ real(WP), dimension(:,:,:), intent(in) :: Vx, Vy, Vz
+ real(WP), dimension(:,:,:), intent(inout) :: vectX, vectY, vectZ
+ ! Local
+ real(WP), dimension(:,:,:), allocatable :: Vx_f, Vy_f, Vz_f
+ integer :: ierr ! Error code.
+
+ allocate(Vx_f(mesh_sc%N_proc(1),mesh_sc%N_proc(2),mesh_sc%N_proc(3)),stat=ierr)
+ if (ierr/=0) write(6,'(a,i0,a)') '[ERROR]Â on cart_rank ', cart_rank, ' - not enough memory for Vx_f'
+ allocate(Vy_f(mesh_sc%N_proc(1),mesh_sc%N_proc(2),mesh_sc%N_proc(3)),stat=ierr)
+ if (ierr/=0) write(6,'(a,i0,a)') '[ERROR]Â on cart_rank ', cart_rank, ' - not enough memory for Vy_f'
+ allocate(Vz_f(mesh_sc%N_proc(1),mesh_sc%N_proc(2),mesh_sc%N_proc(3)),stat=ierr)
+ if (ierr/=0) write(6,'(a,i0,a)') '[ERROR]Â on cart_rank ', cart_rank, ' - not enough memory for Vz_f'
+
+ call Interpol_3D(Vx, mesh_V%dx, Vx_f, mesh_sc%dx)
+ call Interpol_3D(Vy, mesh_V%dx, Vy_f, mesh_sc%dx)
+ call Interpol_3D(Vz, mesh_V%dx, Vz_f, mesh_sc%dx)
+ if (cart_rank==0) write(6,'(a)') ' [INFO PARTICLES] Interpolation done'
+
+ call advec_step_Vect(dt, Vx_f, Vy_f, Vz_f, vectX, vectY, vectZ)
+
+ deallocate(Vx_f)
+ deallocate(Vy_f)
+ deallocate(Vz_f)
+
+end subroutine advec_step_Inter_basic_Vect
!> Solve advection equation - order 2 in time (order 2 dimensional splitting)
!! @param[in] dt = time step
--
GitLab