Add interface to Scales Vector advection. Some bug fix.

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
  

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

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)

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):

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):
         """

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:

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:

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 ..."

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.

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