Merge pull request #320 from xcompact3d/319-redesign-computation-of-r…

…homin-in-lmn-poisson-solver-to-make-it-explicit

Refactor LMN/Poisson solver to compute rho0 explicitly

src/navier.f90

@@ -58,7 +58,7 @@ SUBROUTINE solve_poisson(pp3, px1, py1, pz1, rho1, ux1, uy1, uz1, ep1, drho1, di
     nlock = 1 !! Corresponds to computing div(u*)
     converged = .FALSE.
     poissiter = 0
-    rho0 = one
+    call calc_rho0(rho1(:,:,:,1), rho0)
 #ifdef DOUBLE_PREC
     atol = 1.0e-14_mytype !! Absolute tolerance for Poisson solver
     rtol = 1.0e-14_mytype !! Relative tolerance for Poisson solver
@@ -85,8 +85,7 @@ SUBROUTINE solve_poisson(pp3, px1, py1, pz1, rho1, ux1, uy1, uz1, ep1, drho1, di
 
           IF (.NOT.converged) THEN
              !! Evaluate additional RHS terms
-             CALL calc_varcoeff_rhs(pp3(:,:,:,1), rho1, px1, py1, pz1, dv3, drho1, ep1, divu3, rho0, &
-                  poissiter)
+             CALL calc_varcoeff_rhs(pp3(:,:,:,1), rho1, px1, py1, pz1, dv3, drho1, ep1, divu3, rho0)
           ENDIF
 
        ENDIF
@@ -1172,12 +1171,41 @@ SUBROUTINE test_varcoeff(converged, divup3norm, pp3, dv3, atol, rtol, poissiter)
   ENDSUBROUTINE test_varcoeff
   !############################################################################
   !!
+  !! SUBROUTINE: calc_rho0
+  !! DESCRIPTION: Computes the density used to normalise the variable-coefficient Poisson equation.
+  subroutine calc_rho0(rho1, rho0)
+
+    use mpi
+
+    use param, only: ilmn, ivarcoeff
+
+    implicit none
+
+    real(mytype), dimension(xsize(1), xsize(2), xsize(3)), intent(in) :: rho1 ! The density field
+    real(mytype), intent(out) :: rho0                                         ! The normalising density (globally constant)
+
+    integer :: ierr
+
+    if (ilmn .and. ivarcoeff) then
+       !! Compute rho0
+       rho0 = minval(rho1(:,:,:))
+
+       call mpi_allreduce(MPI_IN_PLACE,rho0,1,real_type,mpi_min,mpi_comm_world,ierr)
+    else
+       ! Ensure a default value
+       rho0 = 1.0_mytype
+    end if
+
+  end subroutine calc_rho0
+  !############################################################################
+  !############################################################################
+  !!
   !!  SUBROUTINE: calc_varcoeff_rhs
   !!      AUTHOR: Paul Bartholomew
   !! DESCRIPTION: Computes RHS of the variable-coefficient Poisson solver
   !!
   !############################################################################
-  SUBROUTINE calc_varcoeff_rhs(pp3, rho1, px1, py1, pz1, dv3, drho1, ep1, divu3, rhomin, poissiter)
+  SUBROUTINE calc_varcoeff_rhs(pp3, rho1, px1, py1, pz1, dv3, drho1, ep1, divu3, rho0)
 
     USE MPI
 
@@ -1190,39 +1218,30 @@ SUBROUTINE calc_varcoeff_rhs(pp3, rho1, px1, py1, pz1, dv3, drho1, ep1, divu3, r
     IMPLICIT NONE
 
     !! INPUTS
-    INTEGER, INTENT(IN) :: poissiter
     REAL(mytype), INTENT(IN), DIMENSION(xsize(1), xsize(2), xsize(3)) :: px1, py1, pz1
     REAL(mytype), INTENT(IN), DIMENSION(xsize(1), xsize(2), xsize(3), nrhotime) :: rho1
     REAL(mytype), INTENT(IN), DIMENSION(xsize(1), xsize(2), xsize(3), ntime) :: drho1
     REAL(mytype), INTENT(IN), DIMENSION(xsize(1), xsize(2), xsize(3)) :: ep1
     REAL(mytype), INTENT(IN), DIMENSION(zsize(1), zsize(2), zsize(3)) :: divu3
     REAL(mytype), INTENT(IN), DIMENSION(ph1%zst(1):ph1%zen(1), ph1%zst(2):ph1%zen(2), nzmsize) :: dv3
-    REAL(mytype) :: rhomin
-    
+    real(mytype), intent(in) :: rho0
+
     !! OUTPUTS
     REAL(mytype), DIMENSION(ph1%zst(1):ph1%zen(1), ph1%zst(2):ph1%zen(2), nzmsize) :: pp3
 
     !! LOCALS
     INTEGER :: nlock, ierr
-
-
-    IF (poissiter.EQ.0) THEN
-       !! Compute rhomin
-       rhomin = MINVAL(rho1(:,:,:,1))
-
-       CALL MPI_ALLREDUCE(MPI_IN_PLACE,rhomin,1,real_type,MPI_MIN,MPI_COMM_WORLD,ierr)
-    ENDIF
 
-    ta1(:,:,:) = (one - rhomin / rho1(:,:,:,1)) * px1(:,:,:)
-    tb1(:,:,:) = (one - rhomin / rho1(:,:,:,1)) * py1(:,:,:)
-    tc1(:,:,:) = (one - rhomin / rho1(:,:,:,1)) * pz1(:,:,:)
+    ta1(:,:,:) = (one - rho0 / rho1(:,:,:,1)) * px1(:,:,:)
+    tb1(:,:,:) = (one - rho0 / rho1(:,:,:,1)) * py1(:,:,:)
+    tc1(:,:,:) = (one - rho0 / rho1(:,:,:,1)) * pz1(:,:,:)
 
     nlock = -1 !! Don't do any funny business with LMN
     call divergence(pp3,rho1,ta1,tb1,tc1,ep1,drho1,divu3,nlock)
 
     !! lapl(p) = div((1 - rhomin/rho) grad(p)) + rhomin(div(u*) - div(u))
     !! dv3 contains div(u*) - div(u)
-    pp3(:,:,:) = pp3(:,:,:) + rhomin * dv3(:,:,:)
+    pp3(:,:,:) = pp3(:,:,:) + rho0 * dv3(:,:,:)
 
   ENDSUBROUTINE calc_varcoeff_rhs
   !############################################################################