From 8e1bc1719ec23b26587e6f482e411e68b0512a81 Mon Sep 17 00:00:00 2001
From: Juan Escobar <escj@aero.obs-mip.fr>
Date: Tue, 19 Feb 2019 15:07:53 +0100
Subject: [PATCH] =?UTF-8?q?Juan=2019/02/2019=20:=20d=C3=A9but=20d=20integr?=
=?UTF-8?q?ation=20de=20MG=20dans=20MesoNH=20=3D>=20mg=5Fmain=5Fmnh.f90?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit
---
.../communication.f90 | 3 +-
.../conjugategradient.f90 | 3 +-
tensorproductmultigrid_Source/datatypes.f90 | 3 +-
.../discretisation.f90 | 3 +-
tensorproductmultigrid_Source/messages.f90 | 3 +-
tensorproductmultigrid_Source/mg_main.f90 | 3 +-
tensorproductmultigrid_Source/mg_main_mnh.f90 | 150 ++++++++++++++++++
tensorproductmultigrid_Source/mode_mg.f90 | 18 ++-
.../mode_mg_read_param.f90 | 28 ++--
tensorproductmultigrid_Source/multigrid.f90 | 4 +-
tensorproductmultigrid_Source/profiles.f90 | 43 ++++-
tensorproductmultigrid_Source/timer.f90 | 3 +-
tridz.f90 | 3 +
13 files changed, 241 insertions(+), 26 deletions(-)
create mode 100644 tensorproductmultigrid_Source/mg_main_mnh.f90
diff --git a/tensorproductmultigrid_Source/communication.f90 b/tensorproductmultigrid_Source/communication.f90
index f0856420b..c0b82d05c 100644
--- a/tensorproductmultigrid_Source/communication.f90
+++ b/tensorproductmultigrid_Source/communication.f90
@@ -38,7 +38,8 @@ module communication
use messages
use datatypes
use parameters
- use mpi
+ !use mpi
+ use modd_mpif
use timer
implicit none
diff --git a/tensorproductmultigrid_Source/conjugategradient.f90 b/tensorproductmultigrid_Source/conjugategradient.f90
index b58adae6b..eaa98b3a0 100644
--- a/tensorproductmultigrid_Source/conjugategradient.f90
+++ b/tensorproductmultigrid_Source/conjugategradient.f90
@@ -40,7 +40,8 @@ module conjugategradient
use discretisation
use messages
use communication
- use mpi
+ !use mpi
+ use modd_mpif
implicit none
diff --git a/tensorproductmultigrid_Source/datatypes.f90 b/tensorproductmultigrid_Source/datatypes.f90
index b0bfbecc3..4705e1201 100644
--- a/tensorproductmultigrid_Source/datatypes.f90
+++ b/tensorproductmultigrid_Source/datatypes.f90
@@ -39,7 +39,8 @@
module datatypes
- use mpi
+ !use mpi
+ use modd_mpif
use parameters
use messages
diff --git a/tensorproductmultigrid_Source/discretisation.f90 b/tensorproductmultigrid_Source/discretisation.f90
index 2a45884d2..35f8b0c89 100644
--- a/tensorproductmultigrid_Source/discretisation.f90
+++ b/tensorproductmultigrid_Source/discretisation.f90
@@ -65,7 +65,8 @@ module discretisation
use messages
use datatypes
use communication
- use mpi
+ !use mpi
+ use modd_mpif
implicit none
diff --git a/tensorproductmultigrid_Source/messages.f90 b/tensorproductmultigrid_Source/messages.f90
index 4ade93172..ce0fd76c3 100644
--- a/tensorproductmultigrid_Source/messages.f90
+++ b/tensorproductmultigrid_Source/messages.f90
@@ -36,7 +36,8 @@
module messages
use parameters
- use mpi
+ !use mpi
+ use modd_mpif
implicit none
diff --git a/tensorproductmultigrid_Source/mg_main.f90 b/tensorproductmultigrid_Source/mg_main.f90
index 4ff091121..0fef58742 100644
--- a/tensorproductmultigrid_Source/mg_main.f90
+++ b/tensorproductmultigrid_Source/mg_main.f90
@@ -51,7 +51,8 @@ program mg_main
use messages
use communication
use timer
- use mpi
+ !use mpi
+ use modd_mpif
use mode_mg_read_param
use mode_mg
diff --git a/tensorproductmultigrid_Source/mg_main_mnh.f90 b/tensorproductmultigrid_Source/mg_main_mnh.f90
new file mode 100644
index 000000000..e99d66fd5
--- /dev/null
+++ b/tensorproductmultigrid_Source/mg_main_mnh.f90
@@ -0,0 +1,150 @@
+!=== COPYRIGHT AND LICENSE STATEMENT ===
+!
+! This file is part of the TensorProductMultigrid code.
+!
+! (c) The copyright relating to this work is owned jointly by the
+! Crown, Met Office and NERC [2014]. However, it has been created
+! with the help of the GungHo Consortium, whose members are identified
+! at https://puma.nerc.ac.uk/trac/GungHo/wiki .
+!
+! Main Developer: Eike Mueller
+!
+! TensorProductMultigrid is free software: you can redistribute it and/or
+! modify it under the terms of the GNU Lesser General Public License as
+! published by the Free Software Foundation, either version 3 of the
+! License, or (at your option) any later version.
+!
+! TensorProductMultigrid is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU Lesser General Public License for more details.
+!
+! You should have received a copy of the GNU Lesser General Public License
+! along with TensorProductMultigrid (see files COPYING and COPYING.LESSER).
+! If not, see <http://www.gnu.org/licenses/>.
+!
+!=== COPYRIGHT AND LICENSE STATEMENT ===
+
+
+!==================================================================
+!
+! Main program for multigrid solver code for Helmholtz/Poisson
+! equation, discretised in the cell centred finite volume scheme
+!
+! Eike Mueller, University of Bath, Feb 2012
+!
+!==================================================================
+
+!==================================================================
+! Main program
+!==================================================================
+
+module mode_mg_main_mnh
+
+ use discretisation
+ use parameters
+ use datatypes
+ use multigrid
+ use conjugategradient
+ use solver
+ use profiles
+ use messages
+ use communication
+ use timer
+ !use mpi
+ use modd_mpif
+
+ use mode_mg_read_param
+ use mode_mg
+
+contains
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ subroutine mg_main_mnh_init()
+ implicit none
+
+ call mg_init_mnh()
+
+ ! Initialise ghosts in initial solution, as mg_solve assumes that they
+ ! are up-to-date
+ call haloswap(mg_param%n_lev,pproc,u)
+
+ end subroutine mg_main_mnh_init
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+subroutine mg_main_initialise_rhs_mnh(xcoef)
+
+implicit none
+
+real(kind=rl) :: xcoef
+
+ call initialise_rhs_mnh(grid_param,model_param,b,xcoef)
+
+end subroutine mg_main_initialise_rhs_mnh
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+subroutine mg_main_initialise_u_mnh(xcoef)
+
+implicit none
+
+real(kind=rl) :: xcoef
+
+ u%s(:,:,:) = xcoef
+
+end subroutine mg_main_initialise_u_mnh
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ subroutine mg_main_mnh_solve()
+
+ implicit none
+
+ ! Solve using multigrid
+ call initialise_timer(t_solve,"t_solve")
+ call start_timer(t_solve)
+ comm_measuretime = .True.
+
+ call mg_solve_mnh(b,u,solver_param)
+
+ comm_measuretime = .False.
+ call finish_timer(t_solve)
+
+ end subroutine mg_main_mnh_solve
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ subroutine mg_main_mnh_finalize()
+
+
+ end subroutine mg_main_mnh_finalize
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ subroutine mg_main_mnh
+
+ implicit none
+
+ integer :: it
+
+ call mg_main_mnh_init()
+
+ DO it=1,10
+
+ call mg_main_initialise_rhs_mnh(1.0*2.0**it)
+ call mg_main_initialise_u_mnh(1.0*it)
+
+ print*,'************************ IT=',it,' ***************************'
+
+ call mg_main_mnh_solve()
+
+ ENDDO
+
+ call mg_finalize()
+
+ end subroutine mg_main_mnh
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+end module mode_mg_main_mnh
+
+program mg_main_mnh_all
+
+ use mode_mg_main_mnh
+
+ call mg_main_mnh()
+
+end program mg_main_mnh_all
diff --git a/tensorproductmultigrid_Source/mode_mg.f90 b/tensorproductmultigrid_Source/mode_mg.f90
index 1ee7c5ad6..a607eecf2 100644
--- a/tensorproductmultigrid_Source/mode_mg.f90
+++ b/tensorproductmultigrid_Source/mode_mg.f90
@@ -8,7 +8,8 @@ module mode_mg
use solver
use parameters
use timer
- use mpi
+ !use mpi
+ use modd_mpif
use profiles
use mode_mg_read_param
@@ -56,8 +57,13 @@ subroutine mg_init_mnh()
implicit none
+logical :: gisinit
+
! Initialise MPI ...
- call mpi_init(ierr)
+ call mpi_initialized(gisinit, ierr )
+ if (.not. gisinit ) then
+ call mpi_init(ierr)
+ end if
! ... and pre initialise communication module
call comm_preinitialise()
@@ -137,11 +143,9 @@ implicit none
call create_scalar3d(MPI_COMM_HORIZ,grid_param,comm_param%halo_size,u)
call create_scalar3d(MPI_COMM_HORIZ,grid_param,comm_param%halo_size,b)
call create_scalar3d(MPI_COMM_HORIZ,grid_param,comm_param%halo_size,r)
- call initialise_rhs(grid_param,model_param,b)
-#ifdef TESTCONVERGENCE
- call create_scalar3d(MPI_COMM_HORIZ,grid_param,comm_param%halo_size,uerror)
- call analytical_solution(grid_param,uerror)
-#endif
+
+!!$ call mg_initialise_rhs_mnh()
+
call finish_timer(t_initialise)
if (i_am_master_mpi) then
write(STDOUT,*) ''
diff --git a/tensorproductmultigrid_Source/mode_mg_read_param.f90 b/tensorproductmultigrid_Source/mode_mg_read_param.f90
index 8e4887ad0..62b8ac478 100644
--- a/tensorproductmultigrid_Source/mode_mg_read_param.f90
+++ b/tensorproductmultigrid_Source/mode_mg_read_param.f90
@@ -36,7 +36,8 @@ subroutine read_general_parameters(filename,savefields_out)
use parameters
use communication
use messages
- use mpi
+ !use mpi
+ use modd_mpif
implicit none
character(*), intent(in) :: filename
logical, intent(out) :: savefields_out
@@ -67,7 +68,8 @@ subroutine read_solver_parameters(filename,solver_param_out)
use parameters
use communication
use messages
- use mpi
+ !use mpi
+ use modd_mpif
implicit none
character(*), intent(in) :: filename
type(solver_parameters), intent(out) :: solver_param_out
@@ -112,7 +114,8 @@ subroutine read_grid_parameters(filename,grid_param)
use datatypes
use communication
use messages
- use mpi
+ !use mpi
+ use modd_mpif
implicit none
character(*), intent(in) :: filename
type(grid_parameters), intent(out) :: grid_param
@@ -133,8 +136,8 @@ subroutine read_grid_parameters(filename,grid_param)
write(STDOUT,'("Grid parameters")')
write(STDOUT,'(" n = ",I15)') n
write(STDOUT,'(" nz = ",I15)') nz
- write(STDOUT,'(" L = ",F8.4)') L
- write(STDOUT,'(" H = ",F8.4)') H
+ write(STDOUT,'(" L = ",E15.6)') L
+ write(STDOUT,'(" H = ",E15.6)') H
if (vertbc == VERTBC_DIRICHLET) then
write(STDOUT,'(" vertbc = DIRICHLET")')
else if (vertbc == VERTBC_NEUMANN) then
@@ -170,7 +173,8 @@ subroutine read_comm_parameters(filename,comm_param)
use parameters
use communication
use messages
- use mpi
+ !use mpi
+ use modd_mpif
implicit none
character(*), intent(in) :: filename
type(comm_parameters), intent(out) :: comm_param
@@ -208,7 +212,8 @@ subroutine read_model_parameters(filename,model_param)
use discretisation
use communication
use messages
- use mpi
+ !use mpi
+ use modd_mpif
implicit none
character(*), intent(in) :: filename
type(model_parameters), intent(out) :: model_param
@@ -245,7 +250,8 @@ subroutine read_smoother_parameters(filename,smoother_param)
use discretisation
use communication
use messages
- use mpi
+ !use mpi
+ use modd_mpif
implicit none
character(*), intent(in) :: filename
type(smoother_parameters), intent(out) :: smoother_param
@@ -309,7 +315,8 @@ subroutine read_multigrid_parameters(filename,mg_param)
use multigrid
use communication
use messages
- use mpi
+ !use mpi
+ use modd_mpif
implicit none
character(*), intent(in) :: filename
type(mg_parameters), intent(out) :: mg_param
@@ -406,7 +413,8 @@ subroutine read_conjugategradient_parameters(filename,cg_param)
use conjugategradient
use communication
use messages
- use mpi
+ !use mpi
+ use modd_mpif
implicit none
character(*), intent(in) :: filename
type(cg_parameters), intent(out) :: cg_param
diff --git a/tensorproductmultigrid_Source/multigrid.f90 b/tensorproductmultigrid_Source/multigrid.f90
index f3017ef94..2ac7232f0 100644
--- a/tensorproductmultigrid_Source/multigrid.f90
+++ b/tensorproductmultigrid_Source/multigrid.f90
@@ -36,7 +36,8 @@
!==================================================================
module multigrid
- use mpi
+ !use mpi
+ use modd_mpif
use parameters
use datatypes
use discretisation
@@ -51,6 +52,7 @@ module multigrid
public::mg_parameters
public::mg_initialise
public::mg_finalise
+public::mg_solve_mnh
public::mg_solve
public::measurehaloswap
public::REST_CELLAVERAGE
diff --git a/tensorproductmultigrid_Source/profiles.f90 b/tensorproductmultigrid_Source/profiles.f90
index 7baf05ed1..f54093c07 100644
--- a/tensorproductmultigrid_Source/profiles.f90
+++ b/tensorproductmultigrid_Source/profiles.f90
@@ -42,6 +42,7 @@ module profiles
implicit none
+ public::initialise_rhs_mnh
public::initialise_rhs
public::analytical_solution
@@ -50,6 +51,47 @@ private
!==================================================================
! Initialise RHS vector
+!==================================================================
+ subroutine initialise_rhs_mnh(grid_param,model_param,b,xcoef)
+ implicit none
+ type(grid_parameters), intent(in) :: grid_param
+ type(model_parameters), intent(in) :: model_param
+ type(scalar3d), intent(inout) :: b
+ integer :: ix, iy, iz, ix_min, ix_max, iy_min, iy_max
+ real(kind=rl) :: x, y, z
+ real(kind=rl) :: rho, sigma, theta, phi, r, b_low, b_up, pi
+ real(kind=rl) :: xcoef
+
+ ix_min = b%ix_min
+ ix_max = b%ix_max
+ iy_min = b%iy_min
+ iy_max = b%iy_max
+ b_low = 1.0_rl+0.25*b%grid_param%H
+ b_up = 1.0_rl+0.75*b%grid_param%H
+ pi = 4.0_rl*atan2(1.0_rl,1.0_rl)
+ ! Initialise RHS
+ do ix=ix_min, ix_max
+ do iy=iy_min, iy_max
+ do iz=1,b%grid_param%nz
+ x = 1.0_rl*((ix-0.5_rl)/(1.0_rl*b%grid_param%n))
+ y = 1.0_rl*((iy-0.5_rl)/(1.0_rl*b%grid_param%n))
+ z = 1.0_rl*((iz-0.5_rl)/(1.0_rl*b%grid_param%nz))
+
+ b%s(iz,iy-iy_min+1,ix-ix_min+1) = 0.0_rl
+
+ if ( ( x .ge. 0.1_rl ) .and. ( x .le. 0.4_rl ) &
+ .and. (y .ge. 0.3_rl ) .and. ( y .le. 0.6_rl ) &
+ .and. (z .ge. 0.2_rl ) .and. ( z .le. 0.7_rl ) ) &
+ then
+ b%s(iz,iy-iy_min+1,ix-ix_min+1) = 1.0_rl * xcoef
+ end if
+
+ end do
+ end do
+ end do
+ end subroutine initialise_rhs_mnh
+!==================================================================
+! Initialise RHS vector
!==================================================================
subroutine initialise_rhs(grid_param,model_param,b)
implicit none
@@ -132,7 +174,6 @@ private
end do
end do
end subroutine initialise_rhs
-
!==================================================================
! Exact solution for test problem
! u(x,y,z) = x*(1-x)*y*(1-y)*z*(1-z)
diff --git a/tensorproductmultigrid_Source/timer.f90 b/tensorproductmultigrid_Source/timer.f90
index 32c980f23..54db8c541 100644
--- a/tensorproductmultigrid_Source/timer.f90
+++ b/tensorproductmultigrid_Source/timer.f90
@@ -36,7 +36,8 @@
module timer
- use mpi
+ !use mpi
+ use modd_mpif
use parameters
implicit none
diff --git a/tridz.f90 b/tridz.f90
index 0073ac7c9..442ce8974 100644
--- a/tridz.f90
+++ b/tridz.f90
@@ -201,6 +201,8 @@ USE MODE_REPRO_SUM
!JUAN
USE MODE_MPPDB
!
+USE mode_mg_main_mnh
+!
IMPLICIT NONE
!
!
@@ -533,6 +535,7 @@ PCF_ZS(IIB:IIE,IJB:IJE,IKB-1) = 0.5 * ( PRHO_ZS(IIB:IIE,IJB:IJE,IKB-1) + PRHO_Z
PAF_ZS(IIB:IIE,IJB:IJE,IKB-1) = 0.0
PCF_ZS(IIB:IIE,IJB:IJE,IKE+1) = 0.0
+call mg_main_mnh_init()
!------------------------------------------------------------------------------
!* 7. COMPUTE THE DIAGONAL ELEMENTS B OF THE MATRIX
--
GitLab