diff --git a/tensorproductmultigrid_Source/communication.f90 b/tensorproductmultigrid_Source/communication.f90
index 2cbb3ee82aaa751e0a5d2c4c8fcd7546bbffb062..5c6a85ab0abf684bed838628d9fc304b24cf7b40 100644
--- a/tensorproductmultigrid_Source/communication.f90
+++ b/tensorproductmultigrid_Source/communication.f90
@@ -582,6 +582,9 @@ contains
integer :: n, ix_min,ix_max,iy_min,iy_max
! Update Real Boundary for Newman case u(0) = u(1) , etc ...
+
+ !return
+
n = a%grid_param%n
ix_min = a%ix_min
ix_max = a%ix_max
diff --git a/tensorproductmultigrid_Source/mode_mg_read_param.f90 b/tensorproductmultigrid_Source/mode_mg_read_param.f90
index 713eb21fef8d9c349382b5d36ab8d51205edcc97..87b1ae7f2b2d8e4c162e386a9f6d29ddb5015a28 100644
--- a/tensorproductmultigrid_Source/mode_mg_read_param.f90
+++ b/tensorproductmultigrid_Source/mode_mg_read_param.f90
@@ -79,13 +79,14 @@ subroutine read_solver_parameters(filename,solver_param_out)
integer, parameter :: file_id = 16
integer :: ierr
namelist /parameters_solver/ solvertype,resreduction, maxiter &
- , LUseO , LUseT
+ , LUseO , LUseT , LMean
if (i_am_master_mpi) then
open(file_id,file=filename)
read(file_id,NML=parameters_solver)
close(file_id)
write(STDOUT,NML=parameters_solver)
write(STDOUT,'("---------------------------------------------- ")')
+ write(STDOUT,'(" LMean = ",L8)') LMean
write(STDOUT,'(" LUseO = ",L8)') LUseO
write(STDOUT,'(" LUseT = ",L8)') LUseT
write(STDOUT,'("Solver parameters ")')
@@ -104,6 +105,7 @@ subroutine read_solver_parameters(filename,solver_param_out)
call mpi_bcast(solvertype,1,MPI_INTEGER,master_rank,MPI_COMM_WORLD,ierr)
call mpi_bcast(maxiter,1,MPI_INTEGER,master_rank,MPI_COMM_WORLD,ierr)
call mpi_bcast(resreduction,1,MPI_DOUBLE_PRECISION,master_rank,MPI_COMM_WORLD,ierr)
+ call mpi_bcast(LMean,1,MPI_LOGICAL,master_rank,MPI_COMM_WORLD,ierr)
call mpi_bcast(LUseO,1,MPI_LOGICAL,master_rank,MPI_COMM_WORLD,ierr)
call mpi_bcast(LUseT,1,MPI_LOGICAL,master_rank,MPI_COMM_WORLD,ierr)
solver_param_out%solvertype = solvertype
diff --git a/tensorproductmultigrid_Source/multigrid.f90 b/tensorproductmultigrid_Source/multigrid.f90
index 9438e254d79ff716acbc1d80b025786dd5e7fd57..ca0f97a64ba462ed78fc18f53bafe33f29569b94 100644
--- a/tensorproductmultigrid_Source/multigrid.f90
+++ b/tensorproductmultigrid_Source/multigrid.f90
@@ -1417,13 +1417,14 @@ contains
integer :: maxiter
integer :: n_gridpoints
integer :: iter, level, finelevel, splitlevel
- real(kind=rl) :: res_old, res_new, res_initial
+ real(kind=rl) :: res_old, res_new, res_initial , mean_initial , norm_initial
logical :: solverconverged = .false.
integer :: nlocalx, nlocaly
integer :: halo_size
type(time) :: t_prec, t_res, t_apply, t_l2norm, t_scalprod, t_mainloop
type(scalar3d) :: pp
type(scalar3d) :: q
+ type(scalar3d) :: zone
real(kind=rl) :: alpha, beta, pq, rz, rz_old
integer :: ierr
@@ -1446,17 +1447,44 @@ contains
call initialise_timer(t_res,"t_residual")
call initialise_timer(t_mainloop,"t_mainloop")
- ! Copy b to b(1)
- ! Copy usolution to u(1)
n_gridpoints = (nlocalx+2*halo_size) &
* (nlocaly+2*halo_size) &
* (usolution%grid_param%nz+2)
+ !
+ ! Init 1 vector = zone
+ call create_scalar3d(MPI_COMM_HORIZ,bRHS%grid_param,halo_size,zone)
+ if (LUseO) then
+ zone%s(:,:,:) = 0.0_rl
+ zone%s(1:zone%grid_param%nz,1:zone%icompy_max,1:zone%icompx_max) = 1.0_rl
+ end if
+ if (LUseT) then
+ zone%st(:,:,:) = 0.0_rl
+ zone%st(1:zone%icompx_max,1:zone%icompy_max,1:zone%grid_param%nz) = 1.0_rl
+ end if
+ ! Mean / Norm of B
+ call scalarprod_mnh(pproc,zone,zone, mean_initial )
+ call scalarprod_mnh(pproc,zone,bRHS, mean_initial )
+ norm_initial = l2norm_mnh(bRHS,.true.)
+ mean_initial = mean_initial / (( zone%grid_param%nz ) * ( zone%grid_param%n )**2)
+ norm_initial = mean_initial / norm_initial
+ if (LMean) then
+ ! b <- b -mean_initial * zone
+ if (LUseO) call daxpy(n_gridpoints,-mean_initial,zone%s,1,bRHS%s,1)
+ if (LUseT) call daxpy(n_gridpoints,-mean_initial,zone%st,1,bRHS%st,1)
+ call scalarprod_mnh(pproc,zone,bRHS, mean_initial )
+ endif
+ !
+ ! Copy b to b(1)
+ ! Copy usolution to u(1)
if (LUseO) call dcopy(n_gridpoints,bRHS%s,1,xb_mg(level,pproc)%s,1)
if (LUseT) call dcopy(n_gridpoints,bRHS%st,1,xb_mg(level,pproc)%st,1)
if (LUseO) call dcopy(n_gridpoints,usolution%s,1,xu_mg(level,pproc)%s,1)
if (LUseT) call dcopy(n_gridpoints,usolution%st,1,xu_mg(level,pproc)%st,1)
+
! Scale with volume of grid cells
call volscale_scalar3d_mnh(xb_mg(level,pproc),1)
+ call scalarprod_mnh(pproc,zone,xb_mg(level,pproc), mean_initial )
+
call start_timer(t_res)
call calculate_residual_mnh(level, pproc, &
xb_mg(level,pproc),xu_mg(level,pproc),xr_mg(level,pproc))
diff --git a/tensorproductmultigrid_Source/parameters.f90 b/tensorproductmultigrid_Source/parameters.f90
index a578dcf177616bcec8234c88097ab736edba2fa7..5195982fdeffc44390f2a2fa9fd116f6e0d76660 100644
--- a/tensorproductmultigrid_Source/parameters.f90
+++ b/tensorproductmultigrid_Source/parameters.f90
@@ -58,5 +58,7 @@ module parameters
! Use Original (K,J,I) or Transposed (I,J,K) array
logical :: LUseO = .false.
logical :: LUseT = .true. ! .false. ! .true.
+! Remove mean to B = Right Hand side
+ logical :: LMean = .false.
end module parameters
diff --git a/tensorproductmultigrid_Source/profiles.f90 b/tensorproductmultigrid_Source/profiles.f90
index f383dbaa7d3bfeb1072604e6b287db1db8dd2775..489210645307e67cfc4c21e87f668963db60cbf4 100644
--- a/tensorproductmultigrid_Source/profiles.f90
+++ b/tensorproductmultigrid_Source/profiles.f90
@@ -75,7 +75,8 @@ private
do iy=iy_min, iy_max
do iz=1,u%grid_param%nz
- u%s(iz,iy-iy_min+1,ix-ix_min+1) = 0.0_rl
+ if (LUseO) u%s(iz,iy-iy_min+1,ix-ix_min+1) = 0.0_rl
+ if (LUseT) u%st(ix-ix_min+1,iy-iy_min+1,iz) = 0.0_rl
end do
end do