Furhter developents made in Delft.

Allmake

@@ -15,8 +15,11 @@
 # Make 'cappuccino' CFD solver:
 make execut=bin/cappuccino makefile_dir=src/cappuccino LDFLAGS=-O2\ -Wall LDLIBS=-llapack 
 #
+# Make Field Operations test:
+make execut=test/testFieldOperations/testFieldOperations makefile_dir=test/testFieldOperations LDFLAGS=-O2\ -Wall LDLIBS=-llapack
+#
 # Make Poisson equation solver:
-make execut=bin/poisson makefile_dir=applications/Poisson LDFLAGS=-O2\ -Wall LDLIBS=-llapack
+#make execut=bin/poisson makefile_dir=applications/Poisson LDFLAGS=-O2\ -Wall LDLIBS=-llapack
 #
 # Make Steady Heat Equation solver:
-make execut=bin/steadyHeatEq makefile_dir=applications/steadyHeatEq LDFLAGS=-O2\ -Wall LDLIBS=-llapack
+#make execut=bin/steadyHeatEq makefile_dir=applications/steadyHeatEq LDFLAGS=-O2\ -Wall LDLIBS=-llapack

src-par/CourantNo.h → src-par/CourantNo.f90

@@ -1,15 +1,27 @@
+subroutine CourantNo
 !
-! Calculate and output the mean and maximum Courant Numbers.
+! Purpose: 
+!   Calculate and output the mean and maximum Courant Numbers.
 !
+  use types
+  use parameters, only: CoNum,meanCoNum, CoNumFixValue, CoNumFix, timestep, ltransient, itime, time, myid
+  use geometry, only: numCells, numInnerFaces, owner, neighbour, numBoundaries, bctype, nfaces, startFace, Vol
+  use sparse_matrix, only: res
+  use variables, only: flmass
 
- if (ltransient) then
- CoNum = 0.0_dp
- meanCoNum = 0.0_dp
+  implicit none
 
- res = 0.0_dp
+  integer :: i, ijp, ijn, inp, ib, iface
+  real(dp):: suma,dt
+
+if (ltransient) then
+  CoNum = 0.0_dp
+  meanCoNum = 0.0_dp
+
+  res = 0.0_dp
 
   !
-  ! Suface sum of magnitude (i.e. absolute value) of mass flux phi, over inner faces only (which includes o- c- faces)
+  ! Suface sum of magnitude (i.e. absolute value) of mass flux phi, over inner faces only
   !
 
   ! Inner faces                                               
@@ -20,12 +32,6 @@
     res(ijn) = res(ijn)+abs(flmass(i))                                                                                                               
   enddo     
 
-  ! Faces on processor boundaries
-  !do i=1,npro
-  !  ijp = owner( iProcFacesStart + i )
-  !  res(ijp) = res(ijp)+abs(fmpro(i))
-  !end do
-
   ! Boundary faces
   do ib=1,numBoundaries
 
@@ -55,6 +61,7 @@
 
   enddo
 
+
   ! Accumulate by looping trough cells
   suma = 0.0_dp
 
@@ -75,16 +82,13 @@
   call global_sum(suma)
 
   CoNum = 0.5*CoNum*timestep
+  meanCoNum = 0.5*meanCoNum/suma*timestep
 
   ! Find global maximum Courant number in the whole field.
   call global_max(CoNum)
 
-
-  ! Now use it to calculate mean Courant number
-  meanCoNum = 0.5*meanCoNum/suma*timestep
-
   !// If we keep the value of Courant Number fixed
-  if( CoNumFix .and. itime.ne.itimes ) then
+  if( CoNumFix ) then
       dt = timestep
       timestep = CoNumFixValue/CoNum * timestep
 
@@ -94,14 +98,15 @@
 
   time = time + timestep
 
-if(myid .eq. 0) then
- write(6,*)
- write(6,'(a,i0,a,es10.3,a,es10.3)') "  Time step no. : ",ITIME," dt : ",timestep," Time = ",time
- write(6,*)
+  if(myid .eq. 0) then
+    write(6,'(a)') ' '
+    write(6,'(a,i0,2(a,es10.3))') "  Time step no. : ",itime," dt : ",timestep," Time = ",time
+    write(6,'(a)') ' '
+    write(6,'(2(a,es10.3))') "  Courant Number mean: ", meanCoNum," max: ", CoNum
+    write(6,'(a)') ' '
+  endif
 
- write(6,'(2(a,es10.3))') "  Courant Number mean: ", meanCoNum," max: ", CoNum
 endif
 
+end subroutine
 
-endif
-!// ************************************************************************* //

src-par/Makefile

@@ -40,7 +40,8 @@ LINEAR_SOLVER_FILES=\
     iccg.f90 \
     bicgstab.f90 \
     dpcg.f90 \
-    jacobi.f90
+    jacobi.f90 \
+    mgmres.f90
 
 TURBULENCE_FILES=\
     temperature.f90 \
@@ -61,23 +62,28 @@ CAFFA_FILES=\
     bpres.f90 \
     fieldManipulation.f90 \
     faceflux_velocity.f90 \
-    facefluxmass.f90 \
+    faceflux_mass.f90 \
     calcheatflux.f90 \
-    calcp-multiple_correction_SIMPLE.f90 \
+    calcp_simple.f90 \
+    calcp_piso.f90\
     calcstress.f90 \
     calc_strain_and_vorticity.f90 \
     calcuvw.f90 \
     calc_statistics.f90 \
+    CourantNo.f90 \
     updateVelocityAtBoundary.f90 \
     correct_turbulence.f90 \
     correct_turbulence_inlet.f90 \
+    continuityErrors.f90\
+    constant_mass_flow_forcing.f90 \
+    recirculate_flow.f90 \
     fvm_laplacian.f90 \
     init.f90 \
-    PISO_multiple_correction.f90 \
     readfiles.f90 \
     read_input.f90 \
     random_seed.f90 \
     time_shift.f90 \
+    vortexIdentification.f90\
     writefiles.f90 \
     write_restart_files.f90 \
     writehistory.f90 \

src-par/adjustMassFlow.f90

@@ -52,6 +52,7 @@ subroutine adjustMassFlow
 
   enddo
 
+  call global_sum(flowo)
 
   ! Correct mass flux to satisfy global mass conservation & add to source
   fac = flomas/(flowo+small)

src-par/allocate.f90

@@ -86,15 +86,12 @@ subroutine allocate_arrays
   allocate( p( numTotal ), stat=ierr) 
     if(ierr /= 0)write(*,*)"allocation error: p"
 
+  allocate(pp(numTotal),stat=ierr) 
+    if(ierr /= 0)write(*,*)"allocation error: pp"
+
   allocate( dPdxi( 3,numTotal ), stat=ierr) 
     if(ierr /= 0)write(*,*)"allocation error: dPdxi"
 
-  if ( simple ) then
-
-    allocate(pp(numTotal),stat=ierr) 
-      if(ierr /= 0)write(*,*)"allocation error: pp"
-
-  endif
 
   if ( piso ) then
 
@@ -193,6 +190,11 @@ subroutine allocate_arrays
     allocate( dCondxi( 3,numTotal ), stat=ierr) 
       if(ierr /= 0)write(*,*)"allocation error: dCondxi"
 
+    if (ltransient) then
+      allocate( con_aver( numTotal ), stat=ierr) 
+        if(ierr /= 0)write(*,*)"allocation error: tt_aver" 
+    endif
+
   endif
 
   ! Temperature variance and dissipation of temperature variance
@@ -276,6 +278,7 @@ subroutine allocate_arrays
   endif
 
 
+  ! Wall things
   allocate(  visw(nwal), stat=ierr) 
     if(ierr /= 0)write(*,*)"allocation error: visw"
 
@@ -285,6 +288,11 @@ subroutine allocate_arrays
   allocate(  tau(nwal), stat=ierr) 
     if(ierr /= 0)write(*,*)"allocation error: tau"
 
+  if (ltransient) then
+    allocate(  wss_aver(nwal), stat=ierr) 
+      if(ierr /= 0)write(*,*)"allocation error: wss_aver"
+  endif
+
   ! Turbulence production
   allocate( gen( numCells ), stat=ierr) 
     if(ierr /= 0)write(*,*)"allocation error: gen"  

src-par/bcin.f90

@@ -4,22 +4,23 @@ subroutine bcin
   use parameters
   use geometry
   use variables
+  use mpi
 
   implicit none
 
   integer :: i,ib,ini,ino,iface
   ! integer :: input_unit
-  integer :: totalNumInlets
+  ! integer :: totalNumInlets
   real(dp) :: uav,outare
-  real(dp) :: flowen, flowte, flowed
+  ! real(dp) :: flowen, flowte, flowed
   real(dp) :: are
 
 
   flomas = 0.0_dp 
-  flomom = 0.0_dp
-  flowen = 0.0_dp
-  flowte = 0.0_dp
-  flowed = 0.0_dp
+  ! flomom = 0.0_dp
+  ! flowen = 0.0_dp
+  ! flowte = 0.0_dp
+  ! flowed = 0.0_dp
 
 
   ! In the case we want to read recirculated profiles
@@ -55,10 +56,10 @@ subroutine bcin
           ! is faced inwards. That means their scalar product will be negative,
           ! so minus signs here is to turn net mass influx - flomas, into positive value.
           flomas = flomas - flmass(iface) 
-          flomom = flomom + abs(flmass(iface))*sqrt(u(ini)**2+v(ini)**2+w(ini)**2)
+          ! flomom = flomom + abs(flmass(iface))*sqrt(u(ini)**2+v(ini)**2+w(ini)**2)
           ! if(lcal(ien)) flowen = flowen + abs(flmass(iface)*t(ini))
-          flowte = flowte + abs(flmass(iface)*te(ini))
-          flowed = flowed + abs(flmass(iface)*ed(ini))
+          ! flowte = flowte + abs(flmass(iface)*te(ini))
+          ! flowed = flowed + abs(flmass(iface)*ed(ini))
 
 
         end do
@@ -67,6 +68,9 @@ subroutine bcin
 
     enddo
 
+    ! Correct turbulence at inlet for appropriate turbulence model
+    if(lturb) call correct_turbulence_inlet()
+
   endif  
 
   ! close(input_unit)
@@ -89,11 +93,10 @@ subroutine bcin
 
   ! Global sum - MPI communication
   call global_sum( flomas )
-  call global_sum( flomom )
-  if(lcal(ien)) call global_sum( flowen )
-  call global_sum( flowte )
-  call global_sum( flowed )
-
+  ! call global_sum( flomom )
+  ! if(lcal(ien)) call global_sum( flowen )
+  ! call global_sum( flowte )
+  ! call global_sum( flowed )
   call global_sum( outare )
 
   ! Average velocity at outlet boundary
@@ -104,7 +107,10 @@ subroutine bcin
     write ( *, '(a)' ) '  Inlet boundary condition information:'
     write ( *, '(a)' ) ' '
     write ( *, '(a,e12.6)' ) '  Mass inflow: ', flomas
-    write ( *, '(a,e12.6)' ) '  Momentum inflow: ', flomom
+    ! write ( *, '(a,e12.6)' ) '  Momentum inflow: ', flomom
+    ! write ( *, '(a,e12.6)' ) '  Temperature inflow: ', flowen  
+    ! write ( *, '(a,e12.6)' ) '  TKE inflow: ', flowte      
+    ! write ( *, '(a,e12.6)' ) '  Dissipation inflow: ', flowed    
   endif
 
 
@@ -138,22 +144,4 @@ subroutine bcin
 
   endif  
 
-  totalNumInlets = ninl 
-  call global_isum( totalNumInlets )
-
-
-  if ( totalNumInlets .eq. 0 ) then
-   ! No inflow into the domain: eg. natural convection case, etc.
-
-    flomas = 1.0_dp
-    flomom = 1.0_dp
-    flowen = 1.0_dp
-    flowte = 1.0_dp
-    flowed = 1.0_dp
-
-  endif
-
-  ! Correct turbulence at inlet for appropriate turbulence model
-  if(lturb) call correct_turbulence_inlet()
-
 end subroutine