From ae062d802d233c99595df92140c3f29ffc11c798 Mon Sep 17 00:00:00 2001
From: Juan Escobar <escj@aero.obs-mip.fr>
Date: Fri, 22 Mar 2019 19:07:24 +0100
Subject: [PATCH] =?UTF-8?q?Juan=2022/03/2019:=20MG=20solve=20fonctionne=20?=
=?UTF-8?q?correctement=20grace=20a=20boundary=5Fmnh=20pour=20les=20condit?=
=?UTF-8?q?ion=20newman=20:=20test=C3=A9=20en=20RICH=20/=20Jacobi=20/=20tr?=
=?UTF-8?q?i-linear=20/=20cell=20average?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit
---
.../communication.f90 | 36 +++++++++++++++++++
.../discretisation.f90 | 5 ++-
tensorproductmultigrid_Source/multigrid.f90 | 2 +-
3 files changed, 41 insertions(+), 2 deletions(-)
diff --git a/tensorproductmultigrid_Source/communication.f90 b/tensorproductmultigrid_Source/communication.f90
index 3f7c44e2f..0739f83a0 100644
--- a/tensorproductmultigrid_Source/communication.f90
+++ b/tensorproductmultigrid_Source/communication.f90
@@ -48,6 +48,7 @@ public::comm_preinitialise
public::comm_initialise
public::comm_finalise
public::scalarprod
+public::boundary_mnh
public::haloswap_mnh
public::haloswap
public::ihaloswap_mnh
@@ -425,6 +426,39 @@ contains
s = global_sum
end subroutine scalarprod
!==================================================================
+! Boundary Neumann
+!==================================================================
+ subroutine boundary_mnh(a) ! data field
+
+ implicit none
+
+ type(scalar3d), intent(inout) :: a
+ !local var
+
+
+ integer :: n, ix_min,ix_max,iy_min,iy_max
+
+
+ ! Update Real Boundary for Newman case u(0) = u(1) , etc ...
+ n = a%grid_param%n
+ ix_min = a%ix_min
+ ix_max = a%ix_max
+ iy_min = a%iy_min
+ iy_max = a%iy_max
+ if ( ix_min == 1 ) then
+ a%s(:,:,0) = a%s(:,:,1)
+ endif
+ if ( ix_max == n ) then
+ a%s(:,:,a%icompx_max+1) = a%s(:,:,a%icompx_max)
+ endif
+ if ( iy_min == 1 ) then
+ a%s(:,0,:) = a%s(:,1,:)
+ endif
+ if ( iy_max == n ) then
+ a%s(:,a%icompy_max+1,:) = a%s(:,a%icompy_max,:)
+ endif
+ end subroutine boundary_mnh
+!==================================================================
! Initiate asynchronous halo exchange
!
! For all processes with horizontal indices that are multiples
@@ -666,6 +700,7 @@ contains
integer, intent(in) :: level
integer, intent(in) :: m
type(scalar3d), intent(inout) :: a
+ !local var
integer :: a_n ! horizontal grid size
integer :: nz ! vertical grid size
integer, dimension(2) :: p_horiz
@@ -769,6 +804,7 @@ contains
call finish_timer(t_haloswap(level,m))
end if
end if
+
end subroutine haloswap_mnh
!==================================================================
subroutine haloswap(level,m, & ! multigrid- and processor- level
diff --git a/tensorproductmultigrid_Source/discretisation.f90 b/tensorproductmultigrid_Source/discretisation.f90
index 4962aeb9e..92a850e13 100644
--- a/tensorproductmultigrid_Source/discretisation.f90
+++ b/tensorproductmultigrid_Source/discretisation.f90
@@ -302,7 +302,7 @@ contains
h = grid_param%L/grid_param%n
! Cartesian coefficients
Tij = h**2
- xcoef = 0.0
+ xcoef = 0.5_rl ! 0.0
l2nd = .false. ! .true. ! .false.
alpha_T(1) = 1.0
alpha_T(2) = 1.0
@@ -531,6 +531,7 @@ end subroutine construct_vertical_coeff
integer :: ix,iy,iz
! r <- A.u
+ call boundary_mnh(u)
call apply_mnh(u,r)
! r <- b - r = b - A.u
do ix=u%icompx_min,u%icompx_max
@@ -1147,6 +1148,8 @@ end subroutine construct_vertical_coeff
integer :: iymin(5), iymax(5)
integer :: iblock, ierr
+ call boundary_mnh(u)
+
! Set optimal smoothing parameter on each level
rho = 2.0_rl/(2.0_rl+4.0_rl*model_param%omega2*u%grid_param%n**2/(1.0_rl+4.0_rl*model_param%omega2*u%grid_param%n**2))
diff --git a/tensorproductmultigrid_Source/multigrid.f90 b/tensorproductmultigrid_Source/multigrid.f90
index 06f61526a..d026a69a7 100644
--- a/tensorproductmultigrid_Source/multigrid.f90
+++ b/tensorproductmultigrid_Source/multigrid.f90
@@ -660,7 +660,7 @@ contains
!local var
integer :: ix,iy,iz
real (kind=rl) :: xcoef
- xcoef = 0.0
+ xcoef = 0.5_rl ! 0.0
do ix=ixmin(iblock),ixmax(iblock)
do iy=iymin(iblock),iymax(iblock)
#ifdef PIECEWISELINEAR
--
GitLab