Merge pull request #248 from loganoz/omarinoDev

Omarino dev

Solver/Makefile.in

@@ -336,6 +336,7 @@ mu: header mkdirs horseslibs-mu horses3d.mu
 ns: header mkdirs horseslibs-ns horses3d.ns
 nssa: header mkdirs horseslibs-nssa horses3d.nssa
 caa: header mkdirs horseslibs-caa horses3d.caa 
+ins: header mkdirs horseslibs-ins horses3d.ins
 
 horses3d.ns: 
 	(cd $(HOME_DIR)/src/NavierStokesSolver/ && $(MAKE) $(HORSES_NS_FLAGS) install) || exit 1; 
@@ -455,6 +456,24 @@ horseslibs-caa:
 		fi \
 	done
 
+horseslibs-ins:
+	@for subdir in $(HORSES_LIBS_GENERAL) ; do \
+		if [ "$$subdir" = "problemfile" ] ; \
+		then  \
+			(cd $(HOME_DIR)/src/libs/problemfile/ && $(MAKE) PLATFORM=$(PLATFORM) $(LIBFLAGS) DYNLIB_FLAG='$(DYNLIB_FLAG)' install) || exit 1 ; \
+	    	else \
+			(cd $(HOME_DIR)/src/libs/$$subdir && $(MAKE) $(LIBFLAGS) install) || exit 1; \
+		fi \
+	done
+	@for subdir in $(HORSES_LIBS_PARTICULAR) ; do \
+		if [ "$$subdir" = "problemfile" ] ; \
+		then  \
+			(cd $(HOME_DIR)/src/libs/problemfile/ && $(MAKE) PLATFORM=$(PLATFORM) $(LIBFLAGS) DYNLIB_FLAG='$(DYNLIB_FLAG)' install-ins) || exit 1 ; \
+	    	else \
+			(cd $(HOME_DIR)/src/libs/$$subdir && $(MAKE) $(LIBFLAGS) install-ins) || exit 1; \
+		fi \
+	done
+
 bench: all
 	(export PYTHONPATH=$(HOME_DIR)/../utils/PythonUtilities \
 	&& export PYTHONISDEFINED=YES \

Solver/src/IncompressibleNSSolver/SpatialDiscretization.f90

@@ -320,7 +320,8 @@ END SUBROUTINE ComputeTimeDerivativeIsolated
 !////////////////////////////////////////////////////////////////////////////////////
 !
       subroutine ComputeNSTimeDerivative( mesh , particles, t )
-         use SpongeClass, only: sponge
+         use SpongeClass, only: sponge, addSourceSponge
+         use ActuatorLine, only: farm, ForcesFarm
          implicit none
          type(HexMesh)              :: mesh
          type(Particles_t)          :: particles
@@ -333,6 +334,7 @@ subroutine ComputeNSTimeDerivative( mesh , particles, t )
 !
          integer     :: eID , i, j, k, ierr, fID
          real(kind=RP)  :: mu_smag, delta
+         real(kind=RP), dimension(NCONS)  :: Source
 !
 !        ****************
 !        Volume integrals
@@ -485,6 +487,7 @@ subroutine ComputeNSTimeDerivative( mesh , particles, t )
 !$omp do schedule(runtime) private(i,j,k)
             do eID = 1, mesh % no_of_elements
                associate ( e => mesh % elements(eID) )
+               e % storage % S_NS = 0.0_RP
                do k = 0, e % Nxyz(3)   ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1)
                   call UserDefinedSourceTermNS(e % geom % x(:,i,j,k), e % storage % Q(:,i,j,k), t, e % storage % S_NS(:,i,j,k), thermodynamics, dimensionless, refValues)
                end do                  ; end do                ; end do
@@ -493,7 +496,8 @@ subroutine ComputeNSTimeDerivative( mesh , particles, t )
 !$omp end do
 
 ! for the sponge, loops are in the internal subroutine as values are precalculated
-            call sponge % addSource(mesh)
+            call addSourceSponge(sponge,mesh)
+            call ForcesFarm(farm, mesh, t)
 
 !$omp do schedule(runtime) private(i,j,k)
             do eID = 1, mesh % no_of_elements
@@ -504,20 +508,43 @@ subroutine ComputeNSTimeDerivative( mesh , particles, t )
                end associate
             end do
 !$omp end do
-
 !
 !           ********************
 !           Add Particles source
 !           ********************
-            if (.not. mesh % child) then
-               if ( particles % active ) then             
-!$omp do schedule(runtime)
-                  do eID = 1, size(mesh % elements)
+            !if (.not. mesh % child) then
+               !if ( particles % active ) then             
+!!$omp do schedule(runtime)
+                  !do eID = 1, size(mesh % elements)
                   !   call particles % AddSource(mesh % elements(eID), t, thermodynamics, dimensionless, refValues)
+                  !end do
+!!$omp end do
+               !endif 
+            !end if
+!
+!        *********************
+!        Add IBM source term
+!        *********************
+! no wall function for INCNS
+         if( mesh% IBM% active ) then
+            if( .not. mesh% IBM% semiImplicit ) then 
+!$omp do schedule(runtime) private(i,j,k,Source)
+                  do eID = 1, mesh % no_of_elements  
+                     associate ( e => mesh % elements(eID) ) 
+                     do k = 0, e % Nxyz(3)   ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1)
+                        if( e% isInsideBody(i,j,k) ) then
+                           ! only without moving for now in INCNS
+                           if( .not. mesh% IBM% stl(e% STL(i,j,k))% move ) then 
+                              call mesh% IBM% SourceTerm( eID = eID, Q = e % storage % Q(:,i,j,k), Source = Source, wallfunction = .false. )
+                           end if 
+                           e % storage % QDot(:,i,j,k) = e % storage % QDot(:,i,j,k) + Source
+                        end if
+                     end do                  ; end do                ; end do
+                     end associate
                   end do
-!$omp end do
-               endif 
-            end if
+!$omp end do       
+            end if 
+         end if
 
       end subroutine ComputeNSTimeDerivative
 !
@@ -1040,4 +1067,4 @@ end subroutine DensityLimiter
 !
 !////////////////////////////////////////////////////////////////////////////////////////
 !
-end module SpatialDiscretization
+end module SpatialDiscretization

Solver/src/MultiphaseSolver/SpatialDiscretization.f90

@@ -601,7 +601,8 @@ END SUBROUTINE ComputeTimeDerivative
 !////////////////////////////////////////////////////////////////////////////////////
 !
       subroutine ComputeNSTimeDerivative( mesh , t )
-         use SpongeClass, only: sponge
+         use SpongeClass, only: sponge, addSourceSponge
+         use ActuatorLine, only: farm, ForcesFarm
          implicit none
          type(HexMesh)              :: mesh
          real(kind=RP)              :: t
@@ -614,6 +615,7 @@ subroutine ComputeNSTimeDerivative( mesh , t )
          integer     :: eID , i, j, k, ierr, fID
          real(kind=RP) :: sqrtRho, invSqrtRho
          real(kind=RP)  :: mu_smag, delta
+         real(kind=RP), dimension(NCONS)  :: Source
 !
 !        ****************
 !        Volume integrals
@@ -677,7 +679,7 @@ subroutine ComputeNSTimeDerivative( mesh , t )
 !
 !        *************************************************************************************
 !        Element without shared faces: Surface integrals, scaling of elements with Jacobian, 
-!                                      sqrt(rho), and add source terms
+!                                      sqrt(rho)
 !        *************************************************************************************
 ! 
 !$omp do schedule(runtime) private(i,j,k,sqrtRho,invSqrtRho)
@@ -698,36 +700,12 @@ subroutine ComputeNSTimeDerivative( mesh , t )
                e % storage % QDot(IMSQRHOU:IMSQRHOW,i,j,k) =   e % storage % QDot(IMSQRHOU:IMSQRHOW,i,j,k) & 
                                                              + sqrtRho * dimensionless % invFr2 * dimensionless % gravity_dir
 
-!
-!            + Add user defined source terms
-               call UserDefinedSourceTermNS(e % geom % x(:,i,j,k), e % storage % Q(:,i,j,k), t, e % storage % S_NS(:,i,j,k), thermodynamics, dimensionless, refValues, multiphase)
-               e % storage % QDot(:,i,j,k) = e % storage % QDot(:,i,j,k) + e % storage % S_NS(:,i,j,k) / [1.0_RP,sqrtRho,sqrtRho,sqrtRho,1.0_RP]
-
             end do         ; end do          ; end do 
 
             end associate 
          end do
 !$omp end do
 
-! Sponges
-!$omp do schedule(runtime)
-         do eID = 1, mesh % no_of_elements
-            associate ( e => mesh % elements(eID) )
-               e % storage % S_NS = 0.0_RP
-            end associate
-         enddo
-!$omp end do
-!for the sponge, loops are in the internal subroutine as values are precalculated
-         call sponge % addSource(mesh)
-!$omp do schedule(runtime) private(i,j,k)
-         do eID = 1, mesh % no_of_elements
-            associate ( e => mesh % elements(eID) )
-            do k = 0, e % Nxyz(3)   ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1)
-               e % storage % QDot(:,i,j,k) = e % storage % QDot(:,i,j,k) + e % storage % S_NS(:,i,j,k)
-            end do                  ; end do                ; end do
-            end associate
-         end do
-!$omp end do  
 !
 !        ******************************************
 !        Do the same for elements with shared faces
@@ -781,11 +759,6 @@ subroutine ComputeNSTimeDerivative( mesh , t )
                   e % storage % QDot(IMSQRHOU:IMSQRHOW,i,j,k) =   e % storage % QDot(IMSQRHOU:IMSQRHOW,i,j,k) & 
                                                                 + sqrtRho * dimensionless % invFr2 * dimensionless % gravity_dir
 
-!   
-!               + Add user defined source terms
-                  call UserDefinedSourceTermNS(e % geom % x(:,i,j,k), e % storage % Q(:,i,j,k), t, e % storage % S_NS(:,i,j,k), thermodynamics, dimensionless, refValues, multiphase)
-                  e % storage % QDot(:,i,j,k) = e % storage % QDot(:,i,j,k) + e % storage % S_NS(:,i,j,k) / [1.0_RP,sqrtRho,sqrtRho,sqrtRho,1.0_RP]
-
                end do         ; end do          ; end do 
 
 
@@ -804,6 +777,62 @@ subroutine ComputeNSTimeDerivative( mesh , t )
          end if
 #endif
 
+!           ***************
+!           Add source term
+!           ***************
+!$omp do schedule(runtime) private(i,j,k)
+            do eID = 1, mesh % no_of_elements
+               associate ( e => mesh % elements(eID) )
+               e % storage % S_NS = 0.0_RP
+               do k = 0, e % Nxyz(3)   ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1)
+                  InvSqrtRho = 1.0_RP / sqrt(e % storage % rho(i,j,k))
+                  call UserDefinedSourceTermNS(e % geom % x(:,i,j,k), e % storage % Q(:,i,j,k), t, e % storage % S_NS(:,i,j,k), thermodynamics, dimensionless, refValues, multiphase)
+                  ! scale UserDefinedSourceTerm momentum with sqrtRho
+                  e % storage % S_NS(:,i,j,k) = e % storage % S_NS(:,i,j,k) * [1.0_RP,InvSqrtRho,InvSqrtRho,InvSqrtRho,1.0_RP]
+               end do                  ; end do                ; end do
+               end associate
+            end do
+!$omp end do
+
+!for the sponge, loops are in the internal subroutine as values are precalculated
+!The scale with sqrtRho is done in the subroutines, not done againg here
+         call addSourceSponge(sponge,mesh)
+         call ForcesFarm(farm, mesh, t)
+
+! Add all the source terms
+!$omp do schedule(runtime) private(i,j,k)
+         do eID = 1, mesh % no_of_elements
+            associate ( e => mesh % elements(eID) )
+            do k = 0, e % Nxyz(3)   ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1)
+               e % storage % QDot(:,i,j,k) = e % storage % QDot(:,i,j,k) + e % storage % S_NS(:,i,j,k)
+            end do                  ; end do                ; end do
+            end associate
+         end do
+!$omp end do  
+!
+!        *********************
+!        Add IBM source term
+!        *********************
+! no wall function for MULTIPHASE
+         if( mesh% IBM% active ) then
+            if( .not. mesh% IBM% semiImplicit ) then 
+!$omp do schedule(runtime) private(i,j,k,Source)
+                  do eID = 1, mesh % no_of_elements  
+                     associate ( e => mesh % elements(eID) ) 
+                     do k = 0, e % Nxyz(3)   ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1)
+                        if( e% isInsideBody(i,j,k) ) then
+                           ! only without moving for now in MULTIPHASE
+                           if( .not. mesh% IBM% stl(e% STL(i,j,k))% move ) then 
+                              call mesh% IBM% SourceTerm( eID = eID, Q = e % storage % Q(:,i,j,k), Source = Source, wallfunction = .false. )
+                           end if 
+                           e % storage % QDot(:,i,j,k) = e % storage % QDot(:,i,j,k) + Source
+                        end if
+                     end do                  ; end do                ; end do
+                     end associate
+                  end do
+!$omp end do       
+            end if 
+         end if
 
       end subroutine ComputeNSTimeDerivative
 !
@@ -1628,4 +1657,4 @@ SUBROUTINE ComputeTimeDerivativeIsolated( mesh, particles, time, mode, HO_Elemen
 
       end subroutine ComputeTimeDerivativeIsolated
 
-end module SpatialDiscretization
+end module SpatialDiscretization

Solver/src/MultiphaseSolver/main.f90

@@ -60,10 +60,10 @@ PROGRAM HORSES3DMainMU
 !     ----------------------------------------------------------------------------------
 !
       if ( MPI_Process % doMPIAction ) then
-         CALL Main_Header("HORSES3D High-Order (DG) Spectral Element Parallel Incompressible Navier-Stokes Solver",__DATE__,__TIME__)
+         CALL Main_Header("HORSES3D High-Order (DG) Spectral Element Parallel Multiphase Inc. Navier-Stokes Solver",__DATE__,__TIME__)
 
       else
-         CALL Main_Header("HORSES3D High-Order (DG) Spectral Element Sequential Incompressible Navier-Stokes Solver",__DATE__,__TIME__)
+         CALL Main_Header("HORSES3D High-Order (DG) Spectral Element Sequential Multiphase Inc. Navier-Stokes Solver",__DATE__,__TIME__)
 
       end if
 

Solver/src/NavierStokesSolver/SpatialDiscretization.f90

@@ -361,8 +361,8 @@ END SUBROUTINE ComputeTimeDerivativeIsolated
       subroutine TimeDerivative_ComputeQDot( mesh , particles, t)
       use WallFunctionConnectivity
          use TripForceClass, only: randomTrip
-         use ActuatorLine, only: farm
-         use SpongeClass, only: sponge
+         use ActuatorLine, only: farm, ForcesFarm
+         use SpongeClass, only: sponge, addSourceSponge
          implicit none
          type(HexMesh)              :: mesh
          type(Particles_t)          :: particles
@@ -550,19 +550,17 @@ subroutine TimeDerivative_ComputeQDot( mesh , particles, t)
 !$omp do schedule(runtime) private(i,j,k)
             do eID = 1, mesh % no_of_elements
                associate ( e => mesh % elements(eID) )
-               ! the source term is reset to 0 each time Qdot is calculated to enable the possibility to add source terms to
-               ! different contributions and not accumulate each call
-               e % storage % S_NS = 0.0_RP
                do k = 0, e % Nxyz(3)   ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1)
+               ! All terms are calculated indepentenly and overwritten in case one gauss point has more than one contribution
                   call UserDefinedSourceTermNS(e % geom % x(:,i,j,k), e % storage % Q(:,i,j,k), t, e % storage % S_NS(:,i,j,k), thermodynamics, dimensionless, refValues)
                   call randomTrip % getTripSource( e % geom % x(:,i,j,k), e % storage % S_NS(:,i,j,k) )
-                  call farm % ForcesFarm(e % geom % x(:,i,j,k), e % storage % Q(:,i,j,k), e % storage % S_NS(:,i,j,k), t)
                end do                  ; end do                ; end do
                end associate
             end do
 !$omp end do
             ! for the sponge, loops are in the internal subroutine as values are precalculated
-            call sponge % addSource(mesh)
+            call addSourceSponge(sponge,mesh)
+            call ForcesFarm(farm, mesh, t)
 !
 !           Add Particles source
 !           ********************
@@ -675,7 +673,8 @@ end subroutine TimeDerivative_ComputeQDot
    subroutine TimeDerivative_ComputeQDotHO( mesh , particles, t)
       use WallFunctionConnectivity
          use TripForceClass, only: randomTrip
-         use ActuatorLine, only: farm
+         use ActuatorLine, only: farm, ForcesFarm
+         use SpongeClass, only: sponge, addSourceSponge
          implicit none
          type(HexMesh)              :: mesh
          type(Particles_t)          :: particles
@@ -865,17 +864,16 @@ subroutine TimeDerivative_ComputeQDotHO( mesh , particles, t)
 !$omp do schedule(runtime) private(i,j,k)
             do eID = 1, size(mesh % HO_Elements)
                associate ( e => mesh % elements(mesh % HO_Elements(eID)) )
-               ! the source term is reset to 0 each time Qdot is calculated to enable the possibility to add source terms to
-               ! different contributions and not accumulate each call
-               e % storage % S_NS = 0.0_RP
                do k = 0, e % Nxyz(3)   ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1)
                   call UserDefinedSourceTermNS(e % geom % x(:,i,j,k), e % storage % Q(:,i,j,k), t, e % storage % S_NS(:,i,j,k), thermodynamics, dimensionless, refValues)
                   call randomTrip % getTripSource( e % geom % x(:,i,j,k), e % storage % S_NS(:,i,j,k) )
-                  call farm % ForcesFarm(e % geom % x(:,i,j,k), e % storage % Q(:,i,j,k), e % storage % S_NS(:,i,j,k), t)
                end do                  ; end do                ; end do
                end associate
             end do
 !$omp end do
+
+            call addSourceSponge(sponge,mesh)
+            call ForcesFarm(farm, mesh, t)
 !
 !           Add Particles source
 !           ********************
@@ -1049,7 +1047,8 @@ end subroutine compute_viscosity_at_faces
 !     -------------------------------------------------------------------------------
       subroutine TimeDerivative_ComputeQDotIsolated( mesh , t )
          use TripForceClass, only: randomTrip
-         use ActuatorLine, only: farm
+         use ActuatorLine, only: farm, ForcesFarm
+         use SpongeClass, only: sponge, addSourceSponge
          implicit none
          type(HexMesh)              :: mesh
          real(kind=RP)              :: t
@@ -1093,15 +1092,15 @@ subroutine TimeDerivative_ComputeQDotIsolated( mesh , t )
 !$omp do schedule(runtime) private(i,j,k)
             do eID = 1, mesh % no_of_elements
                associate ( e => mesh % elements(eID) )
-               e % storage % S_NS = 0.0_RP
                do k = 0, e % Nxyz(3)   ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1)
                   call UserDefinedSourceTermNS(e % geom % x(:,i,j,k), e % storage % Q(:,i,j,k), t, e % storage % S_NS(:,i,j,k), thermodynamics, dimensionless, refValues)
                   call randomTrip % getTripSource( e % geom % x(:,i,j,k), e % storage % S_NS(:,i,j,k) )
-                  call farm % ForcesFarm(e % geom % x(:,i,j,k), e % storage % Q(:,i,j,k), e % storage % S_NS(:,i,j,k), t)
                end do                  ; end do                ; end do
                end associate
             end do
 !$omp end do
+            call ForcesFarm(farm, mesh, t)
+            call addSourceSponge(sponge,mesh)
          end if
 
 !$omp do schedule(runtime) private(i,j,k)

Solver/src/libs/discretization/Makefile

@@ -222,6 +222,12 @@ install-caa: lib$(LIB)_caa.a
 	@echo
 	@echo
 
+install-ins: lib$(LIB)_ins.a
+	cp -p lib$(LIB)_ins.a $(INSTALL_DIR)/lib/
+	cp -p ./include_ins/*.mod $(INSTALL_DIR)/include_ins/
+	@echo
+	@echo
+
 header: FORCE
 	@echo
 	@echo "================================"