Juan add restriction=REST_KHALIL ( en fait Kwak <-> a Corriger )

tensorproductmultigrid_Source/mode_mg_read_param.f90

@@ -419,6 +419,8 @@ subroutine read_multigrid_parameters(filename,mg_param)
     write(STDOUT,'("    n_postsmooth = ",I6)') n_postsmooth
     if (restriction == REST_CELLAVERAGE) then
       write(STDOUT,'("    restriction   = CELLAVERAGE")')
+   else if (restriction == REST_KHALIL) then
+      write(STDOUT,'("    restriction   = KHALIL ")')
     else
       restriction = -1
     endif

tensorproductmultigrid_Source/multigrid.f90

@@ -56,6 +56,7 @@ public::mg_solve_mnh
 public::mg_solve
 public::measurehaloswap
 public::REST_CELLAVERAGE
+public::REST_KHALIL
 public::PROL_CONSTANT
 public::PROL_TRILINEAR
 public::COARSEGRIDSOLVER_SMOOTHER
@@ -66,6 +67,7 @@ private
   ! --- multigrid parameter constants ---
   ! restriction
   integer, parameter :: REST_CELLAVERAGE = 1
+  integer, parameter :: REST_KHALIL      = 2
   ! prolongation method
   integer, parameter :: PROL_CONSTANT = 1
   integer, parameter :: PROL_TRILINEAR = 2
@@ -418,6 +420,66 @@ contains
 
 !==================================================================
 ! Restrict from fine -> coarse
+!==================================================================
+  subroutine restrict_mnh(phifine,phicoarse)
+    implicit none
+    type(scalar3d), intent(in)    :: phifine
+    type(scalar3d), intent(inout) :: phicoarse
+    integer :: ix,iy,iz
+    integer :: ix_min, ix_max, iy_min, iy_max
+    real(kind=rl) :: xn,xs,xw,xe
+
+    ix_min = phicoarse%icompx_min
+    ix_max = phicoarse%icompx_max
+    iy_min = phicoarse%icompy_min
+    iy_max = phicoarse%icompy_max
+    ! three dimensional cell average
+    if (mg_param%restriction == REST_CELLAVERAGE) then
+      ! Do not coarsen in z-direction
+      do ix=ix_min,ix_max
+        do iy=iy_min,iy_max
+          do iz=1,phicoarse%grid_param%nz
+            phicoarse%s(iz,iy,ix) =  &
+              phifine%s(iz  ,2*iy  ,2*ix  ) + &
+              phifine%s(iz  ,2*iy-1,2*ix  ) + &
+              phifine%s(iz  ,2*iy  ,2*ix-1) + &
+              phifine%s(iz  ,2*iy-1,2*ix-1)
+          end do
+        end do
+      end do
+      elseif(mg_param%restriction == REST_KHALIL) then
+     ! Do not coarsen in z-direction
+      do ix=ix_min,ix_max
+         xw=1.0
+         xe=1.0
+         if (ix==phicoarse%ix_min) xw=0.0
+         if (ix==phicoarse%ix_max) xe=0.0
+        do iy=iy_min,iy_max
+           xs=1.0
+           xn=1.0
+           if (iy==phicoarse%iy_min) xs=0.0
+           if (iy==phicoarse%iy_max) xn=0.0
+          do iz=1,phicoarse%grid_param%nz
+            phicoarse%s(iz,iy,ix) = 0.25_rl *         ( &
+              phifine%s(iz,2*iy+1,2*ix-1) * xn        + &
+              phifine%s(iz,2*iy+1,2*ix  ) * xn        + &
+              phifine%s(iz,2*iy  ,2*ix-2) * xw        + &
+              phifine%s(iz,2*iy  ,2*ix-1) * (4-xw-xn) + &
+              phifine%s(iz,2*iy  ,2*ix  ) * (4-xe-xn) + &
+              phifine%s(iz,2*iy  ,2*ix+1) * xe        + &
+              phifine%s(iz,2*iy-1,2*ix-2) * xw        + &
+              phifine%s(iz,2*iy-1,2*ix-1) * (4-xw-xs) + &
+              phifine%s(iz,2*iy-1,2*ix  ) * (4-xe-xs) + &
+              phifine%s(iz,2*iy-2,2*ix-1) * xs        + &
+              phifine%s(iz,2*iy-2,2*ix  ) * xs          &     
+                                                        & )          
+          end do
+        end do
+      end do
+    end if
+  end subroutine restrict_mnh
+!==================================================================
+! Restrict from fine -> coarse
 !==================================================================
   subroutine restrict(phifine,phicoarse)
     implicit none
@@ -1125,7 +1187,7 @@ contains
       ! Restrict residual
       call start_timer(t_restrict(level,m))
       call start_timer(t_total(level,m))
-      call restrict(r(level,m),b(level-1,m))
+      call restrict_mnh(r(level,m),b(level-1,m))
       call finish_timer(t_total(level,m))
       call finish_timer(t_restrict(level,m))
       if ( ((level-1) .le. splitlevel) .and. (m > 0) ) then