From e6ebdee825e4f68335cc40788eb166823c73526e Mon Sep 17 00:00:00 2001 From: Albert Jimenez Ramos Date: Tue, 27 Jan 2026 14:58:29 +0100 Subject: [PATCH 1/4] Qbase now is initialized with either a uniform field set in the control file, or read from a stats.hsol file. --- Solver/src/libs/discretization/DGSEMClass.f90 | 5 +- Solver/src/libs/mesh/HexMesh.f90 | 191 +++++++++++++++++- .../acoustic/VariableConversion_caa.f90 | 2 +- 3 files changed, 191 insertions(+), 7 deletions(-) diff --git a/Solver/src/libs/discretization/DGSEMClass.f90 b/Solver/src/libs/discretization/DGSEMClass.f90 index 19edbb8aaa..5381c1168f 100644 --- a/Solver/src/libs/discretization/DGSEMClass.f90 +++ b/Solver/src/libs/discretization/DGSEMClass.f90 @@ -144,7 +144,6 @@ SUBROUTINE ConstructDGSem( self, meshFileName_, controlVariables, & logical :: useWeightsPartition ! Partitioning mesh using DOF of elements as weights logical :: genMonitor logical :: isReconstruct=.false. - real(kind=RP) :: QbaseUniform(1:NCONS) character(len=*), parameter :: TWOD_OFFSET_DIR_KEY = "2d mesh offset direction" procedure(UserDefinedInitialCondition_f) :: UserDefinedInitialCondition #if (!defined(NAVIERSTOKES)) @@ -399,9 +398,7 @@ SUBROUTINE ConstructDGSem( self, meshFileName_, controlVariables, & ! -------------------- ! #if defined(ACOUSTIC) - ! start by default with no flow conditions - QbaseUniform = [1.0_RP,0.0_RP,0.0_RP,0.0_RP,1.0_RP/(dimensionless % gammaM2)] - call self % mesh % SetUniformBaseFlow(QbaseUniform) + call self % mesh % InitializeBaseFlow(controlVariables) call self % mesh % ProlongBaseSolutionToFaces(NCONS) #ifdef _HAS_MPI_ !$omp single diff --git a/Solver/src/libs/mesh/HexMesh.f90 b/Solver/src/libs/mesh/HexMesh.f90 index c4a2f9d658..2b82bed858 100644 --- a/Solver/src/libs/mesh/HexMesh.f90 +++ b/Solver/src/libs/mesh/HexMesh.f90 @@ -138,8 +138,9 @@ MODULE HexMeshClass procedure :: LoadSolutionForRestart => HexMesh_LoadSolutionForRestart procedure :: WriteCoordFile #if defined(ACOUSTIC) + procedure :: InitializeBaseFlow => HexMesh_InitializeBaseFlow procedure :: SetUniformBaseFlow => HexMesh_SetUniformBaseFlow - ! procedure :: LoadBaseFlowSolution => HexMesh_LoadBaseFlowSolution + procedure :: LoadBaseFlowSolution => HexMesh_LoadBaseFlowSolution procedure :: ProlongBaseSolutionToFaces => HexMesh_ProlongBaseSolutionToFaces procedure :: UpdateMPIFacesBaseSolution => HexMesh_UpdateMPIFacesBaseSolution procedure :: GatherMPIFacesBaseSolution => HexMesh_GatherMPIFacesBaseSolution @@ -3568,7 +3569,7 @@ subroutine HexMesh_SaveStatistics(self, iter, time, name, saveGradients) refs(V_REF) = refValues % V refs(T_REF) = refValues % T refs(MACH_REF) = dimensionless % Mach - refs(RE_REF) = dimensionless % Re + ! refs(RE_REF) = dimensionless % Re ! Create new file ! --------------- @@ -4038,6 +4039,92 @@ END SUBROUTINE HexMesh_LoadSolution !/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! #if defined(ACOUSTIC) + + subroutine HexMesh_InitializeBaseFlow(self, controlVariables) + use FTObjectClass + use mainKeywordsModule + use FileReadingUtilities, only: getRealArrayFromString + Implicit None + CLASS(HexMesh) :: self + class(FTValueDictionary) :: controlVariables + + + class(FTObject), pointer :: obj, obj2 + character(len=LINE_LENGTH) :: qBaseMode + character(len=LINE_LENGTH) :: fileName + real(kind=RP) :: QbaseUniform(1:NCONS) + + CHARACTER(LEN=KEYWORD_LENGTH) :: qBaseKey = "qBase" + CHARACTER(LEN=KEYWORD_LENGTH) :: qBaseFileNameKey = "qBase file name" + CHARACTER(LEN=KEYWORD_LENGTH) :: qBaseVectorKey = "qBase vector" + character(len=LINE_LENGTH) :: qBaseByFile = 'file' + character(len=LINE_LENGTH) :: qBaseByUniformField = 'uniform' + + + ! Check which type of qBase we have: file or uniform + call toLower(qBaseKey) + obj => controlVariables % objectForKey(trim(qBaseKey)) + if ( associated(obj) ) then + + qBaseMode = controlVariables % stringValueForKey(trim(qBaseKey), requestedLength = LINE_LENGTH) + call ToLower(qBaseMode) + + if ( trim(qBaseMode) .eq. trim(qBaseByFile) ) then + ! + ! Read Qbase from file + ! + ! Check that the user has specified the file to read from + call toLower(qBaseFileNameKey) + obj2 => controlVariables % objectForKey(trim(qBaseFileNameKey)) + if ( .not. associated(obj2) ) then + print *, trim(qBaseKey), " = ", trim(qBaseMode), ", but no file specified. Use:" + print *, trim(qBaseFileNameKey), " = filename" + errorMessage(STD_OUT) + error stop + end if + ! Load the base flow from the specified stats file + fileName = controlVariables % stringValueForKey(qBaseFileNameKey,requestedLength = LINE_LENGTH) + call self % LoadBaseFlowSolution(fileName) + elseif ( trim(qBaseMode) .eq. trim(qbaseByUniformField) ) then + ! + ! Read Qbase uniform field from control file + ! + ! Check that the user has specified the field + call toLower(qBaseVectorKey) + obj2 => controlVariables % objectForKey(trim(qBaseVectorKey)) + if ( .not. associated(obj2) ) then + print *, trim(qBaseKey), " = ", trim(qBaseMode), ", but no vector specified. Use:" + print *, trim(qBaseVectorKey), " = [1.0_RP,0.0_RP,0.0_RP,0.0_RP,1.0_RP]" + errorMessage(STD_OUT) + error stop + end if + ! Read the field + QbaseUniform = 0.0_RP + QbaseUniform = GetRealArrayFromString( controlVariables % StringValueForKey(qBaseVectorKey,requestedLength = LINE_LENGTH)) + ! Set uniform field + call self % SetUniformBaseFlow(QbaseUniform) + else + print*, 'Unknown qBase mode "',trim(qBaseMode),'".' + print*, "Implemented modes are:" + print*, " * ", trim(qBaseByFile) + print*, " * ", trim(qbaseByUniformField) + errorMessage(STD_OUT) + error stop + end if + + else + ! Keyword not present: + print*, 'Argument qBase mode is mandatory' + print*, "Implemented modes are:" + print*, " * ", trim(qBaseByFile) + print*, " * ", trim(qbaseByUniformField) + errorMessage(STD_OUT) + error stop + + end if + + end subroutine HexMesh_InitializeBaseFlow + Subroutine HexMesh_SetUniformBaseFlow(self,Q_in) Implicit None CLASS(HexMesh) :: self @@ -4056,6 +4143,106 @@ Subroutine HexMesh_SetUniformBaseFlow(self,Q_in) ! End Subroutine HexMesh_SetUniformBaseFlow ! +!//////////////////////////////////////////////////////////////////////// +! + Subroutine HexMesh_LoadBaseFlowSolution(self, fileName) + use VariableConversion_CAA, only: PressureBaseFlow + Implicit None + CLASS(HexMesh) :: self + character(len=*) :: fileName +! --------------- +! Local variables +! --------------- + INTEGER :: fID, eID, fileType, no_of_elements, flag, nodetype + integer, allocatable :: arrayDimensions(:) + integer :: iter + real(kind=RP) :: time + real(kind=RP), dimension(NO_OF_SAVED_REFS) :: refs + integer :: i, j, k + integer :: padding + integer(kind=AddrInt) :: pos + integer :: no_stat_s + real(kind=RP), allocatable :: Q(:,:,:,:) + character(len=SOLFILE_STR_LEN) :: rstName + logical :: gradients + + gradients = .FALSE. +! +! Get the file title +! ------------------ + rstName = getSolutionFileName(trim(fileName)) +! +! Get the file type +! ----------------- + fileType = getSolutionFileType(trim(fileName)) + + if (fileType .ne. STATS_FILE) then + print*, "The selected file is not a statistics file" + errorMessage(STD_OUT) + error stop + end if +! +! Get the node type +! ----------------- + nodeType = getSolutionFileNodeType(trim(fileName)) + + if ( nodeType .ne. self % nodeType ) then + print*, "WARNING: Solution file uses a different discretization nodes than the mesh." + print*, "Add restart polorder = (Pol order in your restart file) in the control file if you want interpolation routines to be used." + print*, "If restart polorder is not specified the values in the original set of nodes are loaded into the new nodes without interpolation." + errorMessage(STD_OUT) + end if +! +! Read the number of elements +! --------------------------- + no_of_elements = getSolutionFileNoOfElements(trim(fileName)) + if ( no_of_elements .ne. self % no_of_allElements ) then + write(STD_OUT,'(A,A)') "The number of elements stored in the restart file ", & + "do not match that of the mesh file" + errorMessage(STD_OUT) + error stop + end if + allocate(arrayDimensions(4)) +! +! Get time and iteration +! ---------------------- + call getSolutionFileTimeAndIteration(trim(fileName),iter,time) +! +! Read reference values +! --------------------- + refs = getSolutionFileReferenceValues(trim(fileName)) +! +! Read elements data +! ------------------ + fID = putSolutionFileInReadDataMode(trim(fileName)) + do eID = 1, size(self % elements) + associate( e => self % elements(eID) ) + call getSolutionFileArrayDimensions(fID,arrayDimensions) + + no_stat_s = 9 + ! Read and initialize velocity + allocate(Q(1:no_stat_s, 0:e % Nxyz(1), 0:e % Nxyz(2), 0:e % Nxyz(3))) + read(fID) Q + ! WARNING: 1,2,3 should match U,V,W in StatisticsMonitor + self % elements(eID) % storage % Qbase(ICAAU,:,:,:) = Q(1,:,:,:) + self % elements(eID) % storage % Qbase(ICAAV,:,:,:) = Q(2,:,:,:) + self % elements(eID) % storage % Qbase(ICAAW,:,:,:) = Q(3,:,:,:) + deallocate(Q) + ! Read and initialize density + allocate(Q(1:NCONS, 0:e % Nxyz(1), 0:e % Nxyz(2), 0:e % Nxyz(3))) + read(fID) Q + self % elements(eID) % storage % Qbase(ICAARHO,:,:,:) = Q(IRHO,:,:,:) + ! Read and initialize pressure + do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1) + self % elements(eID) % storage % Qbase(ICAAP,i,j,k) = PressureBaseFlow(Q(:,i,j,k)) + end do ; end do ; end do + deallocate(Q) + end associate + end do + + + End Subroutine HexMesh_LoadBaseFlowSolution + !//////////////////////////////////////////////////////////////////////// ! subroutine HexMesh_ProlongBaseSolutionToFaces(self, nEqn) diff --git a/Solver/src/libs/physics/acoustic/VariableConversion_caa.f90 b/Solver/src/libs/physics/acoustic/VariableConversion_caa.f90 index eaf6319c70..25f4f68fce 100644 --- a/Solver/src/libs/physics/acoustic/VariableConversion_caa.f90 +++ b/Solver/src/libs/physics/acoustic/VariableConversion_caa.f90 @@ -6,7 +6,7 @@ module VariableConversion_CAA implicit none private - public Pressure, PressureDot + public Pressure, PressureDot, PressureBaseFlow public NSGradientVariables_STATE ! public getPrimitiveVariables public getVelocityGradients, getTemperatureGradient, getConservativeGradients From 3cb83c0bbba5177f22b9eef699280e5febd56598 Mon Sep 17 00:00:00 2001 From: Albert Jimenez Ramos Date: Wed, 28 Jan 2026 09:30:55 +0100 Subject: [PATCH 2/4] Added explanation of the keywords in the function. --- Solver/src/libs/mesh/HexMesh.f90 | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/Solver/src/libs/mesh/HexMesh.f90 b/Solver/src/libs/mesh/HexMesh.f90 index 2b82bed858..459f252851 100644 --- a/Solver/src/libs/mesh/HexMesh.f90 +++ b/Solver/src/libs/mesh/HexMesh.f90 @@ -4060,6 +4060,13 @@ subroutine HexMesh_InitializeBaseFlow(self, controlVariables) character(len=LINE_LENGTH) :: qBaseByFile = 'file' character(len=LINE_LENGTH) :: qBaseByUniformField = 'uniform' + ! In the control file, the keyword 'qbase' is mandatory: + ! qbase = file/uniform + ! * When qbase is given by a file, the keyword 'qbase file name' is mandatory: + ! qbase file name = path/to/file.stats.hsol + ! * When qbase is given by a uniform field, the keyword 'qbase vector' is mandatory: + ! qbase vector = [1.0_RP,0.0_RP,0.0_RP,0.0_RP,1.0_RP] + ! Check which type of qBase we have: file or uniform call toLower(qBaseKey) From aea82e6fb5b44a7854d554c32633efd6dedd50b5 Mon Sep 17 00:00:00 2001 From: Albert Jimenez Ramos Date: Thu, 5 Feb 2026 10:32:39 +0100 Subject: [PATCH 3/4] Able to read using pos allowing the parallel reading. Added stats_and_gradients file type if we also save the gradients in the stats file. --- Solver/src/libs/io/SolutionFile.f90 | 15 +++++---- Solver/src/libs/mesh/HexMesh.f90 | 47 ++++++++++------------------- 2 files changed, 25 insertions(+), 37 deletions(-) diff --git a/Solver/src/libs/io/SolutionFile.f90 b/Solver/src/libs/io/SolutionFile.f90 index 859c8528db..00df655f5b 100644 --- a/Solver/src/libs/io/SolutionFile.f90 +++ b/Solver/src/libs/io/SolutionFile.f90 @@ -39,7 +39,8 @@ module SolutionFile #endif private - public :: MESH_FILE, SOLUTION_FILE, SOLUTION_AND_GRADIENTS_FILE, STATS_FILE, ZONE_MESH_FILE + public :: MESH_FILE, SOLUTION_FILE, SOLUTION_AND_GRADIENTS_FILE, ZONE_MESH_FILE + public :: STATS_FILE, STATS_AND_GRADIENTS_FILE public :: ZONE_SOLUTION_FILE, ZONE_SOLUTION_AND_DOT_FILE public :: SOLUTION_AND_SENSOR_FILE, SOLUTION_AND_GRADIENTS_AND_SENSOR_FILE public :: BEGINNING_DATA @@ -61,11 +62,12 @@ module SolutionFile integer, parameter :: SOLUTION_FILE = 2 integer, parameter :: SOLUTION_AND_GRADIENTS_FILE = 3 integer, parameter :: STATS_FILE = 4 - integer, parameter :: ZONE_MESH_FILE = 5 - integer, parameter :: ZONE_SOLUTION_FILE = 6 - integer, parameter :: ZONE_SOLUTION_AND_DOT_FILE = 7 - integer, parameter :: SOLUTION_AND_SENSOR_FILE = 8 - integer, parameter :: SOLUTION_AND_GRADIENTS_AND_SENSOR_FILE = 9 + integer, parameter :: STATS_AND_GRADIENTS_FILE = 5 + integer, parameter :: ZONE_MESH_FILE = 6 + integer, parameter :: ZONE_SOLUTION_FILE = 7 + integer, parameter :: ZONE_SOLUTION_AND_DOT_FILE = 8 + integer, parameter :: SOLUTION_AND_SENSOR_FILE = 9 + integer, parameter :: SOLUTION_AND_GRADIENTS_AND_SENSOR_FILE = 10 integer, parameter :: SOLFILE_STR_LEN = 128 integer, parameter :: END_OF_FILE = 99 @@ -133,6 +135,7 @@ subroutine CreateNewSolutionFile(name, type_, nodes, no_of_elements, iter, time, case(SOLUTION_FILE) case(SOLUTION_AND_GRADIENTS_FILE) case(STATS_FILE) + case(STATS_AND_GRADIENTS_FILE) case(ZONE_MESH_FILE) case(ZONE_SOLUTION_FILE) case(ZONE_SOLUTION_AND_DOT_FILE) diff --git a/Solver/src/libs/mesh/HexMesh.f90 b/Solver/src/libs/mesh/HexMesh.f90 index 459f252851..d8fd93fba3 100644 --- a/Solver/src/libs/mesh/HexMesh.f90 +++ b/Solver/src/libs/mesh/HexMesh.f90 @@ -3555,7 +3555,7 @@ subroutine HexMesh_SaveStatistics(self, iter, time, name, saveGradients) ! Local variables ! --------------- ! - integer :: fid, eID + integer :: fid, eID, fileType integer :: no_stat_s integer(kind=AddrInt) :: pos real(kind=RP) :: refs(NO_OF_SAVED_REFS) @@ -3573,7 +3573,9 @@ subroutine HexMesh_SaveStatistics(self, iter, time, name, saveGradients) ! Create new file ! --------------- - call CreateNewSolutionFile(trim(name),STATS_FILE, self % nodeType, self % no_of_allElements, iter, time, refs) + fileType = STATS_FILE + if ( saveGradients .and. computeGradients ) fileType = STATS_AND_GRADIENTS_FILE + call CreateNewSolutionFile(trim(name), fileType, self % nodeType, self % no_of_allElements, iter, time, refs) ! ! Write arrays ! ------------ @@ -3584,7 +3586,6 @@ subroutine HexMesh_SaveStatistics(self, iter, time, name, saveGradients) no_stat_s = 9 call writeArray(fid, e % storage % stats % data(1:no_stat_s,:,:,:), position=pos) allocate(Q(NCONS, 0:e % Nxyz(1), 0:e % Nxyz(2), 0:e % Nxyz(3))) - ! write(fid) e%storage%stats%data(7:,:,:,:) Q(1:NCONS,:,:,:) = e % storage % stats % data(no_stat_s+1:no_stat_s+NCONS,:,:,:) write(fid) Q deallocate(Q) @@ -4160,30 +4161,18 @@ Subroutine HexMesh_LoadBaseFlowSolution(self, fileName) ! --------------- ! Local variables ! --------------- - INTEGER :: fID, eID, fileType, no_of_elements, flag, nodetype - integer, allocatable :: arrayDimensions(:) - integer :: iter - real(kind=RP) :: time - real(kind=RP), dimension(NO_OF_SAVED_REFS) :: refs + INTEGER :: fID, eID, fileType, no_of_elements, nodetype integer :: i, j, k - integer :: padding integer(kind=AddrInt) :: pos - integer :: no_stat_s + integer :: no_stat_s, no_stats_read real(kind=RP), allocatable :: Q(:,:,:,:) - character(len=SOLFILE_STR_LEN) :: rstName - logical :: gradients - gradients = .FALSE. -! -! Get the file title -! ------------------ - rstName = getSolutionFileName(trim(fileName)) ! ! Get the file type ! ----------------- fileType = getSolutionFileType(trim(fileName)) - if (fileType .ne. STATS_FILE) then + if ( (fileType .ne. STATS_FILE) .and. (fileType .ne. STATS_AND_GRADIENTS_FILE) ) then print*, "The selected file is not a statistics file" errorMessage(STD_OUT) error stop @@ -4209,35 +4198,31 @@ Subroutine HexMesh_LoadBaseFlowSolution(self, fileName) errorMessage(STD_OUT) error stop end if - allocate(arrayDimensions(4)) -! -! Get time and iteration -! ---------------------- - call getSolutionFileTimeAndIteration(trim(fileName),iter,time) -! -! Read reference values -! --------------------- - refs = getSolutionFileReferenceValues(trim(fileName)) + ! ! Read elements data ! ------------------ fID = putSolutionFileInReadDataMode(trim(fileName)) + no_stat_s = 9 + no_stats_read = no_stat_s + NCONS + if (fileType .eq. STATS_AND_GRADIENTS_FILE) no_stats_read = no_stats_read + NGRAD*NDIM do eID = 1, size(self % elements) associate( e => self % elements(eID) ) - call getSolutionFileArrayDimensions(fID,arrayDimensions) + pos = POS_INIT_DATA + (e % globID-1)*5_AddrInt*SIZEOF_INT + 1_AddrInt*no_stats_read*e % offsetIO*SIZEOF_RP + pos = pos + 5_AddrInt*SIZEOF_INT ! This is to skip the reading of the dimensions and shape in writeArray - no_stat_s = 9 ! Read and initialize velocity allocate(Q(1:no_stat_s, 0:e % Nxyz(1), 0:e % Nxyz(2), 0:e % Nxyz(3))) - read(fID) Q + read(fID, pos=pos) Q ! WARNING: 1,2,3 should match U,V,W in StatisticsMonitor self % elements(eID) % storage % Qbase(ICAAU,:,:,:) = Q(1,:,:,:) self % elements(eID) % storage % Qbase(ICAAV,:,:,:) = Q(2,:,:,:) self % elements(eID) % storage % Qbase(ICAAW,:,:,:) = Q(3,:,:,:) deallocate(Q) - ! Read and initialize density + ! Read NCONS variables allocate(Q(1:NCONS, 0:e % Nxyz(1), 0:e % Nxyz(2), 0:e % Nxyz(3))) read(fID) Q + ! Initialize density self % elements(eID) % storage % Qbase(ICAARHO,:,:,:) = Q(IRHO,:,:,:) ! Read and initialize pressure do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1) From 2552b6c4a159c9be54eb74c862adf5d82c00232f Mon Sep 17 00:00:00 2001 From: Albert Jimenez Ramos Date: Thu, 5 Feb 2026 16:48:23 +0100 Subject: [PATCH 4/4] Fixed NCONSB missing. --- Solver/src/libs/discretization/DGSEMClass.f90 | 2 +- Solver/src/libs/mesh/HexMesh.f90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Solver/src/libs/discretization/DGSEMClass.f90 b/Solver/src/libs/discretization/DGSEMClass.f90 index 0d6ca4a10c..541b1d28af 100644 --- a/Solver/src/libs/discretization/DGSEMClass.f90 +++ b/Solver/src/libs/discretization/DGSEMClass.f90 @@ -399,7 +399,7 @@ SUBROUTINE ConstructDGSem( self, meshFileName_, controlVariables, & ! #if defined(ACOUSTIC) call self % mesh % InitializeBaseFlow(controlVariables) - call self % mesh % ProlongBaseSolutionToFaces(NCONS) + call self % mesh % ProlongBaseSolutionToFaces(NCONSB) #ifdef _HAS_MPI_ !$omp single call self % mesh % UpdateMPIFacesBaseSolution(NCONSB) diff --git a/Solver/src/libs/mesh/HexMesh.f90 b/Solver/src/libs/mesh/HexMesh.f90 index e764d471ef..aa5d8cf54c 100644 --- a/Solver/src/libs/mesh/HexMesh.f90 +++ b/Solver/src/libs/mesh/HexMesh.f90 @@ -4053,7 +4053,7 @@ subroutine HexMesh_InitializeBaseFlow(self, controlVariables) class(FTObject), pointer :: obj, obj2 character(len=LINE_LENGTH) :: qBaseMode character(len=LINE_LENGTH) :: fileName - real(kind=RP) :: QbaseUniform(1:NCONS) + real(kind=RP) :: QbaseUniform(1:NCONSB) CHARACTER(LEN=KEYWORD_LENGTH) :: qBaseKey = "qBase" CHARACTER(LEN=KEYWORD_LENGTH) :: qBaseFileNameKey = "qBase file name"