Inegrate Scales multi-scales advection with interpolation (still bug in parallel)

HySoP/hysop/f2py/scales2py.f90

 !> Interface to scales-advection solver
 module scales2py
 
-use cart_topology, only : cart_create,set_group_size,discretisation_create,N_proc,coord,cart_rank
-use advec, only : advec_init,advec_step
+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 interpolation_velo, only : interpol_init
 use mpi
 use parmesparam
 
@@ -17,7 +18,7 @@ contains
   !! @param[out] datashape local dimension of the input/output field
   !! @param[out] offset absolute index of the first component of the local field
   subroutine init_advection_solver(ncells,lengths,topodims,datashape,offset,dim,order,stab_coeff,dim_split)
-    integer, intent(in) :: dim 
+    integer, intent(in) :: dim
     integer, dimension(dim),intent(in) :: ncells
     real(pk),dimension(dim), intent(in) :: lengths
     integer, dimension(dim), intent(in) :: topodims
@@ -44,18 +45,37 @@ contains
     !call set_group_size(groupSize)
     ! Create meshes
     call discretisation_create(ncells(1),ncells(2),ncells(3),lengths(1),lengths(2),lengths(3))
-    
+
     ! Init advection solver
     call advec_init(order,stab_coeff,dim_split=dim_split)
-    
+
     ! get the local resolution (saved in scales global variable "N_proc")
     datashape = N_proc + 1
     ! get offset (i.e. global index of the lowest point of the current subdomain == scales global var "coord")
     offset = coord
-    
+
   end subroutine init_advection_solver
 
-  !> Particular solver for the advection of a scalar 
+  !> To change velocity resolution
+  !!    @param[in] Nx   = number of points along X
+  !!    @param[in] Ny   = number of points along Y
+  !!    @param[in] Nz   = number of points along Z
+  !!    @param[in] formula   = interpolation formula to use ('lin', 'L4_4' ,'M4')
+  subroutine init_multiscale(Nx, Ny, Nz, formula)
+
+    integer, intent(in) :: Nx, Ny, Nz
+    character(len=*), optional, intent(in)  ::  formula
+
+    call discretisation_set_mesh_Velo(Nx, Ny, Nz)
+    if(present(formula)) then
+      call interpol_init(formula, .true.)
+    else
+      call interpol_init('L4_4', .true.)
+    end if
+
+  end subroutine init_multiscale
+
+  !> 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
@@ -75,4 +95,41 @@ contains
     !!print *, "inside ...", cart_rank, ":", MPI_Wtime()-t0
   end subroutine solve_advection
 
+
+  !> 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_basic(dt, Vx, Vy, Vz, scal)
+
+    ! 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
+
+    call advec_step_Inter_basic(dt, Vx, Vy, Vz, scal)
+
+  end subroutine solve_advection_inter_basic
+
+  !> Solve advection equation - order 2 - with improved 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)
+
+    ! 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
+
+    call advec_step_Inter_Two(dt, Vx, Vy, Vz, scal)
+
+  end subroutine solve_advection_inter
+
 end module scales2py

HySoP/hysop/operator/advection.py

@@ -233,6 +233,11 @@ class Advection(Operator):
                                                         ghosts=self.ghosts)
                     # ... and the corresponding discrete field
                     self.discreteFields[v] = v.discretize(topo)
+            if self._isMultiScale:
+                self.config['isMultiscale'] = self._isMultiScale
+                v_shape = self.discreteFields[self.velocity].data[0].shape
+                scales.init_multiscale(v_shape[0], v_shape[1], v_shape[2],
+                                       self.method[MultiScale])
             self._my_setUp = self.setUp_Scales
 
         # --- GPU or pure-python advection ---

HySoP/hysop/operator/discrete/scales_advection.py

@@ -21,7 +21,8 @@ class ScalesAdvection(DiscreteOperator):
     """
 
     @debug
-    def __init__(self, velocity, advectedFields, method=''):
+    def __init__(self, velocity, advectedFields, method='',
+                 isMultiscale=False):
         """
         Constructor.
         @param velocity discrete field
@@ -43,6 +44,16 @@ 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
+
     @debug
     def apply(self, simulation=None):
         if simulation is None:
@@ -57,7 +68,7 @@ class ScalesAdvection(DiscreteOperator):
         # Call scales advection
         for adF in self.advectedFields:
             for d in xrange(adF.nbComponents):
-                adF[d][...] = scales2py.solve_advection(
+                adF[d][...] = self.scales_func(
                     dt, self.velocity.data[0],
                     self.velocity.data[1],
                     self.velocity.data[2],

HySoP/src/scalesInterface/particles/advec.f90

@@ -216,8 +216,6 @@ end subroutine advec_setup_alongZ
 !!    @param[in,out]    scal        = scalar field to advect
 subroutine advec_step_Inter_basic(dt, Vx, Vy, Vz, scal)
 
-    use cart_topology, only : mesh_sc, mesh_V, cart_rank
-
     ! Input/Output
     real(WP), intent(in)                        :: dt
     real(WP), dimension(:,:,:), intent(in)      :: Vx, Vy, Vz
@@ -255,8 +253,6 @@ end subroutine advec_step_Inter_basic
 !!    @param[in,out]    scal        = scalar field to advect
 subroutine advec_step_Inter_Two(dt, Vx, Vy, Vz, scal)
 
-    use cart_topology, only : mesh_sc, mesh_V, cart_rank
-
     ! Input/Output
     real(WP), intent(in)                        :: dt
     real(WP), dimension(:,:,:), intent(in)      :: Vx, Vy, Vz

HySoP/src/scalesInterface/particles/interpolation_velo.f90

@@ -619,7 +619,7 @@ subroutine Inter_LastDir_Permut_com(dir, V_coarse, dx_c, V_fine, dx_f)
   com_size(1) = size(V_coarse,1)*size(V_coarse,2)
   N_coarse = size(V_coarse,3)
   N_fine = size(V_fine,2)
-ind_max = ceiling(((N_coarse-stencil_d)*dx_c/dx_f)-1e-6)  ! To avoid numerical error and thus segmentation fault
+  ind_max = ceiling(((N_coarse-stencil_d)*dx_c/dx_f)-1e-6)  ! To avoid numerical error and thus segmentation fault
   ind_min = ceiling((stencil_g)*dx_c/dx_f)+1
 
   ! ==== Communication ====
@@ -1021,7 +1021,7 @@ end subroutine weight_Lambda4_4
 
 
 !> Basic interpolation formula. Be careful, this kernel may create unphysical oscillation
-!! in low frequency. This is rather implemented to show the requirement of "better" 
+!! in low frequency. This is rather implemented to show the requirement of "better"
 !! interpolation kernel.
 subroutine weight_linear(pos,weight)