Juan 24/05/2019: Add management LMean=F/T of substracting Mean value to B...

Juan 24/05/2019: Add management LMean=F/T  of substracting Mean value to B ,for idealised case MG_MAIN_MNH_ALL

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

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

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))

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

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