Skip to content

293 acoustics save lamb vector from flow simulations and read it from acoustics solver #295

New issue

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

Sign up for GitHub

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

Already on GitHub? Sign in to your account

Merged
albertjimenezramos merged 9 commits into from
Feb 11, 2026
Merged

293 acoustics save lamb vector from flow simulations and read it from acoustics solver #295

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
3c2d8e1
Lamb vector is saved into a file. Either a solution file Lamb.hsol or…
albertjimenezramos Jan 29, 2026
cf23787
Lamb vector is read from hsol and stats.hsol file for the acoustics s…
albertjimenezramos Jan 30, 2026
ee68b14
Added write and read in pos for parallel i/o
albertjimenezramos Feb 5, 2026
09d1314
Merge branch 'saveAndReadLambVector' into 293-acoustics-save-lamb-vec…
albertjimenezramos Feb 6, 2026
246d386
Fixed wrong intent in.
albertjimenezramos Feb 9, 2026
8ce4c90
Moved reading to time integration function. [skip ci]
albertjimenezramos Feb 9, 2026
46b480d
Change acosutics by acoustic flag [skip ci]
albertjimenezramos Feb 9, 2026
b8831e4
Dummy commit to run workflow for pr
albertjimenezramos Feb 11, 2026
c001189
Remove saveLambVector keyword from mainKeywords
albertjimenezramos Feb 11, 2026
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
5 changes: 3 additions & 2 deletions Solver/src/NavierStokesSolver/main.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,7 @@ PROGRAM HORSES3DMainNS
TYPE( FTValueDictionary) :: controlVariables
TYPE( DGSem ) :: sem
TYPE( TimeIntegrator_t ) :: timeIntegrator
LOGICAL :: success, saveGradients, saveSensor, saveLES, saveSource
LOGICAL :: success, saveGradients, saveSensor, saveLES, saveSource, saveLambVector
logical :: generateMonitor = .TRUE.
logical :: optimizePartitionLevel=.FALSE.
logical :: isMLRK=.FALSE.
Expand Down Expand Up @@ -204,7 +204,8 @@ PROGRAM HORSES3DMainNS
saveSensor = controlVariables % logicalValueForKey(saveSensorToSolutionKey)
saveLES = controlVariables % logicalValueForKey(saveLESToSolutionKey)
saveSource = controlVariables % logicalValueForKey(saveSourceToSolutionKey)
CALL sem % mesh % SaveSolution(sem % numberOfTimeSteps, timeIntegrator % time, solutionFileName, saveGradients, saveSensor, saveLES, saveSource)
saveLambVector = controlVariables % logicalValueForKey(saveLambVectorToSolutionKey)
CALL sem % mesh % SaveSolution(sem % numberOfTimeSteps, timeIntegrator % time, solutionFileName, saveGradients, saveSensor, saveLES, saveSource, saveLambVector)
if ( sem % particles % active ) then
call sem % particles % ExportToVTK ( sem % numberOfTimeSteps, sem % monitors % solution_file )
end if
Expand Down
1 change: 1 addition & 0 deletions Solver/src/libs/foundation/Setup.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@ Module mainKeywordsModule
CHARACTER(LEN=KEYWORD_LENGTH), PARAMETER :: restartFileNameKey = "restart file name"
CHARACTER(LEN=KEYWORD_LENGTH), PARAMETER :: solutionFileNameKey = "solution file name"
CHARACTER(LEN=KEYWORD_LENGTH), PARAMETER :: saveGradientsToSolutionKey = "save gradients with solution"
CHARACTER(LEN=KEYWORD_LENGTH), PARAMETER :: saveLambVectorToSolutionKey = "save lamb vector"
CHARACTER(LEN=KEYWORD_LENGTH), PARAMETER :: saveSensorToSolutionKey = "save sensor with solution"
CHARACTER(LEN=KEYWORD_LENGTH), PARAMETER :: saveLESToSolutionKey = "save les with solution"
CHARACTER(LEN=KEYWORD_LENGTH), PARAMETER :: saveSourceToSolutionKey = "save source with solution"
Expand Down
296 changes: 285 additions & 11 deletions Solver/src/libs/mesh/HexMesh.f90
View file
Open in desktop

Large diffs are not rendered by default.

16 changes: 15 additions & 1 deletion Solver/src/libs/mesh/MappedGeometry.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@ Module MappedGeometryClass
private
public MappedGeometry, MappedGeometryFace
public SetMappingsToCrossProduct
public vDot, vCross
public vDot, vCross, ComputeVorticity
!
! ---------
! Constants
Expand Down Expand Up @@ -860,5 +860,19 @@ end subroutine SetMappingsToCrossProduct
!
!///////////////////////////////////////////////////////////////////////////////
!
pure subroutine ComputeVorticity(U_x, U_y, U_z, vorticity)

real(kind=RP), intent(in) :: U_x(3)
real(kind=RP), intent(in) :: U_y(3)
real(kind=RP), intent(in) :: U_z(3)
real(kind=RP), intent(out) :: vorticity(3)

vorticity(1) = U_y(3) - U_z(2)
vorticity(2) = U_z(1) - U_x(3)
vorticity(3) = U_x(2) - U_y(1)

end subroutine ComputeVorticity
!
!///////////////////////////////////////////////////////////////////////////////
!
END Module MappedGeometryClass
20 changes: 20 additions & 0 deletions Solver/src/libs/mesh/StorageClass.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -94,6 +94,8 @@ module StorageClass
#endif
#ifdef ACOUSTIC
real(kind=RP), allocatable :: Qbase(:,:,:,:) ! Base flow State vector
real(kind=RP), allocatable :: Lambbase(:,:,:,:) ! Lamb vector base
real(kind=RP), allocatable :: Lamb_NS(:,:,:,:) ! Lamb vector from NS
#endif
#ifdef CAHNHILLIARD
real(kind=RP), dimension(:,:,:,:), allocatable :: c ! CHE concentration
Expand Down Expand Up @@ -195,6 +197,8 @@ module StorageClass
real(kind=RP), dimension(:,:), allocatable :: wallNodeDistance ! for BC walls, distance to the first fluid node
#ifdef ACOUSTIC
real(kind=RP), dimension(:,:,:), allocatable :: Qbase ! Base flow State vector
real(kind=RP), dimension(:,:,:), allocatable :: Lambbase ! Lamb vector base
real(kind=RP), dimension(:,:,:), allocatable :: Lamb_NS ! Lamb vector
#endif
#ifdef MULTIPHASE
real(kind=RP), dimension(:,:), allocatable :: invMa2
Expand Down Expand Up @@ -813,6 +817,8 @@ elemental subroutine ElementStorage_Construct(self, Nx, Ny, Nz, computeGradients
#endif
#if defined (ACOUSTIC)
ALLOCATE( self % Qbase (NCONSB,0:Nx,0:Ny,0:Nz) )
ALLOCATE( self % Lambbase (NDIM,0:Nx,0:Ny,0:Nz) )
ALLOCATE( self % Lamb_NS (NDIM,0:Nx,0:Ny,0:Nz) )
#endif
if (computeGradients) then
ALLOCATE( self % U_xNS (NGRAD,0:Nx,0:Ny,0:Nz) )
Expand Down Expand Up @@ -893,6 +899,8 @@ elemental subroutine ElementStorage_Construct(self, Nx, Ny, Nz, computeGradients
#endif
#if defined (ACOUSTIC)
self % Qbase = 0.0_RP
self % Lambbase = 0.0_RP
self % Lamb_NS = 0.0_RP
#endif
if (computeGradients) then
self % U_xNS = 0.0_RP
Expand Down Expand Up @@ -1002,6 +1010,8 @@ elemental subroutine ElementStorage_Assign(to, from)
#endif
#if defined (ACOUSTIC)
to % Qbase = from % Qbase
to % Lambbase = from % Lambbase
to % Lamb_NS = from % Lamb_NS
#endif

#ifndef ACOUSTIC
Expand Down Expand Up @@ -1098,6 +1108,8 @@ elemental subroutine ElementStorage_Destruct(self)
#endif
#if defined (ACOUSTIC)
safedeallocate(self % Qbase)
safedeallocate(self % Lambbase)
safedeallocate(self % Lamb_NS)
#endif

if (self % computeGradients) then
Expand Down Expand Up @@ -1407,6 +1419,8 @@ pure subroutine FaceStorage_Construct(self, NDIM, Nf, Nel, computeGradients, ana
end if
#if defined (ACOUSTIC)
ALLOCATE( self % Qbase (NCONSB,0:Nf(1),0:Nf(2)) )
ALLOCATE( self % Lambbase (NDIM,0:Nf(1),0:Nf(2)) )
ALLOCATE( self % Lamb_NS (NDIM,0:Nf(1),0:Nf(2)) )
#endif
! Biggest Interface flux memory size is u\vec{n}
! ----------------------------------------------
Expand Down Expand Up @@ -1470,6 +1484,8 @@ pure subroutine FaceStorage_Construct(self, NDIM, Nf, Nel, computeGradients, ana
end if
#if defined (ACOUSTIC)
self % Qbase = 0.0_RP
self % Lambbase = 0.0_RP
self % Lamb_NS = 0.0_RP
#endif

self % rho = 0.0_RP
Expand Down Expand Up @@ -1568,6 +1584,8 @@ elemental subroutine FaceStorage_Destruct(self)
safedeallocate(self % rho )
#if defined (ACOUSTIC)
safedeallocate(self % Qbase )
safedeallocate(self % Lambbase )
safedeallocate(self % Lamb_NS )
#endif

self % anJacobian = .FALSE.
Expand Down Expand Up @@ -1732,6 +1750,8 @@ elemental subroutine FaceStorage_Assign(to, from)
to % QNS = from % QNS
#if defined (ACOUSTIC)
to % Qbase = from % Qbase
to % Lambbase = from % Lambbase
to % Lamb_NS = from % Lamb_NS
#endif
if (to % computeGradients) then
to % U_xNS = from % U_xNS
Expand Down
5 changes: 3 additions & 2 deletions Solver/src/libs/monitors/Monitors.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -84,7 +84,7 @@ subroutine Monitors_Construct( Monitors, mesh, controlVariables )
character(len=STR_LEN_MONITORS) :: line
character(len=STR_LEN_MONITORS) :: solution_file
logical, save :: FirstCall = .TRUE.
logical :: saveGradients
logical :: saveGradients, saveLambVector
logical :: monitorMem =.FALSE.
!
! Setup the buffer
Expand Down Expand Up @@ -147,7 +147,8 @@ subroutine Monitors_Construct( Monitors, mesh, controlVariables )

#if defined(NAVIERSTOKES) || defined(INCNS)
saveGradients = controlVariables % logicalValueForKey(saveGradientsToSolutionKey)
call Monitors % stats % Construct(mesh, saveGradients)
saveLambVector = controlVariables % logicalValueForKey(saveLambVectorToSolutionKey)
call Monitors % stats % Construct(mesh, saveGradients, saveLambVector)


allocate ( Monitors % surfaceMonitors ( Monitors % no_of_surfaceMonitors ) )
Expand Down
66 changes: 53 additions & 13 deletions Solver/src/libs/monitors/StatisticsMonitor.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -51,6 +51,7 @@ module StatisticsMonitor
real(kind=RP) :: starting_time
integer :: no_of_samples
logical :: saveGradients
logical :: saveLambVector
contains
procedure :: Construct => StatisticsMonitor_Construct
procedure :: Update => StatisticsMonitor_Update
Expand All @@ -67,23 +68,23 @@ module StatisticsMonitor
!
!//////////////////////////////////////////////////////////////////////////////////
!
subroutine StatisticsMonitor_Construct(self, mesh, saveGradients)
subroutine StatisticsMonitor_Construct(self, mesh, saveGradients, saveLambVector)
use ParamfileRegions
use PhysicsStorage, only: NCONS, NGRAD
use HexMeshClass, only: no_of_stats_variables
implicit none
class(StatisticsMonitor_t) :: self
class(HexMesh) :: mesh
logical, intent(in) :: saveGradients
logical, intent(in) :: saveGradients, saveLambVector
integer, allocatable :: Nsample, i0, Ndump, Nreset
real(kind=RP), allocatable :: t0
integer :: eID
character(len=LINE_LENGTH) :: paramFile

NO_OF_VARIABLES = NO_OF_VARIABLES_Sij + NCONS
self % saveGradients = saveGradients
self % saveLambVector = saveLambVector
if (saveGradients) NO_OF_VARIABLES = NO_OF_VARIABLES + NGRAD * NDIM
no_of_stats_variables = NO_OF_VARIABLES
if (saveLambVector) NO_OF_VARIABLES = NO_OF_VARIABLES + NDIM
!
! Search for the parameters in the case file
! ------------------------------------------
Expand Down Expand Up @@ -168,7 +169,11 @@ subroutine StatisticsMonitor_WriteFile(self, mesh, iter, t, solution_file)
character(len=LINE_LENGTH) :: fileName

write(fileName,'(A,A)') trim(solution_file),'.stats.hsol'
if ( self % state .ne. OFF) call mesh % SaveStatistics(iter, t, trim(fileName), self % saveGradients)
if ( self % state .ne. OFF) then
call mesh % SaveStatistics(iter, t, trim(fileName), self % saveGradients)
write(fileName,'(A,A)') trim(solution_file),'.Lamb.stats.hsol'
if (self % saveLambVector) call mesh % SaveLambVectorStatistics(iter, t, trim(fileName), self % saveGradients)
end if

end subroutine StatisticsMonitor_WriteFile

Expand Down Expand Up @@ -225,6 +230,7 @@ end subroutine StatisticsMonitor_Update

subroutine StatisticsMonitor_UpdateValues(self, mesh)
use PhysicsStorage
use MappedGeometryClass, only: vCross, ComputeVorticity
implicit none
class(StatisticsMonitor_t) :: self
class(HexMesh) :: mesh
Expand All @@ -237,16 +243,23 @@ subroutine StatisticsMonitor_UpdateValues(self, mesh)
integer :: i, j, k
real(RP) :: ratio, inv_nsamples_plus_1
real(RP) :: rfactor1, rfactor2
integer, dimension(5) :: limits
integer, dimension(6) :: limits
real(kind=RP), dimension(NDIM) :: vorticity, velocity, LambVector
real(kind=RP), pointer, contiguous :: U_x(:), U_y(:), U_z(:)


#ifdef NAVIERSTOKES
! if gradients are not saved, limits(2) is equal to limits(5), the latter wont be used
limits(1) = NO_OF_VARIABLES_Sij + IRHO
limits(2) = NO_OF_VARIABLES_Sij + NCONS
limits(3) = NO_OF_VARIABLES_Sij + NCONS + NGRAD
limits(4) = NO_OF_VARIABLES_Sij + NCONS + 2*NGRAD
limits(5) = NO_OF_VARIABLES
limits(3:6) = limits(2)
if (self % saveGradients) then
limits(3) = NO_OF_VARIABLES_Sij + NCONS + NGRAD ! U_x
limits(4) = NO_OF_VARIABLES_Sij + NCONS + 2*NGRAD ! U_y
limits(5) = NO_OF_VARIABLES_Sij + NCONS + NDIM * NGRAD ! U_z
limits(6) = limits(5)
end if
if (self % saveLambVector) limits(6) = limits(5) + NDIM

inv_nsamples_plus_1 = 1.0_RP / (self % no_of_samples + 1)
ratio = self % no_of_samples * inv_nsamples_plus_1
Expand All @@ -272,7 +285,18 @@ subroutine StatisticsMonitor_UpdateValues(self, mesh)
data(limits(2)+1:limits(3),i,j,k) = data(limits(2)+1:limits(3),i,j,k) * ratio + e % storage % U_x(:,i,j,k) * inv_nsamples_plus_1
data(limits(3)+1:limits(4),i,j,k) = data(limits(3)+1:limits(4),i,j,k) * ratio + e % storage % U_y(:,i,j,k) * inv_nsamples_plus_1
data(limits(4)+1:limits(5),i,j,k) = data(limits(4)+1:limits(5),i,j,k) * ratio + e % storage % U_z(:,i,j,k) * inv_nsamples_plus_1
end if
end if
! Save Lamb vector: u x w
if (self % saveLambVector) then
! Attention: qm in APE is -(w x u)'
U_x => e % storage % U_x(:,i,j,k)
U_y => e % storage % U_y(:,i,j,k)
U_z => e % storage % U_z(:,i,j,k)
call ComputeVorticity(U_x, U_y, U_z, vorticity)
velocity = e % storage % Q(IRHOU:IRHOW,i,j,k) / e % storage % Q(IRHO,i,j,k)
call vCross(velocity, vorticity, LambVector)
data(limits(5)+1:limits(6),i,j,k) = data(limits(5)+1:limits(6),i,j,k) * ratio + LambVector * inv_nsamples_plus_1
end if
end do ; end do ; end do

end associate
Expand All @@ -283,9 +307,14 @@ subroutine StatisticsMonitor_UpdateValues(self, mesh)
! if gradients are not saved, limits(2) is equal to limits(5), the latter wont be used
limits(1) = NO_OF_VARIABLES_Sij + INSRHO
limits(2) = NO_OF_VARIABLES_Sij + NCONS
limits(3) = NO_OF_VARIABLES_Sij + NCONS + NGRAD
limits(4) = NO_OF_VARIABLES_Sij + NCONS + 2*NGRAD
limits(5) = NO_OF_VARIABLES
limits(3:6) = limits(2)
if (self % saveGradients) then
limits(3) = NO_OF_VARIABLES_Sij + NCONS + NGRAD ! U_x
limits(4) = NO_OF_VARIABLES_Sij + NCONS + 2*NGRAD ! U_y
limits(5) = NO_OF_VARIABLES_Sij + NCONS + NDIM * NGRAD ! U_z
limits(6) = limits(5)
end if
if (self % saveLambVector) limits(6) = limits(5) + NDIM

inv_nsamples_plus_1 = 1.0_RP / (self % no_of_samples + 1)
ratio = self % no_of_samples * inv_nsamples_plus_1
Expand All @@ -312,6 +341,17 @@ subroutine StatisticsMonitor_UpdateValues(self, mesh)
data(limits(3)+1:limits(4),i,j,k) = data(limits(3)+1:limits(4),i,j,k) * ratio + e % storage % U_y(:,i,j,k) * inv_nsamples_plus_1
data(limits(4)+1:limits(5),i,j,k) = data(limits(4)+1:limits(5),i,j,k) * ratio + e % storage % U_z(:,i,j,k) * inv_nsamples_plus_1
end if
! Save Lamb vector: u x w
if (self % saveLambVector) then
! Attention: qm in APE is -(w x u)'
U_x => e % storage % U_x(:,i,j,k)
U_y => e % storage % U_y(:,i,j,k)
U_z => e % storage % U_z(:,i,j,k)
call ComputeVorticity(U_x, U_y, U_z, vorticity)
velocity = e % storage % Q(INSRHOU:INSRHOW,i,j,k) / e % storage % Q(INSRHO,i,j,k)
call vCross(velocity, vorticity, LambVector)
data(limits(5)+1:limits(6),i,j,k) = data(limits(5)+1:limits(6),i,j,k) * ratio + LambVector * inv_nsamples_plus_1
end if
end do ; end do ; end do

end associate
Expand Down
9 changes: 9 additions & 0 deletions Solver/src/libs/physics/navierstokes/PhysicsStorage_NS.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -145,6 +145,7 @@ MODULE PhysicsStorage_NS
SUBROUTINE ConstructPhysicsStorage_NS( controlVariables, Lref, timeref, success )
USE FTValueDictionaryClass
USE Physics_NSKeywordsModule
use mainKeywordsModule, only: saveLambVectorToSolutionKey
use Utilities, only: toLower, almostEqual
!
! ---------
Expand Down Expand Up @@ -273,6 +274,14 @@ SUBROUTINE ConstructPhysicsStorage_NS( controlVariables, Lref, timeref, success
end if
end if

! If Lamb vector is requested, computeGradient should be enabled
if (controlVariables % logicalValueForKey(saveLambVectorToSolutionKey)) then
if (.not. computeGradients) then
print *, trim(COMPUTE_GRADIENTS_KEY), " should be enabled to compute the Lamb vector."
error stop
end if
end if

dimensionless_ % gammaM2 = thermodynamics_ % gamma * POW2( dimensionless_ % Mach )
!
! ********************
Expand Down
Loading
Loading