diff --git a/tensorproductmultigrid_Source/communication.f90 b/tensorproductmultigrid_Source/communication.f90
index f0856420bb9b7209750330fab75ec3691882e3af..c0b82d05c6752be4fad908f66468b0527f4aca82 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 b58adae6bb572e192aee480f35b94b2137388c20..eaa98b3a05f564ec7d48c3fbe11f7196e3379c8a 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 b0bfbecc3ba82cb8dd7efd828e55af5c9e901219..4705e12019fa94b3924ec61bfc7184d14024fcff 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 2a45884d2504aed2fbedc5ff08b78621e984c57f..35f8b0c89d0240ed055df7bfa6da744829042009 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 4ade93172d586a6ab3e890f46f0aadb279cd0453..ce0fd76c35e8517fa73c433258e2b3f2d23fe0e2 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 4ff09112121ed11aafe1e03b80aa20b943594d27..0fef587424479d8c9b93dd5a81b6723203962f7e 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 0000000000000000000000000000000000000000..e99d66fd516a6ec29587387003956aad16ab509b
--- /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 1ee7c5ad60c10ae82c26f66890f70b7e2203c1f9..a607eecf284b6c06672b5fd194bfc6268960d581 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 8e4887ad01b98d0214413f76b746272c250ea3a3..62b8ac478da3e4c4a05e4a2e5d413be59c480036 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 f3017ef94447c1a9f778a1baa9ac277eb97ac016..2ac7232f0dd01590ab5f44645ec548e209942e15 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 7baf05ed192b5f677b4bbf031ea5c189494dc915..f54093c076fe45b89e2efcea94f4fb90bc8eab4b 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 32c980f236a898c3f9a45d048ca2851b76e818a9..54db8c541451934a7da6cda7e4f71c1fce824f36 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 0073ac7c93c35b9a7fed835d7e52c4dca1a6d92a..442ce8974f672afc37a7e1906b22e0e171b8fd04 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