Skip to content
New issue

Have a question about this project? # for a free GitHub account to open an issue and contact its maintainers and the community.

#

By clicking “#”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? # to your account

Huge IBM change by Stefano #133

Merged
merged 1 commit into from
Feb 2, 2023
Merged

Huge IBM change by Stefano #133

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions Solver/configure
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -51,6 +51,7 @@ TEST_CASES="./Euler/BoxAroundCircle \
./NavierStokes/CylinderDucros \
./NavierStokes/CylinderDifferentOrders \
./NavierStokes/CylinderLES \
./NavierStokes/IBM_Cylinder \
./NavierStokes/DualTimeStepping \
./NavierStokes/EntropyConservingTest \
./NavierStokes/EnergyConservingTest \
Expand Down
95 changes: 79 additions & 16 deletions Solver/src/NavierStokesSolver/SpatialDiscretization.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -357,8 +357,8 @@ subroutine TimeDerivative_ComputeQDot( mesh , particles, t)
! Local variables
! ---------------
!
integer :: eID , i, j, k, ierr, fID, iFace, iEl
real(kind=RP) :: mu_smag, delta, Source(NCONS)
integer :: eID , i, j, k, ierr, fID, iFace, iEl, iP
real(kind=RP) :: mu_smag, delta, Source(NCONS), TurbulentSource(NCONS)
!
! ***********************************************
! Compute the viscosity at the elements and faces
Expand Down Expand Up @@ -602,19 +602,40 @@ subroutine TimeDerivative_ComputeQDot( mesh , particles, t)
! *********************
! Add IBM source term
! *********************
if( mesh% IBM% active .and. .not. mesh% IBM% semiImplicit ) then
!$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)
if( e% isInsideBody(i,j,k) ) then
call mesh% IBM% SourceTerm( eID = eID, Q = e % storage % Q(:,i,j,k), Source = Source )
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
if( mesh% IBM% active ) then
if( .not. mesh% IBM% semiImplicit .and. t .gt. 0.0_RP ) 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
call mesh% IBM% SourceTerm( eID = eID, Q = e % storage % Q(:,i,j,k), Source = Source, wallfunction = .false. )
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
if( mesh% IBM% Wallfunction ) then
!$omp single
call mesh% IBM% GetBandRegionStates( mesh% elements )
!$omp end single
!$omp do schedule(runtime) private(i,j,k,TurbulentSource)
do iP = 1, mesh% IBM% NumOfForcingPoints
associate( e => mesh% elements(mesh% IBM% ImagePoints(iP)% element_index), &
e_in => mesh% elements(mesh% IBM% ImagePoints(iP)% element_in) )
i = mesh% IBM% ImagePoints(iP)% local_position(1)
j = mesh% IBM% ImagePoints(iP)% local_position(2)
k = mesh% IBM% ImagePoints(iP)% local_position(3)
call mesh % IBM % SourceTermTurbulence( mesh% IBM% ImagePoints(iP), e% storage% Q(:,i,j,k), &
e% geom% normal(:,i,j,k), e% geom% dWall(i,j,k), &
e% STL(i,j,k), TurbulentSource )
e% storage% QDot(:,i,j,k) = e % storage % QDot(:,i,j,k) + TurbulentSource
end associate
end do
!$omp end do
end if
end if
end if

end subroutine TimeDerivative_ComputeQDot
Expand Down Expand Up @@ -671,6 +692,7 @@ subroutine compute_viscosity_at_faces(no_of_faces, no_of_sides, face_ids, mesh)




end subroutine compute_viscosity_at_faces
!
!////////////////////////////////////////////////////////////////////////
Expand All @@ -690,8 +712,9 @@ subroutine TimeDerivative_ComputeQDotIsolated( mesh , t )
! Local variables
! ---------------
!
integer :: eID , i, j, k, fID
integer :: eID , i, j, k, fID, iP
procedure(UserDefinedSourceTermNS_f) :: UserDefinedSourceTermNS
real(kind=rp) :: Source(NCONS), TurbulentSource(NCONS)
!
! ****************
! Volume integrals
Expand Down Expand Up @@ -745,6 +768,46 @@ subroutine TimeDerivative_ComputeQDotIsolated( mesh , t )
end do
!$omp end do

!
! *********************
! Add IBM source term
! *********************

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
call mesh% IBM% SourceTerm( eID = eID, Q = e % storage % Q(:,i,j,k), Source = Source, wallfunction = .false. )
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
if( mesh% IBM% Wallfunction ) then
!$omp single
call mesh% IBM% GetBandRegionStates( mesh% elements )
!$omp end single
!$omp do schedule(runtime) private(i,j,k,TurbulentSource)
do iP = 1, mesh% IBM% NumOfForcingPoints
associate( e => mesh% elements(mesh% IBM% ImagePoints(iP)% element_index) )
i = mesh% IBM% ImagePoints(iP)% local_position(1)
j = mesh% IBM% ImagePoints(iP)% local_position(2)
k = mesh% IBM% ImagePoints(iP)% local_position(3)
call mesh % IBM % SourceTermTurbulence( mesh% IBM% ImagePoints(iP), e% storage% Q(:,i,j,k), &
e% geom% normal(:,i,j,k), e% geom% dWall(i,j,k), &
e% STL(i,j,k), TurbulentSource )
e% storage% QDot(:,i,j,k) = e % storage % QDot(:,i,j,k) + TurbulentSource
end associate
end do
!$omp end do
end if
end if
endif

end subroutine TimeDerivative_ComputeQDotIsolated
!
!///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
Expand Down Expand Up @@ -1261,4 +1324,4 @@ SUBROUTINE computeBoundaryFlux(f, time, mesh)

end subroutine computeBoundaryFlux

end module SpatialDiscretization
end module SpatialDiscretization
93 changes: 77 additions & 16 deletions Solver/src/NavierStokesSolverRANS/SpatialDiscretization.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -386,8 +386,9 @@ subroutine TimeDerivative_ComputeQDot( mesh , particles, t)
! Local variables
! ---------------
!
integer :: eID , i, j, k, ierr, fID, iFace, iEl
real(kind=RP) :: mu_smag, delta, mu_t, eta, kinematic_viscocity, mu_dim, Source(NCONS)
integer :: eID , i, j, k, ierr, fID, iFace, iEl, iP
real(kind=RP) :: mu_smag, delta, mu_t, eta, kinematic_viscocity, mu_dim, &
Source(NCONS), TurbulentSource(NCONS)
logical :: isfirst = .TRUE.
!
! ***********************************************
Expand Down Expand Up @@ -649,19 +650,39 @@ subroutine TimeDerivative_ComputeQDot( mesh , particles, t)
! *********************
! Add IBM source term
! *********************
if( mesh% IBM% active .and. .not. mesh% IBM% semiImplicit ) then
!$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)
if( e% isInsideBody(i,j,k) ) then
call mesh% IBM% SourceTerm( eID = eID, Q = e % storage % Q(:,i,j,k), Source = Source )
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
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
call mesh% IBM% SourceTerm( eID = eID, Q = e % storage % Q(:,i,j,k), Source = Source, wallfunction = .false. )
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
if( mesh% IBM% Wallfunction ) then
!$omp single
call mesh% IBM% GetBandRegionStates( mesh% elements )
!$omp end single
!$omp do schedule(runtime) private(i,j,k,TurbulentSource)
do iP = 1, mesh% IBM% NumOfForcingPoints
associate( e => mesh% elements(mesh% IBM% ImagePoints(iP)% element_index) )
i = mesh% IBM% ImagePoints(iP)% local_position(1)
j = mesh% IBM% ImagePoints(iP)% local_position(2)
k = mesh% IBM% ImagePoints(iP)% local_position(3)
call mesh % IBM % SourceTermTurbulence( mesh% IBM% ImagePoints(iP), e% storage% Q(:,i,j,k), &
e% geom% normal(:,i,j,k), e% geom% dWall(i,j,k), &
e% STL(i,j,k), TurbulentSource )
e% storage% QDot(:,i,j,k) = e % storage % QDot(:,i,j,k) + TurbulentSource
end associate
end do
!$omp end do
end if
end if
end if

end subroutine TimeDerivative_ComputeQDot
Expand Down Expand Up @@ -733,8 +754,9 @@ subroutine TimeDerivative_ComputeQDotIsolated( mesh , t )
! Local variables
! ---------------
!
integer :: eID , i, j, k, fID
integer :: eID , i, j, k, fID, iP
procedure(UserDefinedSourceTermNS_f) :: UserDefinedSourceTermNS
real(kind=rp) :: Source(NCONS), TurbulentSource(NCONS)
!
! ****************
! Volume integrals
Expand Down Expand Up @@ -785,6 +807,45 @@ subroutine TimeDerivative_ComputeQDotIsolated( mesh , t )
end do
!$omp end do

!
! *********************
! Add IBM source term
! *********************
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
call mesh% IBM% SourceTerm( eID = eID, Q = e % storage % Q(:,i,j,k), Source = Source, wallfunction = .false. )
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
if( mesh% IBM% Wallfunction ) then
!$omp single
call mesh% IBM% GetBandRegionStates( mesh% elements )
!$omp end single
!$omp do schedule(runtime) private(i,j,k,TurbulentSource)
do iP = 1, mesh% IBM% NumOfForcingPoints
associate( e => mesh% elements(mesh% IBM% ImagePoints(iP)% element_index) )
i = mesh% IBM% ImagePoints(iP)% local_position(1)
j = mesh% IBM% ImagePoints(iP)% local_position(2)
k = mesh% IBM% ImagePoints(iP)% local_position(3)
call mesh % IBM % SourceTermTurbulence( mesh% IBM% ImagePoints(iP), e% storage% Q(:,i,j,k), &
e% geom% normal(:,i,j,k), e% geom% dWall(i,j,k), &
e% STL(i,j,k), TurbulentSource )
e% storage% QDot(:,i,j,k) = e % storage % QDot(:,i,j,k) + TurbulentSource
end associate
end do
!$omp end do
end if
end if
end if

end subroutine TimeDerivative_ComputeQDotIsolated
!
!///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
Expand Down Expand Up @@ -1320,4 +1381,4 @@ SUBROUTINE computeBoundaryFlux_NSSA(f, time, mesh)

END SUBROUTINE computeBoundaryFlux_NSSA

end module SpatialDiscretization
end module SpatialDiscretization
5 changes: 1 addition & 4 deletions Solver/src/addons/horses2tecplot/Storage.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -249,9 +249,6 @@ subroutine Mesh_ReadSolution(self,solutionName)
self % hasSensor = .false.
self % hasTimeDeriv = .true.

! TODO?
case (IBM_MESH)

case (SOLUTION_AND_SENSOR_FILE)
dimensionsSize = 4
self % hasGradients = .false.
Expand Down Expand Up @@ -551,4 +548,4 @@ Subroutine getNVARS(originalDim, useFlowEq)
end select

End Subroutine getNVARS
end module Storage
end module Storage
Loading