diff --git a/Solver/src/AcousticSolver/main.f90 b/Solver/src/AcousticSolver/main.f90 index 3b40600e67..aa62ee5938 100644 --- a/Solver/src/AcousticSolver/main.f90 +++ b/Solver/src/AcousticSolver/main.f90 @@ -114,7 +114,7 @@ PROGRAM HORSES3DMainCAA ! Integrate in time ! ----------------- ! - CALL timeIntegrator % integrate(sem, controlVariables, sem % monitors, ComputeTimeDerivative, ComputeTimeDerivativeIsolated) + CALL timeIntegrator % integrate(sem, controlVariables, sem % monitors, sem % samplings, ComputeTimeDerivative, ComputeTimeDerivativeIsolated) ! ! ------------------------------------------ ! Finish measuring the total simulation time diff --git a/Solver/src/CahnHilliardSolver/main.f90 b/Solver/src/CahnHilliardSolver/main.f90 index 367849c566..1ef02aee94 100644 --- a/Solver/src/CahnHilliardSolver/main.f90 +++ b/Solver/src/CahnHilliardSolver/main.f90 @@ -122,7 +122,7 @@ PROGRAM HORSES3DMainCH ! Integrate in time ! ----------------- ! - CALL timeIntegrator % integrate(sem, controlVariables, sem % monitors, ComputeTimeDerivative, ComputeTimeDerivativeIsolated) + CALL timeIntegrator % integrate(sem, controlVariables, sem % monitors, sem % samplings, ComputeTimeDerivative, ComputeTimeDerivativeIsolated) ! ! ------------------------------------------ ! Finish measuring the total simulation time diff --git a/Solver/src/IncompressibleNSSolver/main.f90 b/Solver/src/IncompressibleNSSolver/main.f90 index 5570df2ae9..9090286e88 100644 --- a/Solver/src/IncompressibleNSSolver/main.f90 +++ b/Solver/src/IncompressibleNSSolver/main.f90 @@ -115,7 +115,7 @@ PROGRAM HORSES3DMainiNS ! Integrate in time ! ----------------- ! - CALL timeIntegrator % integrate(sem, controlVariables, sem % monitors, ComputeTimeDerivative, ComputeTimeDerivativeIsolated) + CALL timeIntegrator % integrate(sem, controlVariables, sem % monitors, sem % samplings, ComputeTimeDerivative, ComputeTimeDerivativeIsolated) ! ! ---------------------------------- ! Export particles to VTK (temporal) diff --git a/Solver/src/MultiphaseSolver/SpatialDiscretization.f90 b/Solver/src/MultiphaseSolver/SpatialDiscretization.f90 index 2de5ef3bd5..00d76518b4 100644 --- a/Solver/src/MultiphaseSolver/SpatialDiscretization.f90 +++ b/Solver/src/MultiphaseSolver/SpatialDiscretization.f90 @@ -192,15 +192,17 @@ subroutine Initialize_SpaceAndTimeMethods(controlVariables, mesh) end select use_non_constant_speed_of_sound = controlVariables % ContainsKey(FLUID1_COMPRESSIBILITY_KEY) + + call CHDiscretization % Construct(controlVariables, ELLIPTIC_CH) + call CHDiscretization % Describe + + if ( .not. MPI_Process % isRoot ) return + if(use_non_constant_speed_of_sound) then write(STD_OUT,'(A)') " Implementing artificial compressibility with a non-constant speed of sound in each phase" else - write(STD_OUT,'(A)') " Implementing artificial compressibility with a constant speed of sound in each phase" + write(STD_OUT,'(A)') " Implementing artificial compressibility with a constant ACM factor in each phase" endif - - call CHDiscretization % Construct(controlVariables, ELLIPTIC_CH) - call CHDiscretization % Describe - end if diff --git a/Solver/src/MultiphaseSolver/main.f90 b/Solver/src/MultiphaseSolver/main.f90 index 4b78d55062..2081f0f3aa 100644 --- a/Solver/src/MultiphaseSolver/main.f90 +++ b/Solver/src/MultiphaseSolver/main.f90 @@ -115,7 +115,7 @@ PROGRAM HORSES3DMainMU ! Integrate in time ! ----------------- ! - CALL timeIntegrator % integrate(sem, controlVariables, sem % monitors, ComputeTimeDerivative, ComputeTimeDerivativeIsolated) + CALL timeIntegrator % integrate(sem, controlVariables, sem % monitors, sem % samplings, ComputeTimeDerivative, ComputeTimeDerivativeIsolated) ! ! ---------------------------------- ! Export particles to VTK (temporal) diff --git a/Solver/src/NavierStokesSolver/SpatialDiscretization.f90 b/Solver/src/NavierStokesSolver/SpatialDiscretization.f90 index ac44a61c29..2a46836ae7 100644 --- a/Solver/src/NavierStokesSolver/SpatialDiscretization.f90 +++ b/Solver/src/NavierStokesSolver/SpatialDiscretization.f90 @@ -243,7 +243,7 @@ SUBROUTINE ComputeTimeDerivative( mesh, particles, time, mode, HO_Elements, elem ! Local variables ! --------------- ! - INTEGER :: k, locLevel + INTEGER :: k logical :: HOElements if (present(HO_Elements)) then @@ -308,7 +308,7 @@ SUBROUTINE ComputeTimeDerivative( mesh, particles, time, mode, HO_Elements, elem if (present(Level)) then call TimeDerivative_ComputeQDotMLRK(mesh = mesh , & particles = particles, & - t = time, Level=locLevel) + t = time, Level=Level) else call TimeDerivative_ComputeQDot(mesh = mesh , & particles = particles, & @@ -713,10 +713,6 @@ subroutine TimeDerivative_ComputeQDotMLRK( mesh , particles, t, Level) locLevel = Level else locLevel = 1 - - if (.not.allocated(mesh % MLRK % MLIter)) then - call mesh % MLRK % construct(mesh, 1) ! default 1 level - end if end if associate ( MLRK => mesh % MLRK) diff --git a/Solver/src/NavierStokesSolver/main.f90 b/Solver/src/NavierStokesSolver/main.f90 index 7c8c5fc260..b0743c601e 100644 --- a/Solver/src/NavierStokesSolver/main.f90 +++ b/Solver/src/NavierStokesSolver/main.f90 @@ -131,7 +131,7 @@ PROGRAM HORSES3DMainNS ! Integrate in time ! ----------------- ! - CALL timeIntegrator % integrate(sem, controlVariables, sem % monitors, ComputeTimeDerivative, ComputeTimeDerivativeIsolated) + CALL timeIntegrator % integrate(sem, controlVariables, sem % monitors, sem % samplings, ComputeTimeDerivative, ComputeTimeDerivativeIsolated) ! ! ------------------------------------------ ! Finish measuring the total simulation time diff --git a/Solver/src/NavierStokesSolverRANS/main.f90 b/Solver/src/NavierStokesSolverRANS/main.f90 index ad317c6708..7fa8bbba79 100644 --- a/Solver/src/NavierStokesSolverRANS/main.f90 +++ b/Solver/src/NavierStokesSolverRANS/main.f90 @@ -131,7 +131,7 @@ PROGRAM HORSES3DMainNSSA ! Integrate in time ! ----------------- ! - CALL timeIntegrator % integrate(sem, controlVariables, sem % monitors, ComputeTimeDerivative, ComputeTimeDerivativeIsolated) + CALL timeIntegrator % integrate(sem, controlVariables, sem % monitors, sem % samplings, ComputeTimeDerivative, ComputeTimeDerivativeIsolated) ! ! ------------------------------------------ ! Finish measuring the total simulation time diff --git a/Solver/src/libs/discretization/DGSEMClass.f90 b/Solver/src/libs/discretization/DGSEMClass.f90 index 136738ebdb..e3cba32522 100644 --- a/Solver/src/libs/discretization/DGSEMClass.f90 +++ b/Solver/src/libs/discretization/DGSEMClass.f90 @@ -25,6 +25,7 @@ Module DGSEMClass use SpallartAlmarasTurbulence , only: Spalart_Almaras_t #endif use MonitorsClass + use Samplings use ParticlesClass use Physics use FluidData @@ -50,6 +51,7 @@ Module DGSEMClass TYPE(HexMesh) :: mesh LOGICAL :: ManufacturedSol = .FALSE. ! Use manufactured solutions? default .FALSE. type(Monitor_t) :: monitors + type(Sampling_t) :: samplings #if defined(NAVIERSTOKES) && (!(SPALARTALMARAS)) type(FWHClass) :: fwh #endif @@ -432,11 +434,12 @@ SUBROUTINE ConstructDGSem( self, meshFileName_, controlVariables, & END IF #endif ! -! ------------------ -! Build the monitors -! ------------------ +! -------------------------------- +! Build the monitors and samplings +! -------------------------------- ! call self % monitors % construct (self % mesh, controlVariables) + call self % samplings % construct (self % mesh, controlVariables) ! ! ------------------ ! Build the FWH general class @@ -483,6 +486,7 @@ SUBROUTINE DestructDGSem( self ) CALL self % mesh % destruct call self % monitors % destruct + call self % samplings % destruct #if defined(NAVIERSTOKES) && (!(SPALARTALMARAS)) IF (flowIsNavierStokes) call self % fwh % destruct diff --git a/Solver/src/libs/monitors/Makefile b/Solver/src/libs/monitors/Makefile index 61bd6f46da..9a3075f332 100644 --- a/Solver/src/libs/monitors/Makefile +++ b/Solver/src/libs/monitors/Makefile @@ -20,6 +20,11 @@ LIBOBJS_NS = ./build_ns/MonitorDefinitions.o \ ./build_ns/VolumeMonitor.o \ ./build_ns/LoadBalancingMonitor.o \ ./build_ns/Monitors.o \ + ./build_ns/SamplingOperator.o \ + ./build_ns/SurfaceSampling.o \ + ./build_ns/PlaneSampling.o \ + ./build_ns/SpatialMeanNode.o \ + ./build_ns/Samplings.o \ ./build_ns/FWHDefinitions.o \ ./build_ns/FWHObseverClass.o \ ./build_ns/FWHGeneralClass.o @@ -33,7 +38,12 @@ LIBOBJS_NSSA = ./build_nssa/MonitorDefinitions.o \ ./build_nssa/SurfaceMonitor.o \ ./build_nssa/VolumeMonitor.o \ ./build_nssa/LoadBalancingMonitor.o \ - ./build_nssa/Monitors.o + ./build_nssa/Monitors.o \ + ./build_nssa/SamplingOperator.o \ + ./build_nssa/SurfaceSampling.o \ + ./build_nssa/PlaneSampling.o \ + ./build_nssa/SpatialMeanNode.o \ + ./build_nssa/Samplings.o LIBOBJS_iNS = ./build_ins/MonitorDefinitions.o \ ./build_ins/Probe.o \ @@ -44,7 +54,12 @@ LIBOBJS_iNS = ./build_ins/MonitorDefinitions.o \ ./build_ins/SurfaceMonitor.o \ ./build_ins/VolumeMonitor.o \ ./build_ins/LoadBalancingMonitor.o \ - ./build_ins/Monitors.o + ./build_ins/Monitors.o \ + ./build_ins/SamplingOperator.o \ + ./build_ins/SurfaceSampling.o \ + ./build_ins/PlaneSampling.o \ + ./build_ins/SpatialMeanNode.o \ + ./build_ins/Samplings.o LIBOBJS_CH = ./build_ch/MonitorDefinitions.o \ ./build_ch/Probe.o\ @@ -55,7 +70,12 @@ LIBOBJS_CH = ./build_ch/MonitorDefinitions.o \ ./build_ch/SurfaceMonitor.o \ ./build_ch/VolumeMonitor.o \ ./build_ch/LoadBalancingMonitor.o \ - ./build_ch/Monitors.o + ./build_ch/Monitors.o \ + ./build_ch/SamplingOperator.o \ + ./build_ch/SurfaceSampling.o \ + ./build_ch/PlaneSampling.o \ + ./build_ch/SpatialMeanNode.o \ + ./build_ch/Samplings.o LIBOBJS_MU = ./build_mu/MonitorDefinitions.o \ ./build_mu/Probe.o\ @@ -66,7 +86,12 @@ LIBOBJS_MU = ./build_mu/MonitorDefinitions.o \ ./build_mu/SurfaceMonitor.o \ ./build_mu/VolumeMonitor.o \ ./build_mu/LoadBalancingMonitor.o \ - ./build_mu/Monitors.o + ./build_mu/Monitors.o \ + ./build_mu/SamplingOperator.o \ + ./build_mu/SurfaceSampling.o \ + ./build_mu/PlaneSampling.o \ + ./build_mu/SpatialMeanNode.o \ + ./build_mu/Samplings.o LIBOBJS_CAA = ./build_caa/MonitorDefinitions.o \ ./build_caa/Probe.o\ @@ -77,7 +102,12 @@ LIBOBJS_CAA = ./build_caa/MonitorDefinitions.o \ ./build_caa/SurfaceMonitor.o \ ./build_caa/VolumeMonitor.o \ ./build_caa/LoadBalancingMonitor.o \ - ./build_caa/Monitors.o + ./build_caa/Monitors.o \ + ./build_caa/SamplingOperator.o \ + ./build_caa/SurfaceSampling.o \ + ./build_caa/PlaneSampling.o \ + ./build_caa/SpatialMeanNode.o \ + ./build_caa/Samplings.o LIB = monitors diff --git a/Solver/src/libs/monitors/PlaneSampling.f90 b/Solver/src/libs/monitors/PlaneSampling.f90 new file mode 100644 index 0000000000..e5ae70dd5e --- /dev/null +++ b/Solver/src/libs/monitors/PlaneSampling.f90 @@ -0,0 +1,995 @@ +#include "Includes.h" + +module PlaneSampling + use SMConstants + use HexMeshClass + use MonitorDefinitions + use PhysicsStorage + use VariableConversion + use MPI_Process_Info + use FluidData + use FileReadingUtilities, only: getRealArrayFromString, GetRealValue, getCharArrayFromString + use NodalStorageClass , only: NodalStorage +#ifdef _HAS_MPI_ + use mpi +#endif + implicit none + + private + public PlaneSampling_t +! +! ******************************* +! PlaneSamplings class definition +! ******************************* +! + + type PlaneSampling_t + logical , allocatable :: active(:) + integer , allocatable :: rank (:) + integer :: ID + integer :: nVariables + integer :: interval + integer :: bufferSize + integer :: bufferLine + integer :: intervalCount + integer :: N(2) + integer :: nNodes + integer , allocatable :: eID (:) + real(kind=RP), allocatable :: lxi(:,:) , leta(:,:), lzeta(:,:) + real(kind=RP), allocatable :: values(:,:,:) + real(kind=RP), allocatable :: x(:,:) + real(kind=RP), allocatable :: xi(:,:) + logical :: disturbanceData =.false. + character(len=STR_LEN_MONITORS), allocatable :: fileName (:) + character(len=STR_LEN_MONITORS) :: planeName + character(len=STR_LEN_MONITORS) :: fileInput + character(len=STR_LEN_MONITORS), allocatable :: variable (:) + contains + procedure :: Initialization => Plane_Initialization + procedure :: Update => Plane_Update + procedure :: UpdateInterp => Plane_UpdateLagrangeInterp + procedure :: WriteToFile => Plane_WriteToFile + procedure :: LookInOtherPartitions => Plane_LookInOtherPartitions + procedure :: destruct => Plane_Destruct + procedure :: copy => Plane_Assign + generic :: assignment(=) => copy + end type PlaneSampling_t + + contains + + subroutine Plane_Initialization(self, mesh, ID, solution_file, FirstCall) + use Headers + use ParamfileRegions + use MPI_Process_Info + use Utilities, only: toLower + implicit none + class(PlaneSampling_t) :: self + class(HexMesh) :: mesh + integer :: ID + character(len=*) :: solution_file + logical, intent(in) :: FirstCall +! +! --------------- +! Local variables +! --------------- +! + integer :: i, j, k, l, fid, nNodes,inputType, ierr + logical :: firstOutDomain + real(kind=RP) :: point_1(NDIM,1),point_2(NDIM,1),point_3(NDIM,1) + real(kind=RP), allocatable :: pStart(:,:), pEnd(:,:) + real(kind=RP) :: x(NDIM), xi(NDIM) + real(kind=RP) :: delta(NDIM,2), invNodes(2) + real(kind=RP), allocatable :: deltaArray1(:,:) + character(len=STR_LEN_MONITORS) :: in_label + character(len=STR_LEN_MONITORS) :: fileName + character(len=STR_LEN_MONITORS) :: paramFile + character(len=STR_LEN_MONITORS) :: interval + character(len=STR_LEN_MONITORS) :: point1_char,point2_char,point3_char + character(len=STR_LEN_MONITORS) :: N12,N23, inputType_char + character(len=STR_LEN_MONITORS) :: variables + character(len=STR_LEN_MONITORS) :: fileFormat + character(len=STR_LEN_MONITORS) :: writeInterval + character(len=STR_LEN_MONITORS) :: distData + + firstOutDomain = .FALSE. + + if (FirstCall) then +! +! Get monitor ID, assign zero to bufferLine and intervalCount +! -------------- + self % ID = ID + self % bufferLine = 0 + self % intervalCount = 0 +! +! Search for the parameters in the case file +! ------------------------------------------ + write(in_label , '(A,I0)') "#define plane sampling " , self % ID + + call get_command_argument(1, paramFile) + call readValueInRegion(trim(paramFile), "name" , self % planeName , in_label, "#end" ) + call readValueInRegion(trim(paramFile), "variables" , variables , in_label, "#end" ) + call readValueInRegion(trim(paramFile), "node input type", inputType_char , in_label, "#end" ) + if (len(trim(inputType_char))>0) then + inputType = GetRealValue(inputType_char) + else + inputType = 0 + end if + if (inputType.ne.0) then + call readValueInRegion(trim(paramFile), "input file" , self % fileInput , in_label, "#end" ) + else + call readValueInRegion(trim(paramFile), "point 1 [x,y,z]", point1_char , in_label, "#end" ) + call readValueInRegion(trim(paramFile), "point 2 [x,y,z]", point2_char , in_label, "#end" ) + call readValueInRegion(trim(paramFile), "point 3 [x,y,z]", point3_char , in_label, "#end" ) + call readValueInRegion(trim(paramFile), "n12" , N12 , in_label, "#end" ) + call readValueInRegion(trim(paramFile), "n23" , N23 , in_label, "#end" ) + end if + call readValueInRegion(trim(paramFile), "sampling interval", interval , in_label, "#end" ) + call readValueInRegion(trim(paramFile), "write interval", writeInterval , in_label, "#end" ) + call readValueInRegion(trim(paramFile), "data type", distData, in_label, "#end" ) + +! +! Get the variables, points, N discretization, and interval +! ---------------------------------------------- + call getCharArrayFromString(variables, STR_LEN_MONITORS, self % variable) + self % nVariables = size(self % variable) + if (inputType.eq.0) then + point_1(:,1) = getRealArrayFromString(point1_char) + point_2(:,1) = getRealArrayFromString(point2_char) + point_3(:,1) = getRealArrayFromString(point3_char) + self % N(1) = GetRealValue(N12) + self % N(2) = GetRealValue(N23) + else +! Open file input for nodes + open(newunit=fid, file= trim(self % fileInput),status="old",action="read", iostat=ierr) + if (ierr /= 0) then + write(*,*) "ERROR-opening input file for plane sampling: ", trim(self % fileInput) + stop + end if + read(fid, *, iostat=ierr) self % N(1), self % N(2) + end if + self % interval = GetRealValue(interval) + + + interval=TRIM(ADJUSTL(interval)) + writeInterval=TRIM(ADJUSTL(writeInterval)) + +! +! Failsafe if the number of node is less than 2 +! --------------------------------------------- + IF ((self % N(1).LT.2)) THEN + self % N(1)=2 + END IF + IF ((self % N(2).LT.2)) THEN + self % N(2)=2 + END IF + + self % nNodes = self % N(1) * self % N(2) + nNodes = self % nNodes +! +! Get the max. number of timestep in the buffer file before being written +! ----------------------------------------------------------------------- + IF (LEN(TRIM(writeInterval)) .EQ. 0) THEN + self % bufferSize = 1; + ELSE + self % bufferSize = GetRealValue(writeInterval) +! +! Failsafe to prevent too many data being written at one time +! ----------------------------------------------------------- + IF (nNodes * self % bufferSize .GT. 500000) THEN + self % bufferSize = 1; + END IF + END IF +! +! Allocate Variables +! ------------------ + ALLOCATE(self % x (NDIM, nNodes), self % xi (NDIM, nNodes), self % active (nNodes) & + , self % eID (nNodes), self % rank (nNodes), self % fileName (self % nVariables) & + , self % values(nNodes+1,self % bufferSize,self % nVariables), deltaArray1(NDIM,self % N(1)) & + , pStart(NDIM, self % N(1))) +! +! Allocate Lagrange interpolants - Assumed identical polynomial for all elements +! ----------------------------------------------------------------------------- + ALLOCATE(self % lxi( 0 : mesh % elements (1) % Nxyz(1), 1:nNodes), & + self % leta( 0 : mesh % elements (1) % Nxyz(2), 1:nNodes), & + self % lzeta( 0 : mesh % elements (1) % Nxyz(3), 1:nNodes)) +! +! Nodes spacing +! ------------- + invNodes(1)=1.0_RP/(self % N(1)-1) + invNodes(2)=1.0_RP/(self % N(2)-1) + if (inputType.eq.0) then + delta(1,1)=(point_2(1,1)-point_1(1,1))*invNodes(1) + delta(2,1)=(point_2(2,1)-point_1(2,1))*invNodes(1) + delta(3,1)=(point_2(3,1)-point_1(3,1))*invNodes(1) + delta(1,2)=(point_3(1,1)-point_2(1,1))*invNodes(2) + delta(2,2)=(point_3(2,1)-point_2(2,1))*invNodes(2) + delta(3,2)=(point_3(3,1)-point_2(3,1))*invNodes(2) +! +! Create spacing Array for Nodes in direction 1 +! --------------------------------------------- + DO l=1,self % N(1) + deltaArray1(1,l)=(delta(1,1))*(l-1) + deltaArray1(2,l)=(delta(2,1))*(l-1) + deltaArray1(3,l)=(delta(3,1))*(l-1) + END DO +! +! Get Nodes location of the plane +! ------------------------------- + DO l=1,self % N(2) + DO k=1, self % N(1) + DO j=1,3 + pStart(j,k) = (point_1(j,1)+delta(j,2)*(l-1)) + self % x(j,(l-1)*self % N(1)+k) = pStart(j,k) + deltaArray1(j,k) + END DO + END DO + END DO + else + ALLOCATE(pEnd(NDIM, self % N(1))) +! +! Read first points location from input file +! ---------------------------------- + DO l=1,self % N(1) + read(fid, *, iostat=ierr) pStart(1,l), pStart(2,l), pStart(3,l) + if (ierr /= 0) exit + END DO +! +! Read last points location from input file +! ---------------------------------- + DO l=1,self % N(1) + read(fid, *, iostat=ierr) pEnd(1,l), pEnd(2,l), pEnd(3,l) + if (ierr /= 0) exit + END DO + close(fid) +! +! Create spacing Array for Nodes in sweep direction +! ------------------------------------------------- + DO l=1,self % N(1) + deltaArray1(1,l)=(pEnd(1,l)-pStart(1,l))*invNodes(2) + deltaArray1(2,l)=(pEnd(2,l)-pStart(2,l))*invNodes(2) + deltaArray1(3,l)=(pEnd(3,l)-pStart(3,l))*invNodes(2) + END DO +! +! Get Nodes location of the plane +! ------------------------------- + DO l=1,self % N(2) + DO j=1,NDIM + DO k=1, self % N(1) + self % x(j,(l-1)*self % N(1)+k) = pStart(j,k) + deltaArray1(j,k)*(l-1) + END DO + END DO + END DO + DEALLOCATE(pEnd) + end if + + + DEALLOCATE (deltaArray1,pStart) + +! +! Find node location in the element ID +! ------------------------------------ + DO i=1,self % nNodes + + x(1)=self % x(1,i) + x(2)=self % x(2,i) + x(3)=self % x(3,i) +! +! Find the requested point in the mesh +! ------------------------------------ + self % active (i) = mesh % FindPointWithCoords(x, self % eID(i), xi) + + self % xi(1,i) = xi(1) + self % xi(2,i) = xi(2) + self % xi(3,i) = xi(3) +! +! Check whether the Plane is located in other partition +! ----------------------------------------------------- + call self % LookInOtherPartitions (i) +! +! Disable the nodes if the point is not found +! ------------------------------------------- + IF ( .not. self % active (i) ) then + IF ( .not. firstOutDomain) then + IF ( MPI_Process % isRoot ) then + firstOutDomain = .TRUE. + END IF + END IF + END IF +! +! Get the Lagrange interpolants +! ----------------------------- + if ( (MPI_Process % rank .ne. self % rank (i)).OR.( .not. self % active (i) ) ) then + self % lxi (:,i)= 0.0_RP + self % leta (:,i)= 0.0_RP + self % lzeta (:,i)= 0.0_RP + ELSE + associate(e => mesh % elements(self % eID(i))) + associate( spAxi => NodalStorage(e % Nxyz(1)), & + spAeta => NodalStorage(e % Nxyz(2)), & + spAzeta => NodalStorage(e % Nxyz(3)) ) + self % lxi (:,i) = spAxi % lj(self % xi(1,i)) + self % leta (:,i) = spAeta % lj(self % xi(2,i)) + self % lzeta (:,i)= spAzeta % lj(self % xi(3,i)) + end associate + end associate + END IF + + END DO + + end if +! +! Check Variables, Create Files, and Write Header Files +! ----------------------------------------------------- + do i=1,self % nVariables +! +! Prepare the file in which the Plane is exported +! ----------------------------------------------- + write( self % fileName (i) , '(A,A,A,A,A,A,I0,A)') trim(solution_file) ,"_", trim(self % planeName) ,"_", trim(self % variable(i)) & + , "_plane_" , self % ID,".sampling" +! +! Check the variable +! ------------------ + call tolower(self % variable (i)) + + select case ( trim(self % variable (i)) ) +#ifdef NAVIERSTOKES + case ("density") + case ("pressure") + case ("ptotal") + case ("velocity") + case ("viscosity") + case ("u") + case ("v") + case ("w") + case ("mach") + case ("k") + case ("omegax") + case ("omegay") + case ("omegaz") + case ("q1") + case ("q2") + case ("q3") + case ("q4") + case ("q5") + case default + if ( MPI_Process % isRoot ) then + print*, 'Plane variable "',trim(self % variable(i)),'" not implemented.' + print*, "Options available are:" + print*, " * density" + print*, " * pressure" + print*, " * velocity" + print*, " * viscosity" + print*, " * u" + print*, " * v" + print*, " * w" + print*, " * Mach" + print*, " * K" + print*, " * q1" + print*, " * q2" + print*, " * q3" + print*, " * q4" + print*, " * q5" + end if +#endif +#ifdef INCNS + case default + print*, "Planes are not implemented for the incompressible NSE" +#endif +#ifdef MULTIPHASE + case ("density") + case ("pressure") + case ("static-pressure") + case ("concentration") + case ("velocity") + case ("u") + case ("v") + case ("w") + case ("omegax") + case ("omegay") + case ("omegaz") + case ("q1") + case ("q2") + case ("q3") + case ("q4") + case ("q5") + case default + if ( MPI_Process % isRoot ) then + print*, 'Plane variable "',trim(self % variable(i)),'" not implemented.' + print*, "Options available are:" + print*, " * density" + print*, " * pressure" + print*, " * static-pressure" + print*, " * velocity" + print*, " * concentration" + print*, " * u" + print*, " * v" + print*, " * w" + print*, " * omegax" + print*, " * omegay" + print*, " * omegaz" + print*, " * q1" + print*, " * q2" + print*, " * q3" + print*, " * q4" + print*, " * q5" + end if +#endif + end select + + if ( .not. MPI_Process % isRoot ) return + + +! +! Create file +! ----------- + if (FirstCall) then + if (MPI_Process % isRoot ) then + open ( newunit = fID , file = trim(self % fileName (i)) , action = "write" , access = "stream" , status = "replace", position='append' ) +! +! Write the file headers +! ---------------------- + write( fID) self % ID + write( fID) self % N(1) + write( fID) self % N(2) + write( fID) self % nNodes + write( fID ) self % interval +#if defined(NAVIERSTOKES) && (!(INCNS)) + write( fID ) refValues % rho + write( fID ) refValues % V + write( fID ) refValues % p + write( fID ) refValues % T + write( fID ) refValues % mu + write( fID ) refValues % AoATheta + write( fID ) refValues % AoAPhi +#elif defined(MULTIPHASE) + write( fID ) refValues % rho + write( fID ) refValues % V + write( fID ) refValues % p + write( fID ) thermodynamics % rho(2) + write( fID ) thermodynamics % mu(1) + write( fID ) thermodynamics % mu(2) + write( fID ) refValues % g0 +#endif + write( fID ) transpose(self % x) + close ( fID ) + end if + end if + + end do +! +! File Format +! ----------------------------------------------- + write( fileFormat , '(A,A,A,A,I0,A)') trim(solution_file) ,"_", trim(self % planeName) ,"_'variable'_plane_", self % ID, ".sampling" +! +! Write Information +! ----------------------------------------------- + write(STD_OUT,'(/)') + call SubSection_Header("Plane Samplings") + write(STD_OUT,'(30X,A,A27,I4)') "->" , "Plane ID: " , self % ID + write(STD_OUT,'(30X,A,A27,A128)') "->" , "Variables: " , variables + if (inputType.eq.0) then + write(STD_OUT,'(30X,A,A27,A36)') "->" , "Input Type: " , "Specified 3 points location" + write(STD_OUT,'(30X,A,A27,A,F7.2,A,F7.2,A,F7.2,A)') "->" , "Point x1 [m]: ","[", & + point_1(1,1), ", ", point_1(2,1), ", ", point_1(3,1), "]" + write(STD_OUT,'(30X,A,A27,A,F7.2,A,F7.2,A,F7.2,A)') "->" , "Point x2 [m]: ","[", & + point_2(1,1), ", ", point_2(2,1), ", ", point_2(3,1), "]" + write(STD_OUT,'(30X,A,A27,A,F7.2,A,F7.2,A,F7.2,A)') "->" , "Point x3 [m]: ","[", & + point_3(1,1), ", ", point_3(2,1), ", ", point_3(3,1), "]" + else + write(STD_OUT,'(30X,A,A27,A27)') "->" , "Input Type: " , "Specified Input File" + write(STD_OUT,'(30X,A,A27,A27)') "->" , "Input File: " , self % fileInput + write(STD_OUT,'(30X,A,A27,A,F7.2,A,F7.2,A,F7.2,A)') "->" , "First Point [m]: ","[", & + self % x(1,1), ", ", self % x(2,1), ", ", self % x(3,1), "]" + write(STD_OUT,'(30X,A,A27,A,F7.2,A,F7.2,A,F7.2,A)') "->" , "Last Point [m]: ","[", & + self % x(1,self % nNodes), ", ", self % x(2,self % nNodes), ", ", self % x(3,self % nNodes), "]" + end if + write(STD_OUT,'(30X,A,A27,A,I4,A,I4,A)') "->" , "Nodes [N12,N23]: ","[", & + self % N(1), ", ", self % N(2), "]" + write(STD_OUT,'(30X,A,A27,I4)') "->" , "Samplings Interval: ", self % interval + write(STD_OUT,'(30X,A,A27,A128)') "->" , "Filename: ", fileFormat + if (firstOutDomain) then + write(STD_OUT,'(30X,A,A27,A61)') "->" ,"Note: ","Node/Nodes are located out of domain. A NaN will be assigned" + end if + + end subroutine Plane_Initialization +! +! Update Lagrange Interpolants After pAdaptation +! ----------------------------------------------------- + subroutine Plane_UpdateLagrangeInterp(self, mesh) + use Physics + use MPI_Process_Info + use, intrinsic :: ieee_arithmetic, only: IEEE_Value, IEEE_QUIET_NAN + implicit none + class(PlaneSampling_t) :: self + class(HexMesh) :: mesh + +! +! --------------- +! Local variables +! --------------- +! + integer :: i +! +! Find node location in the element ID +! ------------------------------------ + DO i=1,self % nNodes +! +! Get the Lagrange interpolants +! ----------------------------- + if ( (MPI_Process % rank .ne. self % rank (i)).OR.( .not. self % active (i) ) ) then + self % lxi (:,i)= 0.0_RP + self % leta (:,i)= 0.0_RP + self % lzeta (:,i)= 0.0_RP + ELSE + associate(e => mesh % elements(self % eID(i))) + associate( spAxi => NodalStorage(e % Nxyz(1)), & + spAeta => NodalStorage(e % Nxyz(2)), & + spAzeta => NodalStorage(e % Nxyz(3)) ) + self % lxi (:,i) = spAxi % lj(self % xi(1,i)) + self % leta (:,i) = spAeta % lj(self % xi(2,i)) + self % lzeta (:,i)= spAzeta % lj(self % xi(3,i)) + end associate + end associate + END IF + + END DO + end subroutine Plane_UpdateLagrangeInterp + + subroutine Plane_Update(self, mesh, bufferPosition, t) + use Physics + use MPI_Process_Info + use, intrinsic :: ieee_arithmetic, only: IEEE_Value, IEEE_QUIET_NAN + implicit none + class(PlaneSampling_t) :: self + class(HexMesh) :: mesh + integer :: bufferPosition + real(kind=RP) :: t + +! +! --------------- +! Local variables +! --------------- +! + integer :: i, j, k, l, m, ierr + real(kind=RP) :: value, kappa + real(kind=RP) :: c, sqrtRho + real(kind=RP) :: NaN + real(kind=RP) :: Sym, Asym + real(kind=RP) :: invRhoBase, invRhoF + real(kind=RP), allocatable :: var(:,:,:) + + if (self % intervalCount .EQ. 0 ) then + self % bufferLine = self % bufferLine + 1 + + DO m=1,self % nVariables + + self % values(1, bufferPosition, m) = t + + DO l=1, self % nNodes +! +! Assign NaN if the node located outside of domain +! ------------------------------------------------ + IF ( .not. self % active (l) ) THEN + self % values(l, bufferPosition,m) = IEEE_VALUE(NaN, IEEE_QUIET_NAN) + ELSE + if ( MPI_Process % rank .eq. self % rank (l) ) then +! +! Update the Node +! ---------------- + associate( e => mesh % elements(self % eID(l)) ) + allocate (var(0:e % Nxyz(1),0:e % Nxyz(2),0:e % Nxyz(3) )) + associate( Q => e % storage % Q, S => e % storage ) + +#ifdef NAVIERSTOKES + select case (trim(self % variable(m))) + case("density") + do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1) + var(i,j,k) = Q(IRHO,i,j,k) + end do ; end do ; end do + case("pressure") + do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1) + var(i,j,k) = Pressure(Q(:,i,j,k)) + end do ; end do ; end do + + case("viscosity") + do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1) + call get_laminar_mu_kappa(Q(:,i,j,k),var(i,j,k),kappa) + end do ; end do ; end do + + case("ptotal") + do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1) + var(i,j,k) = POW2(Q(IRHOU,i,j,k)) + POW2(Q(IRHOV,i,j,k)) + POW2(Q(IRHOW,i,j,k))/POW2(Q(IRHO,i,j,k)) ! Vabs**2 + var(i,j,k) = var(i,j,k) / ( thermodynamics % gamma*(thermodynamics % gamma-1.0_RP)*(Q(IRHOE,i,j,k)/Q(IRHO,i,j,k)-0.5_RP * var(i,j,k)) ) ! Mach ^2 + var(i,j,k) = Pressure(Q(:,i,j,k))*(1.0_RP+0.5_RP*(thermodynamics % gamma-1.0_RP)*var(i,j,k))**( thermodynamics % gamma/(thermodynamics % gamma-1.0_RP)) + end do ; end do ; end do + + case("velocity") + do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1) + var(i,j,k) = sqrt(POW2(Q(IRHOU,i,j,k)) + POW2(Q(IRHOV,i,j,k)) + POW2(Q(IRHOW,i,j,k)))/Q(IRHO,i,j,k) + end do ; end do ; end do + + case("omegax") + do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1) + var(i,j,k) = (1/Q(IRHO,i,j,k) * S % U_y(IRHOW,i,j,k) - Q(IRHOW,i,j,k)/(Q(IRHO,i,j,k)**2) * S % U_y(IRHO,i,j,k)) & + - (1/Q(IRHO,i,j,k) * S % U_z(IRHOV,i,j,k) - Q(IRHOV,i,j,k)/(Q(IRHO,i,j,k)**2) * S % U_z(IRHO,i,j,k)) + end do ; end do ; end do + + case("omegay") + do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1) + var(i,j,k) = (1/Q(IRHO,i,j,k) * S % U_z(IRHOU,i,j,k) - Q(IRHOU,i,j,k)/(Q(IRHO,i,j,k)**2) * S % U_z(IRHO,i,j,k)) & + - (1/Q(IRHO,i,j,k) * S % U_x(IRHOW,i,j,k) - Q(IRHOW,i,j,k)/(Q(IRHO,i,j,k)**2) * S % U_x(IRHO,i,j,k)) + end do ; end do ; end do + + case("omegaz") + do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1) + var(i,j,k) = (1/Q(IRHO,i,j,k) * S % U_x(IRHOV,i,j,k) - Q(IRHOV,i,j,k)/(Q(IRHO,i,j,k)**2) * S % U_x(IRHO,i,j,k)) & + - (1/Q(IRHO,i,j,k) * S % U_y(IRHOU,i,j,k) - Q(IRHOU,i,j,k)/(Q(IRHO,i,j,k)**2) * S % U_y(IRHO,i,j,k)) + end do ; end do ; end do + + case("u") + do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1) + var(i,j,k) = Q(IRHOU,i,j,k) / Q(IRHO,i,j,k) + end do ; end do ; end do + + case("v") + do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1) + var(i,j,k) = Q(IRHOV,i,j,k) / Q(IRHO,i,j,k) + end do ; end do ; end do + + case("w") + do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1) + var(i,j,k) = Q(IRHOW,i,j,k) / Q(IRHO,i,j,k) + end do ; end do ; end do + + case("mach") + do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1) + var(i,j,k) = POW2(Q(IRHOU,i,j,k)) + POW2(Q(IRHOV,i,j,k)) + POW2(Q(IRHOW,i,j,k))/POW2(Q(IRHO,i,j,k)) ! Vabs**2 + var(i,j,k) = sqrt( var(i,j,k) / ( thermodynamics % gamma*(thermodynamics % gamma-1.0_RP)*(Q(IRHOE,i,j,k)/Q(IRHO,i,j,k)-0.5_RP * var(i,j,k)) ) ) + end do ; end do ; end do + + case("k") + do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1) + var(i,j,k) = 0.5_RP * (POW2(Q(IRHOU,i,j,k)) + POW2(Q(IRHOV,i,j,k)) + POW2(Q(IRHOW,i,j,k)))/Q(IRHO,i,j,k) + end do ; end do ; end do + + case("q1") + do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1) + var(i,j,k) = Q(IRHO,i,j,k) + end do ; end do ; end do + + case("q2") + do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1) + var(i,j,k) = Q(IRHOU,i,j,k) + end do ; end do ; end do + + case("q3") + do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1) + var(i,j,k) = Q(IRHOV,i,j,k) + end do ; end do ; end do + + case("q4") + do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1) + var(i,j,k) = Q(IRHOW,i,j,k) + end do ; end do ; end do + + case("q5") + do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1) + var(i,j,k) = Q(IRHOE,i,j,k) + end do ; end do ; end do + + end select +#endif +#ifdef MULTIPHASE + select case (trim(self % variable(m))) + case("static-pressure") + do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1) + var(i,j,k) = Q(IMP,i,j,k) + Q(IMC,i,j,k)*e % storage % mu(1,i,j,k) - 12.0_RP*multiphase%sigma*multiphase%invEps*(POW2(Q(IMC,i,j,k)*(1.0_RP-Q(IMC,i,j,k)))) & + - 0.25_RP*3.0_RP*multiphase % sigma * multiphase % eps * (POW2(e % storage % c_x(1,i,j,k))+POW2(e % storage % c_y(1,i,j,k))+POW2(e % storage % c_z(1,i,j,k))) + end do ; end do ; end do + + case("concentration") + do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1) + var(i,j,k) = min(max(Q(1,i,j,k),0.0_RP),1.0_RP) + end do ; end do ; end do + + case("density") + do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1) + c = min(max(Q(1,i,j,k),0.0_RP),1.0_RP) + var(i,j,k) = c * dimensionless % rho (1) + (1.0_RP-c) * dimensionless % rho(2) + end do ; end do ; end do + + case("pressure") + do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1) + var(i,j,k) = Q(5,i,j,k) + end do ; end do ; end do + + case("velocity") + do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1) + c = min(max(Q(1,i,j,k),0.0_RP),1.0_RP) + sqrtRho = sqrt(c * dimensionless % rho (1) + (1.0_RP-c) * dimensionless % rho(2)) + var(i,j,k) = sqrt(POW2(Q(2,i,j,k)) + POW2(Q(3,i,j,k)) + POW2(Q(4,i,j,k)))/sqrtRho + end do ; end do ; end do + + case("omegax") + do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1) + c = min(max(Q(1,i,j,k),0.0_RP),1.0_RP) + sqrtRho = sqrt(c * dimensionless % rho (1) + (1.0_RP-c) * dimensionless % rho(2)) + var(i,j,k) = (1/sqrtRho * S % U_y(IMSQRHOW,i,j,k) - Q(IMSQRHOW,i,j,k)/(sqrtRho**2) * S % U_y(IMC,i,j,k)) & + - (1/sqrtRho * S % U_z(IMSQRHOV,i,j,k) - Q(IMSQRHOV,i,j,k)/(sqrtRho**2) * S % U_z(IMC,i,j,k)) + end do ; end do ; end do + + case("omegay") + do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1) + c = min(max(Q(1,i,j,k),0.0_RP),1.0_RP) + sqrtRho = sqrt(c * dimensionless % rho (1) + (1.0_RP-c) * dimensionless % rho(2)) + var(i,j,k) = (1/sqrtRho * S % U_z(IMSQRHOU,i,j,k) - Q(IMSQRHOU,i,j,k)/(sqrtRho**2) * S % U_z(IMC,i,j,k)) & + - (1/sqrtRho * S % U_x(IMSQRHOW,i,j,k) - Q(IMSQRHOW,i,j,k)/(sqrtRho**2) * S % U_x(IMC,i,j,k)) + end do ; end do ; end do + + case("omegaz") + do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1) + c = min(max(Q(1,i,j,k),0.0_RP),1.0_RP) + sqrtRho = sqrt(c * dimensionless % rho (1) + (1.0_RP-c) * dimensionless % rho(2)) + var(i,j,k) = (1/sqrtRho * S % U_x(IMSQRHOV,i,j,k) - Q(IMSQRHOV,i,j,k)/(sqrtRho**2) * S % U_x(IMC,i,j,k)) & + - (1/sqrtRho * S % U_y(IMSQRHOU,i,j,k) - Q(IMSQRHOU,i,j,k)/(sqrtRho**2) * S % U_y(IMC,i,j,k)) + end do ; end do ; end do + + case("u") + do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1) + c = min(max(Q(1,i,j,k),0.0_RP),1.0_RP) + sqrtRho = sqrt(c * dimensionless % rho (1) + (1.0_RP-c) * dimensionless % rho(2)) + var(i,j,k) = Q(IMSQRHOU,i,j,k) / sqrtRho + end do ; end do ; end do + + case("v") + do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1) + c = min(max(Q(1,i,j,k),0.0_RP),1.0_RP) + sqrtRho = sqrt(c * dimensionless % rho (1) + (1.0_RP-c) * dimensionless % rho(2)) + var(i,j,k) = Q(IMSQRHOV,i,j,k) / sqrtRho + end do ; end do ; end do + + case("w") + do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1) + c = min(max(Q(1,i,j,k),0.0_RP),1.0_RP) + sqrtRho = sqrt(c * dimensionless % rho (1) + (1.0_RP-c) * dimensionless % rho(2)) + var(i,j,k) = Q(IMSQRHOW,i,j,k) / sqrtRho + end do ; end do ; end do + + case("q1") + do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1) + var(i,j,k) = Q(IMC,i,j,k) + end do ; end do ; end do + + case("q2") + do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1) + var(i,j,k) = Q(IMSQRHOU,i,j,k) + end do ; end do ; end do + + case("q3") + do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1) + var(i,j,k) = Q(IMSQRHOV,i,j,k) + end do ; end do ; end do + + case("q4") + do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1) + var(i,j,k) = Q(IMSQRHOW,i,j,k) + end do ; end do ; end do + + case("q5") + do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1) + var(i,j,k) = Q(IMP,i,j,k) + end do ; end do ; end do + + end select +#endif + + value = 0.0_RP + do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1) + value = value + var(i,j,k) * self % lxi(i,l) * self % leta(j,l) * self % lzeta(k,l) + end do ; end do ; end do + + self % values(l+1, bufferPosition, m) = value + + end associate + deallocate(var) + end associate +#ifdef _HAS_MPI_ + if ( MPI_Process % doMPIAction ) then + +! Share the result with the rest of the processes +! ----------------------------------------------- + call mpi_bcast(value, 1, MPI_DOUBLE, self % rank (l), MPI_COMM_WORLD, ierr) + end if +#endif + else + +! Receive the result from the rank that contains the Plane +! -------------------------------------------------------- +#ifdef _HAS_MPI_ + if ( MPI_Process % doMPIAction ) then + call mpi_bcast(self % values(l+1, bufferPosition, m), 1, MPI_DOUBLE, self % rank (l), MPI_COMM_WORLD, ierr) + end if +#endif + end if + END IF + END DO + END DO + end if + self % intervalCount = self % intervalCount + 1 + + + if (self % intervalCount .EQ. self % interval) then + self % intervalCount = 0 + end if + + end subroutine Plane_Update + + subroutine Plane_WriteToFile ( self, no_of_lines) +! +! ************************************************************* +! This subroutine writes the buffer to the file. +! ************************************************************* +! + implicit none + class(PlaneSampling_t) :: self + integer :: no_of_lines +! +! --------------- +! Local variables +! --------------- +! + integer :: i, j, ierr + integer :: fID + + if ( MPI_Process % isRoot ) then + DO j=1, self % nVariables + open( newunit = fID , file = trim ( self % fileName (j) ) , action = "write" , access = "stream" , status = "old", position='append' ) + do i = 1 , no_of_lines + write( fID ) self % values(:,i,j) + end do + close ( fID ) + END DO + end if + if ( no_of_lines .ne. 0 ) self % values(:,1,:) = self % values(:,no_of_lines,:) + + self % bufferLine = 0 + + end subroutine Plane_WriteToFile + + subroutine Plane_LookInOtherPartitions(self, i) + use MPI_Process_Info + implicit none + class(PlaneSampling_t) :: self + integer ,intent(in) :: i + integer :: allActives(MPI_Process % nProcs) + integer :: active, ierr + + if ( MPI_Process % doMPIAction ) then +#ifdef _HAS_MPI_ +! +! Cast the logicals onto integers +! ------------------------------- + if ( self % active (i) ) then + active = 1 + else + active = 0 + end if + allActives=0 +! +! Gather all data from all processes +! ---------------------------------- + call mpi_allgather(active, 1, MPI_INT, allActives, 1, MPI_INT, MPI_COMM_WORLD, ierr) +! +! Check if any of them found the Plane +! ------------------------------------ + if ( any(allActives .eq. 1) ) then +! +! Assign the domain of the partition that contains the Plane +! ---------------------------------------------------------- + self % active (i) = .true. + self % rank (i) = maxloc(allActives, dim = 1) - 1 + + else +! +! Disable the Plane +! ----------------- + self % active (i) = .false. + self % rank (i) = -1 + self % eID (i) = 1 + + end if +#endif + else +! +! Without MPI select the rank 0 as default +! ---------------------------------------- + self % rank (i)= 0 + + end if + + end subroutine Plane_LookInOtherPartitions + + elemental subroutine Plane_Destruct (self) + implicit none + class(PlaneSampling_t), intent(inout) :: self + + safedeallocate (self % values) + safedeallocate (self % lxi) + safedeallocate (self % leta) + safedeallocate (self % lzeta) + safedeallocate (self % active) + safedeallocate (self % rank) + safedeallocate (self % eID) + safedeallocate (self % x) + safedeallocate (self % xi) + safedeallocate (self % fileName) + safedeallocate (self % variable) + + end subroutine Plane_Destruct + + elemental subroutine Plane_Assign (to, from) + implicit none + class(PlaneSampling_t), intent(inout) :: to + type(PlaneSampling_t) , intent(in) :: from + + safedeallocate ( to % active ) + allocate ( to % active ( size(from % active) ) ) + to % active = from % active + + safedeallocate ( to % rank ) + allocate ( to % rank ( size(from % rank) ) ) + to % rank = from % rank + + to % ID = from % ID + to % nVariables = from % nVariables + to % interval = from % interval + to % bufferSize = from % bufferSize + to % bufferLine = from % bufferLine + to % intervalCount = from % intervalCount + to % N = from % N + to % nNodes = from % nNodes + to % disturbanceData = from % disturbanceData + + safedeallocate ( to % eID ) + allocate ( to % eID ( size(from % eID) ) ) + to % eID = from % eID + + safedeallocate ( to % lxi ) + allocate ( to % lxi ( size(from % lxi,1), size(from % lxi,2) ) ) + to % lxi = from % lxi + + safedeallocate ( to % leta ) + allocate ( to % leta ( size(from % leta,1), size(from % leta,2) ) ) + to % leta = from % leta + + safedeallocate ( to % lzeta ) + allocate ( to % lzeta ( size(from % lzeta,1), size(from % lzeta,2) ) ) + to % lzeta = from % lzeta + + safedeallocate ( to % values ) + allocate ( to % values( size(from % values,1),size(from % values,2),size(from % values,3) ) ) + to % values = from % values + + safedeallocate ( to % x ) + allocate ( to % x( size(from % x,1),size(from % x,2) ) ) + to % x = from % x + + safedeallocate ( to % xi ) + allocate ( to % xi( size(from % xi,1),size(from % xi,2) ) ) + to % xi = from % xi + + safedeallocate ( to % fileName ) + allocate ( to % fileName ( size(from % fileName) ) ) + to % fileName = from % fileName + + to % planeName = from % planeName + + safedeallocate ( to % variable ) + allocate ( to % variable ( size(from % variable) ) ) + to % variable = from % variable + + + + end subroutine Plane_Assign + +end module PlaneSampling + diff --git a/Solver/src/libs/monitors/SamplingOperator.f90 b/Solver/src/libs/monitors/SamplingOperator.f90 new file mode 100644 index 0000000000..1072ab147c --- /dev/null +++ b/Solver/src/libs/monitors/SamplingOperator.f90 @@ -0,0 +1,295 @@ +#include "Includes.h" + +module SamplingOperator + use SMConstants + use PhysicsStorage + use Physics + use FaceClass + use ElementClass + use HexMeshClass + use FluidData + use VariableConversion + use NodalStorageClass +#ifdef _HAS_MPI_ + use mpi +#endif + implicit none + + private + public SHEAR_STRESS_TANGENT, SHEAR_STRESS_X, SHEAR_STRESS_Y, SHEAR_STRESS_Z, PRESSURE_SURF, Q1, Q2, Q3, Q4, Q5 + public VectorSurfaceSampling + + integer, parameter :: SHEAR_STRESS_TANGENT = 1 + integer, parameter :: SHEAR_STRESS_X = 2 + integer, parameter :: SHEAR_STRESS_Y = 3 + integer, parameter :: SHEAR_STRESS_Z = 4 + integer, parameter :: PRESSURE_SURF = 5 + integer, parameter :: Q1 = 6 + integer, parameter :: Q2 = 7 + integer, parameter :: Q3 = 8 + integer, parameter :: Q4 = 9 + integer, parameter :: Q5 = 10 + integer, parameter :: USER_DEFINED = 99 +! +! ======== + contains +! ======== +! +! SUBROUTINE TO GET CONSTRUCT THE DATA +!//////////////////////////////////////////////////////////////////////////////////////// +! + subroutine VectorSurfaceSampling(mesh, zoneID, integralType, monitorName, data_out) + use MonitorDefinitions +#ifdef _HAS_MPI_ + use mpi +#endif + implicit none + class(HexMesh), intent(inout), target :: mesh + integer, intent(in) :: zoneID + integer, intent(in) :: integralType + real(kind=RP), allocatable, intent(out) :: data_out(:) + +! +! --------------- +! Local variables +! --------------- +! + integer :: zonefID, fID, eID, fIDs(6), ierr, Nx, Ny, fsID + class(Element), pointer :: elements(:) + real(kind=RP) , allocatable :: data_proc(:) + logical :: file_exists + character(len=STR_LEN_MONITORS) :: fileName + character(len=STR_LEN_MONITORS) :: monitorName + + ! +! Get the number of Order and the sampling interval +! ------------------------------------------------- + if (mesh % zones(zoneID) % no_of_faces .gt.0) then + Nx = mesh % faces (mesh % zones(zoneID) % faces(1))%Nf(1)+1 + Ny = mesh % faces (mesh % zones(zoneID) % faces(1))%Nf(2)+1 + else + Nx=0 + Ny=0 + end if + + ALLOCATE (data_out(Nx*Ny*mesh % zones(zoneID) % no_of_faces), data_proc(Nx*Ny)) + +! +! ************************* +! Perform the interpolation +! ************************* +! +#if defined(NAVIERSTOKES) && (!(INCNS)) + elements => mesh % elements +!$omp parallel private(fID, eID, fIDs,data_proc) shared(elements,mesh,NodalStorage,zoneID,integralType,& +!$omp& computeGradients,data_out,Nx) +!$omp single + do zonefID = 1, mesh % zones(zoneID) % no_of_faces + fID = mesh % zones(zoneID) % faces(zonefID) + + eID = mesh % faces(fID) % elementIDs(1) + fIDs = mesh % elements(eID) % faceIDs + +!$omp task depend(inout:elements(eID)) + call elements(eID) % ProlongSolutionToFaces(NCONS, mesh % faces(fIDs(1)),& + mesh % faces(fIDs(2)),& + mesh % faces(fIDs(3)),& + mesh % faces(fIDs(4)),& + mesh % faces(fIDs(5)),& + mesh % faces(fIDs(6)) ) + if ( computeGradients ) then + call elements(eID) % ProlongGradientsToFaces(NGRAD, mesh % faces(fIDs(1)),& + mesh % faces(fIDs(2)),& + mesh % faces(fIDs(3)),& + mesh % faces(fIDs(4)),& + mesh % faces(fIDs(5)),& + mesh % faces(fIDs(6)) ) + end if +!$omp end task + end do +!$omp end single + +! +! Loop the zone to get faces and elements +! --------------------------------------- +!$omp do private(fID,data_proc) schedule(runtime) + do zonefID = 1, mesh % zones(zoneID) % no_of_faces +! +! Face global ID +! -------------- + fID = mesh % zones(zoneID) % faces(zonefID) +! +! Compute the integral +! -------------------- + + CALL VectorSurfaceSampling_Face(mesh % faces(fID), integralType, fID, data_proc) + + data_out((zonefID-1)*Nx*Ny+1:zonefID*Nx*Ny)=data_proc + + end do + +!$omp end do +!$omp end parallel + +#endif + + end subroutine VectorSurfaceSampling +! +! SUBROUTINE TO GET CONSTRUCT THE DATA in a FACE +!//////////////////////////////////////////////////////////////////////////////////////// +! + subroutine VectorSurfaceSampling_Face(f, integralType, fID, data_out) + implicit none + class(Face), intent(in) :: f + integer, intent(in) :: integralType, fID + real(kind=RP) ,intent(out) :: data_out(1:((f%Nf(1)+1)*(f%Nf(2)+1))) +! +! --------------- +! Local variables +! --------------- +! + integer :: i, j ! Face indices + integer :: k + real(kind=RP) :: p, tau(NDIM,NDIM) + type(NodalStorage_t), pointer :: spAxi, spAeta + real(kind=RP) :: shear(NDIM) + real(kind=RP) :: shearTangent + +! +! Initialization +! -------------- + spAxi => NodalStorage(f % Nf(1)) + spAeta => NodalStorage(f % Nf(2)) +! +! Perform the numerical integration +! --------------------------------- + associate( Q => f % storage(1) % Q, & + U_x => f % storage(1) % U_x, & + U_y => f % storage(1) % U_y, & + U_z => f % storage(1) % U_z ) + select case ( integralType ) +#if defined(NAVIERSTOKES) && (!(INCNS)) +! +! ************************************************* +! Computes the shear stress, tangent to the surface +! ************************************************* +! + case ( SHEAR_STRESS_TANGENT ) + k=0 + do j = 0, f % Nf(2); do i = 0, f % Nf(1) + k=k+1 +! +! Get the stress tensor and the tangent component +! ---------------------------------------------- + call getStressTensor(Q(:,i,j),U_x(:,i,j),U_y(:,i,j),U_z(:,i,j), tau) + shear=- matmul(tau,f % geom % normal(:,i,j))* POW2(refValues % V) *refValues % rho + shearTangent=dot_product(shear,f % geom % t1(:,i,j)) + data_out(k)=shearTangent + end do ; end do +! +! ********************************************* +! Computes the shear stress, in the x direction +! ********************************************* +! + case ( SHEAR_STRESS_X ) + k=0 + do j = 0, f % Nf(2); do i = 0, f % Nf(1) + k=k+1 +! +! Get the stress tensor and the tangent component +! ---------------------------------------------- + call getStressTensor(Q(:,i,j),U_x(:,i,j),U_y(:,i,j),U_z(:,i,j), tau) + shear=- matmul(tau,f % geom % normal(:,i,j))* POW2(refValues % V) *refValues % rho + data_out(k)=shear(1) + end do ; end do +! +! ********************************************* +! Computes the shear stress, in the x direction +! ********************************************* +! + case ( SHEAR_STRESS_Y ) + k=0 + do j = 0, f % Nf(2); do i = 0, f % Nf(1) + k=k+1 +! +! Get the stress tensor +! --------------------- + call getStressTensor(Q(:,i,j),U_x(:,i,j),U_y(:,i,j),U_z(:,i,j), tau) + shear=- matmul(tau,f % geom % normal(:,i,j))* POW2(refValues % V) *refValues % rho + data_out(k)=shear(2) + end do ; end do +! +! ********************************************* +! Computes the shear stress, in the x direction +! ********************************************* +! + case ( SHEAR_STRESS_Z ) + k=0 + do j = 0, f % Nf(2); do i = 0, f % Nf(1) + k=k+1 +! +! Get the stress tensor +! --------------------- + call getStressTensor(Q(:,i,j),U_x(:,i,j),U_y(:,i,j),U_z(:,i,j), tau) + shear=- matmul(tau,f % geom % normal(:,i,j))* POW2(refValues % V) *refValues % rho + data_out(k)=shear(3) + end do ; end do +! +! ********************************************* +! Computes the shear stress, in the x direction +! ********************************************* +! + case ( PRESSURE_SURF ) + k=0 + do j = 0, f % Nf(2); do i = 0, f % Nf(1) + k=k+1 +! +! Get the pressure +! --------------------- + data_out(k)=Pressure(Q(:,i,j))* POW2(refValues % V) *refValues % rho + end do ; end do + + case ( Q1 ) + k=0 + do j = 0, f % Nf(2); do i = 0, f % Nf(1) + k=k+1 + data_out(k)=Q(1,i,j) + end do ; end do + + case ( Q2 ) + k=0 + do j = 0, f % Nf(2); do i = 0, f % Nf(1) + k=k+1 + data_out(k)=Q(2,i,j) + end do ; end do + + case ( Q3 ) + k=0 + do j = 0, f % Nf(2); do i = 0, f % Nf(1) + k=k+1 + data_out(k)=Q(3,i,j) + end do ; end do + + case ( Q4 ) + k=0 + do j = 0, f % Nf(2); do i = 0, f % Nf(1) + k=k+1 + data_out(k)=Q(4,i,j) + end do ; end do + + case ( Q5 ) + k=0 + do j = 0, f % Nf(2); do i = 0, f % Nf(1) + k=k+1 + data_out(k)=Q(5,i,j) + end do ; end do +#endif + end select + end associate + nullify (spAxi, spAeta) + end subroutine VectorSurfaceSampling_Face + + +end module SamplingOperator + + diff --git a/Solver/src/libs/monitors/Samplings.f90 b/Solver/src/libs/monitors/Samplings.f90 new file mode 100644 index 0000000000..e69d7c41ac --- /dev/null +++ b/Solver/src/libs/monitors/Samplings.f90 @@ -0,0 +1,390 @@ +#include "Includes.h" +module Samplings + use SMConstants + use NodalStorageClass + use HexMeshClass + use MonitorDefinitions + use FileReadingUtilities , only: getFileName + use SurfaceSampling + use PlaneSampling + use SpatialMeanNode + implicit none +! + + private + public Sampling_t +! +! ***************************** +! Main sampling class definition +! ***************************** +! + type Sampling_t + character(len=LINE_LENGTH) :: solution_file + integer :: no_of_surfaceSamplings =0 + integer :: no_of_planeSamplings =0 + integer :: no_of_spatialMeanNodes =0 + integer :: dt_restriction + logical :: write_dt_restriction + class(SurfaceSampling_t), allocatable :: surfaceSamplings(:) + class(PlaneSampling_t), allocatable :: planeSamplings(:) + class(SpatialMeanNode_t), allocatable :: spatialMeanNodes(:) + contains + procedure :: Construct => Samplings_Construct + procedure :: UpdateInterp => Samplings_UpdateLagrangeInterp + procedure :: UpdateValues => Sampling_UpdateValues + procedure :: WriteToFile => Sampling_WriteToFile + procedure :: destruct => Sampling_Destruct + procedure :: copy => Sampling_Assign + generic :: assignment(=) => copy + end type Sampling_t +! +! ======== + contains +! ======== +! +!/////////////////////////////////////////////////////////////////////////////////////// +! + subroutine Samplings_Construct( Samplings, mesh, controlVariables ) + use Headers + use FTValueDictionaryClass + use mainKeywordsModule + implicit none + class(Sampling_t) :: Samplings + class(HexMesh), intent(in) :: mesh + class(FTValueDictionary), intent(in) :: controlVariables + +! +! --------------- +! Local variables +! --------------- +! + integer :: fID , io + integer :: i + character(len=STR_LEN_MONITORS) :: line + character(len=STR_LEN_MONITORS) :: solution_file + logical, save :: FirstCall = .TRUE. +! +! Setup the buffer +! ---------------- + if (controlVariables % containsKey("monitors flush interval") ) then + BUFFER_SIZE = controlVariables % integerValueForKey("monitors flush interval") + end if +! +! Get the solution file name +! -------------------------- + solution_file = controlVariables % stringValueForKey( solutionFileNameKey, requestedLength = STR_LEN_MONITORS ) +! +! Remove the *.hsol termination +! ----------------------------- + solution_file = trim(getFileName(solution_file)) + Samplings % solution_file = trim(solution_file) +! +! Search in case file for samplings +! --------------------------------------------------------------------- + if (mesh % child) then ! Return doing nothing if this is a child mesh + Samplings % no_of_surfaceSamplings = 0 + Samplings % no_of_planeSamplings = 0 + Samplings % no_of_spatialMeanNodes = 0 + else + call getNoOfSamplings( Samplings % no_of_surfaceSamplings, Samplings % no_of_planeSamplings, Samplings % no_of_spatialMeanNodes) + end if +! +! Initialize +! ---------- + + + allocate ( Samplings % surfaceSamplings ( Samplings % no_of_surfaceSamplings ) ) + allocate ( Samplings % planeSamplings ( Samplings % no_of_planeSamplings ) ) + allocate ( Samplings % spatialMeanNodes ( Samplings % no_of_spatialMeanNodes ) ) + + if (Samplings % no_of_surfaceSamplings .GT. 0) then + call Section_Header("Initialize Surface Samplings") + end if + do i = 1 , Samplings % no_of_surfaceSamplings + call Samplings % surfaceSamplings(i) % Initialization ( mesh , i, solution_file , FirstCall ) + end do + + if (Samplings % no_of_planeSamplings .GT. 0) then + call Section_Header("Initialize Plane Samplings") + end if + do i = 1 , Samplings % no_of_planeSamplings + call Samplings % planeSamplings(i) % Initialization ( mesh , i, solution_file , FirstCall ) + end do + + if (Samplings % no_of_spatialMeanNodes .GT. 0) then + call Section_Header("Initialize Spatial Mean Node") + end if + do i = 1 , Samplings % no_of_spatialMeanNodes + call Samplings % spatialMeanNodes(i) % Initialization ( mesh , i, solution_file , FirstCall ) + end do + + FirstCall = .FALSE. + + end subroutine Samplings_Construct + + subroutine Samplings_UpdateLagrangeInterp (self, mesh) +! +! ******************************************************************* +! This subroutine updates the Lagrange interpolants after pAdaptation +! ******************************************************************* +! + implicit none + class(Sampling_t) :: self + class(HexMesh) :: mesh +! +! --------------- +! Local variables +! --------------- +! + integer :: i + +! +! Update interpolants plane Samplings +! ----------------------------------- + do i = 1 , self % no_of_planeSamplings + call self % planeSamplings(i) % UpdateInterp( mesh ) + end do + + end subroutine Samplings_UpdateLagrangeInterp + + subroutine Sampling_UpdateValues (self, mesh, t) +! +! *************************************************************** +! This subroutine updates the values for the Samplings. +! *************************************************************** +! + use PhysicsStorage + use StopwatchClass + implicit none + class(Sampling_t) :: self + class(HexMesh) :: mesh + real(kind=RP) :: t +! +! --------------- +! Local variables +! --------------- +! + integer :: i + +! +! Update surface Samplings +! ------------------------ + do i = 1 , self % no_of_surfaceSamplings + call self % surfaceSamplings(i) % Update( mesh , self % surfaceSamplings(i) % bufferLine, t ) + end do +! +! Update plane Samplings +! ---------------------- + do i = 1 , self % no_of_planeSamplings + call self % planeSamplings(i) % Update( mesh , self % planeSamplings(i) % bufferLine, t ) + end do +! +! Update spatial mean node +! ------------------------ + do i = 1 , self % no_of_spatialMeanNodes + call self % spatialMeanNodes(i) % Update( mesh , self % spatialMeanNodes(i) % bufferLine, t ) + end do + + end subroutine Sampling_UpdateValues + + subroutine Sampling_WriteToFile ( self , mesh, force) +! +! ****************************************************************** +! This routine has a double behaviour: +! force = .true. -> Writes to file and resets buffers +! force = .false. -> Just writes to file if the buffer is full +! ****************************************************************** +! + use MPI_Process_Info + implicit none + class(Sampling_t) :: self + class(HexMesh) :: mesh + logical, optional :: force +! ------------------------------------------------ + integer :: i + logical :: forceVal + if ( present ( force ) ) then + forceVal = force + + else + forceVal = .false. + + end if + + if ( forceVal ) then +! +! In this case the Samplings are exported to their files and the buffer is reseted +! ------------------------------------------------------------------------------- + + do i = 1 , self % no_of_surfaceSamplings + call self % surfaceSamplings(i) % WriteToFile ( self % surfaceSamplings(i) % bufferLine ) + end do + + do i = 1 , self % no_of_planeSamplings + call self % planeSamplings(i) % WriteToFile ( self % planeSamplings(i) % bufferLine ) + end do + + do i = 1 , self % no_of_spatialMeanNodes + call self % spatialMeanNodes(i) % WriteToFile ( self % spatialMeanNodes(i) % bufferLine ) + end do + + else +! +! The Samplings are exported just if the buffer is full +! ---------------------------------------------------- + + do i = 1 , self % no_of_surfaceSamplings + if ( self % surfaceSamplings(i) % bufferLine .eq. self % surfaceSamplings(i) % bufferSize) then + call self % surfaceSamplings(i) % WriteToFile ( self % surfaceSamplings(i) % bufferLine ) + end if + end do + do i = 1 , self % no_of_planeSamplings + if ( self % planeSamplings(i) % bufferLine .eq. self % planeSamplings(i) % bufferSize) then + call self % planeSamplings(i) % WriteToFile ( self % planeSamplings(i) % bufferLine ) + end if + end do + + do i = 1 , self % no_of_spatialMeanNodes + if ( self % spatialMeanNodes(i) % bufferLine .eq. self % spatialMeanNodes(i) % bufferSize) then + call self % spatialMeanNodes(i) % WriteToFile ( self % spatialMeanNodes(i) % bufferLine ) + end if + end do + + end if + end subroutine Sampling_WriteToFile + + subroutine Sampling_Destruct (self) + implicit none + class(Sampling_t) :: self + + if ( self % no_of_surfaceSamplings .gt. 0) call self % surfaceSamplings % destruct + if ( self % no_of_planeSamplings .gt. 0) call self % planeSamplings % destruct + if ( self % no_of_spatialMeanNodes .gt. 0) call self % spatialMeanNodes % destruct + + safedeallocate (self % surfaceSamplings) + safedeallocate (self % planeSamplings) + safedeallocate (self % spatialMeanNodes) + + end subroutine + + impure elemental subroutine Sampling_Assign ( to, from ) + implicit none + !-arguments-------------------------------------- + class(Sampling_t), intent(inout) :: to + type(Sampling_t) , intent(in) :: from + !-local-variables-------------------------------- + !------------------------------------------------ + + to % solution_file = from % solution_file + to % no_of_surfaceSamplings = from % no_of_surfaceSamplings + to % no_of_planeSamplings = from % no_of_planeSamplings + to % no_of_spatialMeanNodes = from % no_of_spatialMeanNodes + + to % dt_restriction = from % dt_restriction + to % write_dt_restriction = from % write_dt_restriction + + if ( to % no_of_surfaceSamplings .gt. 0) call to % surfaceSamplings % destruct + if ( to % no_of_planeSamplings .gt. 0) call to % planeSamplings % destruct + if ( to % no_of_spatialMeanNodes .gt. 0) call to % spatialMeanNodes % destruct + + safedeallocate ( to % surfaceSamplings ) + safedeallocate ( to % planeSamplings ) + safedeallocate ( to % spatialMeanNodes ) + + allocate ( to % surfaceSamplings ( size(from % surfaceSamplings) ) ) + allocate ( to % planeSamplings ( size(from % planeSamplings ) ) ) + allocate ( to % spatialMeanNodes ( size(from % spatialMeanNodes ) ) ) + to % surfaceSamplings = from % surfaceSamplings + to % planeSamplings = from % planeSamplings + to % spatialMeanNodes = from % spatialMeanNodes + + end subroutine Sampling_Assign + +! +!////////////////////////////////////////////////////////////////////////////// +! +! Auxiliars +! +!////////////////////////////////////////////////////////////////////////////// +! + subroutine getNoOfSamplings(no_of_surfaceSamplings, no_of_planeSamplings, no_of_spatialMeanNodes) + use ParamfileRegions + implicit none + integer, intent(out) :: no_of_surfaceSamplings + integer, intent(out) :: no_of_planeSamplings + integer, intent(out) :: no_of_spatialMeanNodes +! +! --------------- +! Local variables +! --------------- +! + character(len=LINE_LENGTH) :: case_name, line + integer :: fID + integer :: io +! +! Initialize +! ---------- + no_of_surfaceSamplings = 0 + no_of_planeSamplings = 0 + no_of_spatialMeanNodes = 0 +! +! Get case file name +! ------------------ + call get_command_argument(1, case_name) + +! +! Open case file +! -------------- + open ( newunit = fID , file = case_name , status = "old" , action = "read" ) + +! +! Read the whole file to find Samplings +! ------------------------------------ +readloop:do + read ( fID , '(A)' , iostat = io ) line + + if ( io .lt. 0 ) then +! +! End of file +! ----------- + line = "" + exit readloop + + elseif ( io .gt. 0 ) then +! +! Error +! ----- + errorMessage(STD_OUT) + stop "Stopped." + + else +! +! Succeeded +! --------- + line = getSquashedLine( line ) + + if ( index ( line , '#definesurfacesampling' ) .gt. 0 ) then + no_of_surfaceSamplings = no_of_surfaceSamplings + 1 + + else if ( index ( line , '#defineplanesampling' ) .gt. 0 ) then + no_of_planeSamplings = no_of_planeSamplings + 1 + + else if ( index ( line , '#definespatialmeannode' ) .gt. 0 ) then + no_of_spatialMeanNodes = no_of_spatialMeanNodes + 1 + + end if + + end if + + end do readloop +! +! Close case file +! --------------- + close(fID) + +end subroutine getNoOfSamplings + +end module Samplings +! +!/////////////////////////////////////////////////////////////////////////////////// +! \ No newline at end of file diff --git a/Solver/src/libs/monitors/SpatialMeanNode.f90 b/Solver/src/libs/monitors/SpatialMeanNode.f90 new file mode 100644 index 0000000000..30835380c7 --- /dev/null +++ b/Solver/src/libs/monitors/SpatialMeanNode.f90 @@ -0,0 +1,829 @@ +#include "Includes.h" +module SpatialMeanNode + use SMConstants + use HexMeshClass + use MonitorDefinitions + use PhysicsStorage + use VariableConversion + use MPI_Process_Info + use FluidData + use FileReadingUtilities, only: getRealArrayFromString, GetRealValue, getCharArrayFromString, GetIntValue, GetLogicalValue +#ifdef _HAS_MPI_ + use mpi +#endif + implicit none + + private + public SpatialMeanNode_t +! +! ********************************* +! SpatialMeanNodes class definition +! ********************************* +! + + type SpatialMeanNode_t + integer :: ID + integer :: nVariables + integer :: interval + integer :: bufferSize + integer :: bufferLine + integer :: intervalCount + integer :: nActive + integer :: dirAxis + integer :: nUniqueAll + integer :: iVarU, iVarV, iVarW + integer , allocatable :: activeLoc(:,:) + integer , allocatable :: nMultiply(:) + integer , allocatable :: nMultiplyAll(:) + logical :: meanData = .false. + real(kind=RP) :: pmin(3), pmax(3) + real(kind=RP) :: error=0.000001 ! tolerance of coordinate + real(kind=RP), allocatable :: geom(:) ! size nUnique + real(kind=RP), allocatable :: meanU(:), meanV(:), meanW(:) + real(kind=RP), allocatable :: values(:,:,:) ! (nUnique, bufferSize, nVariables) + character(len=STR_LEN_MONITORS), allocatable :: fileName (:) + character(len=STR_LEN_MONITORS) :: spatialMeanName + character(len=STR_LEN_MONITORS), allocatable :: variable (:) + contains + procedure :: Initialization => SpatialMeanNode_Initialization + procedure :: Update => SpatialMeanNode_Update + procedure :: WriteToFile => SpatialMeanNode_WriteToFile + procedure :: LookForUniqueCoordinate => SpatialMeanNode_LookForUniqueCoordinate + procedure :: destruct => SpatialMeanNode_Destruct + procedure :: copy => SpatialMeanNode_Assign + generic :: assignment(=) => copy + end type SpatialMeanNode_t + + contains + + subroutine SpatialMeanNode_Initialization(self, mesh, ID, solution_file, FirstCall) + use Headers + use ParamfileRegions + use MPI_Process_Info + use Utilities, only: toLower + implicit none + class(SpatialMeanNode_t) :: self + class(HexMesh) :: mesh + integer :: ID + character(len=*) :: solution_file + logical, intent(in) :: FirstCall +! +! --------------- +! Local variables +! --------------- +! + integer :: i, j, k, fID + real(kind=RP) :: point(2,3) + character(len=STR_LEN_MONITORS) :: in_label + character(len=STR_LEN_MONITORS) :: fileName + character(len=STR_LEN_MONITORS) :: paramFile + character(len=STR_LEN_MONITORS) :: interval, direction + character(len=STR_LEN_MONITORS) :: point1_char,point2_char,point3_char + character(len=STR_LEN_MONITORS) :: variables + character(len=STR_LEN_MONITORS) :: fileFormat + character(len=STR_LEN_MONITORS) :: writeInterval + character(len=STR_LEN_MONITORS) :: meanData +#if defined(NAVIERSTOKES) && (!(INCNS)) + if (FirstCall) then +! +! Get monitor ID, assign zero to bufferLine and intervalCount +! -------------- + self % ID = ID + self % bufferLine = 0 + self % intervalCount = 0 +! +! Search for the parameters in the case file +! ------------------------------------------ + write(in_label , '(A,I0)') "#define spatial mean node " , self % ID + + call get_command_argument(1, paramFile) + call readValueInRegion(trim(paramFile), "name" , self % spatialMeanName , in_label, "#end" ) + call readValueInRegion(trim(paramFile), "variables" , variables , in_label, "#end" ) + call readValueInRegion(trim(paramFile), "direction axis" , direction , in_label, "#end" ) + call readValueInRegion(trim(paramFile), "xrange [xmin,xmax]", point1_char , in_label, "#end" ) + call readValueInRegion(trim(paramFile), "yrange [ymin,ymax]", point2_char , in_label, "#end" ) + call readValueInRegion(trim(paramFile), "zrange [zmin,zmax]", point3_char , in_label, "#end" ) + call readValueInRegion(trim(paramFile), "sampling interval" , interval , in_label, "#end" ) + call readValueInRegion(trim(paramFile), "write interval" , writeInterval , in_label, "#end" ) + call readValueInRegion(trim(paramFile), "mean data" , meanData , in_label, "#end" ) +! +! Get the variables, points, N discretization, and interval +! ---------------------------------------------- + call getCharArrayFromString(variables, STR_LEN_MONITORS, self % variable) + self % nVariables = size(self % variable) + point(:,1) = getRealArrayFromString(point1_char) + point(:,2) = getRealArrayFromString(point2_char) + point(:,3) = getRealArrayFromString(point3_char) + + self % pmin(1) = point(1,1) + self % pmin(2) = point(1,2) + self % pmin(3) = point(1,3) + + self % pmax(1) = point(2,1) + self % pmax(2) = point(2,2) + self % pmax(3) = point(2,3) + + self % dirAxis = GetIntValue(direction) + self % interval = GetIntValue(interval) + if (len(trim(meanData)).lt.3) then + self % meanData = .false. + else + self % meanData = GetLogicalValue(meanData) + end if + + if ( .not. allocated(mesh % elements(1) % storage % stats % data ) ) then + self % meanData = .false. + end if + + interval=TRIM(ADJUSTL(interval)) +! +! Get the max. number of timestep in the buffer file before being written +! ----------------------------------------------------------------------- + IF (LEN(TRIM(writeInterval)) .EQ. 0) THEN + self % bufferSize = 1; + ELSE + self % bufferSize = GetIntValue(writeInterval) +! +! Failsafe to prevent too many data being written at one time +! ----------------------------------------------------------- + IF (self % bufferSize .GT. 10000) THEN + self % bufferSize = 10000; + END IF + END IF +! +! Look for unique data point in the range +! --------------------------------------- + call self % LookForUniqueCoordinate (mesh) +! +! Allocate Variables +! ------------------ + ALLOCATE(self % fileName (self % nVariables)) + self % values = 0_RP + end if + if (self % meanData) then + allocate(self % meanU(self % nUniqueAll), self % meanV(self % nUniqueAll), self % meanW(self % nUniqueAll)) + self % meanU=0_RP + self % meanV=0_RP + self % meanW=0_RP + end if +! +! Check Variables, Create Files, and Write Header Files +! ----------------------------------------------------- + do i=1,self % nVariables + if (self % meanData) then +! +! Prepare the file in which the SpatialMeanNode is exported +! --------------------------------------------------------- + write( self % fileName (i) , '(A,A,A,A,A,A,I0,A)') trim(solution_file) ,"_", trim(self % spatialMeanName) ,"_", trim(self % variable(i)) & + , "_spatialTemporalMean_" , self % ID,".node" +! +! Check the variable +! ------------------ + call tolower(self % variable (i)) +#ifdef NAVIERSTOKES + select case ( trim(self % variable (i)) ) + case ("density") + case ("pressure") + case ("ptotal") + case ("velocity") + case ("viscosity") + case ("u") + self % iVarU=i + case ("v") + self % iVarV=i + case ("w") + self % iVarW=i + case ("uu") + case ("vv") + case ("ww") + case ("uv") + case ("uw") + case ("vw") + case ("uprime2") + case ("vprime2") + case ("wprime2") + case ("uvprime") + case ("uwprime") + case ("vwprime") + case ("mach") + case ("k") + case ("q1") + case ("q2") + case ("q3") + case ("q4") + case ("q5") + case default + if ( MPI_Process % isRoot ) then + print*, 'SpatialMeanNode from temporal mean data, variable "',trim(self % variable(i)),'" not implemented.' + print*, "Options available are:" + print*, " * density" + print*, " * pressure" + print*, " * velocity" + print*, " * viscosity" + print*, " * u" + print*, " * v" + print*, " * w" + print*, " * uu" + print*, " * vv" + print*, " * ww" + print*, " * uv" + print*, " * uw" + print*, " * vw" + print*, " * uprime2" + print*, " * vprime2" + print*, " * wprime2" + print*, " * uvprime" + print*, " * uwprime" + print*, " * vwprime" + print*, " * Mach" + print*, " * K" + print*, " * q1" + print*, " * q2" + print*, " * q3" + print*, " * q4" + print*, " * q5" + end if + end select +#endif + else +! +! Prepare the file in which the SpatialMeanNode is exported +! --------------------------------------------------------- + write( self % fileName (i) , '(A,A,A,A,A,A,I0,A)') trim(solution_file) ,"_", trim(self % spatialMeanName) ,"_", trim(self % variable(i)) & + , "_spatialmean_" , self % ID,".node" +! +! Check the variable +! ------------------ + call tolower(self % variable (i)) + + select case ( trim(self % variable (i)) ) +#ifdef NAVIERSTOKES + case ("density") + case ("pressure") + case ("ptotal") + case ("velocity") + case ("viscosity") + case ("u") + case ("v") + case ("w") + case ("mach") + case ("k") + case ("omegax") + case ("omegay") + case ("omegaz") + case ("q1") + case ("q2") + case ("q3") + case ("q4") + case ("q5") + case default + if ( MPI_Process % isRoot ) then + print*, 'SpatialMeanNode variable "',trim(self % variable(i)),'" not implemented.' + print*, "Options available are:" + print*, " * density" + print*, " * pressure" + print*, " * velocity" + print*, " * viscosity" + print*, " * u" + print*, " * v" + print*, " * w" + print*, " * Mach" + print*, " * K" + print*, " * q1" + print*, " * q2" + print*, " * q3" + print*, " * q4" + print*, " * q5" + end if +#endif +#ifdef INCNS + case default + print*, "SpatialMeanNodes are not implemented for the incompressible NSE" +#endif +#ifdef MULTIPHASE + case default + print*, 'SpatialMeanNodes are not implemented.' +#endif + end select + end if + end do + + + if ( .not. MPI_Process % isRoot ) return + + do i=1,self % nVariables +#if defined(NAVIERSTOKES) && (!(INCNS)) +! +! Create file +! ----------- + if (FirstCall) then + if (MPI_Process % isRoot ) then + open ( newunit = fID , file = trim(self % fileName (i)) , action = "write" , access = "stream" , status = "replace", position='append' ) +! +! Write the file headers +! ---------------------- + write( fID) self % ID + write( fID) sum(self % nMultiplyAll) + write( fID) self % nUniqueAll + write( fID) self % dirAxis + write( fID ) self % interval + write( fID ) refValues % rho + write( fID ) refValues % V + write( fID ) refValues % p + write( fID ) refValues % T + write( fID ) refValues % mu + write( fID ) refValues % AoATheta + write( fID ) refValues % AoAPhi + write( fID ) self % geom + close ( fID ) + end if + end if +#endif + end do +#if defined(NAVIERSTOKES) && (!(INCNS)) +! +! File Format +! ----------------------------------------------- + write( fileFormat , '(A,A,A,A,I0,A)') trim(solution_file) ,"_", trim(self % spatialMeanName) ,"_'variable'_spatialmean_", self % ID, ".node" +! +! Write Information +! ----------------------------------------------- + write(STD_OUT,'(/)') + call SubSection_Header("SpatialMean Nodes") + write(STD_OUT,'(30X,A,A27,I4)') "->" , "SpatialMeanNode ID: " , self % ID + write(STD_OUT,'(30X,A,A27,A128)') "->" , "Variables: " , variables + write(STD_OUT,'(30X,A,A27,A,F6.4,A,F6.4,A)') "->" , "xrange [xmin,xmax] (m): ","[", & + self % pmin(1), ", ", self % pmax(1), "]" + write(STD_OUT,'(30X,A,A27,A,F6.4,A,F6.4,A)') "->" , "yrange [ymin,ymax] (m): ","[", & + self % pmin(2), ", ", self % pmax(2), "]" + write(STD_OUT,'(30X,A,A27,A,F6.4,A,F6.4,A)') "->" , "zrange [zmin,zmax] (m): ","[", & + self % pmin(3), ", ", self % pmax(3), "]" + write(STD_OUT,'(30X,A,A27,I5)') "->" , "Number of unique nodes: ", self % nUniqueAll + write(STD_OUT,'(30X,A,A27,I4)') "->" , "Samplings Interval: ", self % interval + write(STD_OUT,'(30X,A,A27,I4)') "->" , "Write Interval: ", self % bufferSize + write(STD_OUT,'(30X,A,A27,L1)') "->" , "Meanflow data: ", self % meanData + write(STD_OUT,'(30X,A,A27,A128)') "->" , "Filename: ", fileFormat +#endif + +#endif + end subroutine SpatialMeanNode_Initialization + + subroutine SpatialMeanNode_Update(self, mesh, bufferPosition, t) + use Physics + use MPI_Process_Info + use, intrinsic :: ieee_arithmetic, only: IEEE_Value, IEEE_QUIET_NAN + implicit none + class(SpatialMeanNode_t) :: self + class(HexMesh) :: mesh + integer :: bufferPosition + real(kind=RP) :: t +! +! --------------- +! Local variables +! --------------- +! + integer :: i, j, k, l, m, ierr, eID, pID + integer, parameter :: NO_OF_VARIABLES_Sij = 9 + integer, parameter :: U = 1 + integer, parameter :: V = 2 + integer, parameter :: W = 3 + integer, parameter :: UU = 4 + integer, parameter :: VV = 5 + integer, parameter :: WW = 6 + integer, parameter :: UV = 7 + integer, parameter :: UW = 8 + integer, parameter :: VW = 9 + real(kind=RP) :: value, rhoInv, kappa + real(kind=RP) , allocatable :: buff(:,:) +#ifdef NAVIERSTOKES +! +! Update data based on interval +! ----------------------------- + if (self % intervalCount .EQ. 0 ) then + self % bufferLine = self % bufferLine + 1 + DO m=1,self % nVariables + self % values(:, bufferPosition, m) = 0_RP + self % values(1, bufferPosition, m) = t + DO l=1, self % nActive +! +! Update the Node +! ---------------- + associate( e => mesh % elements(self % activeLoc(1,l)) ) + associate( Q => e % storage % Q, S => e % storage, meanQ => e % storage % stats % data) + eID=self % activeLoc(1,l) + i = self % activeLoc(2,l) + j = self % activeLoc(3,l) + k = self % activeLoc(4,l) + pID=self % activeLoc(5,l) + + select case (trim(self % variable(m))) + case("density") + self % values ( pID+1,bufferPosition, m) = self % values (pID+1,bufferPosition, m)+ Q(IRHO,i,j,k) + case("pressure") + self % values ( pID+1,bufferPosition, m) = self % values (pID+1,bufferPosition, m)+ Pressure(Q(:,i,j,k)) + case("ptotal") + value = POW2(Q(IRHOU,i,j,k)) + POW2(Q(IRHOV,i,j,k)) + POW2(Q(IRHOW,i,j,k))/POW2(Q(IRHO,i,j,k)) ! Vabs**2 + value = value / ( thermodynamics % gamma*(thermodynamics % gamma-1.0_RP)*(Q(IRHOE,i,j,k)/Q(IRHO,i,j,k)-0.5_RP * value) ) ! Mach ^2 + value = Pressure(Q(:,i,j,k))*(1.0_RP+0.5_RP*(thermodynamics % gamma-1.0_RP)*value)**( thermodynamics % gamma/(thermodynamics % gamma-1.0_RP)) + self % values ( pID+1,bufferPosition, m) = self % values (pID+1,bufferPosition, m)+ value + case("velocity") + value = sqrt(POW2(Q(IRHOU,i,j,k)) + POW2(Q(IRHOV,i,j,k)) + POW2(Q(IRHOW,i,j,k)))/Q(IRHO,i,j,k) + self % values ( pID+1,bufferPosition, m) = self % values (pID+1,bufferPosition, m)+ value + case("viscosity") + call get_laminar_mu_kappa(Q(:,i,j,k),value,kappa) + self % values ( pID+1,bufferPosition, m) = self % values (pID+1,bufferPosition, m)+ value + case("omegax") + value = (1/Q(IRHO,i,j,k) * S % U_y(IRHOW,i,j,k) - Q(IRHOW,i,j,k)/(Q(IRHO,i,j,k)**2) * S % U_y(IRHO,i,j,k)) & + - (1/Q(IRHO,i,j,k) * S % U_z(IRHOV,i,j,k) - Q(IRHOV,i,j,k)/(Q(IRHO,i,j,k)**2) * S % U_z(IRHO,i,j,k)) + self % values ( pID+1,bufferPosition, m) = self % values (pID+1,bufferPosition, m)+ value + case("omegay") + value = (1/Q(IRHO,i,j,k) * S % U_z(IRHOU,i,j,k) - Q(IRHOU,i,j,k)/(Q(IRHO,i,j,k)**2) * S % U_z(IRHO,i,j,k)) & + - (1/Q(IRHO,i,j,k) * S % U_x(IRHOW,i,j,k) - Q(IRHOW,i,j,k)/(Q(IRHO,i,j,k)**2) * S % U_x(IRHO,i,j,k)) + self % values ( pID+1,bufferPosition, m) = self % values (pID+1,bufferPosition, m)+ value + case("omegaz") + value = (1/Q(IRHO,i,j,k) * S % U_x(IRHOV,i,j,k) - Q(IRHOV,i,j,k)/(Q(IRHO,i,j,k)**2) * S % U_x(IRHO,i,j,k)) & + - (1/Q(IRHO,i,j,k) * S % U_y(IRHOU,i,j,k) - Q(IRHOU,i,j,k)/(Q(IRHO,i,j,k)**2) * S % U_y(IRHO,i,j,k)) + self % values ( pID+1,bufferPosition, m) = self % values (pID+1,bufferPosition, m)+ value + case("u") + self % values ( pID+1,bufferPosition, m) = self % values (pID+1,bufferPosition, m)+ Q(IRHOU,i,j,k) / Q(IRHO,i,j,k) + case("v") + self % values ( pID+1,bufferPosition, m) = self % values (pID+1,bufferPosition, m)+ Q(IRHOV,i,j,k) / Q(IRHO,i,j,k) + case("w") + self % values ( pID+1,bufferPosition, m) = self % values (pID+1,bufferPosition, m)+ Q(IRHOW,i,j,k) / Q(IRHO,i,j,k) + case("mach") + value = POW2(Q(IRHOU,i,j,k)) + POW2(Q(IRHOV,i,j,k)) + POW2(Q(IRHOW,i,j,k))/POW2(Q(IRHO,i,j,k)) ! Vabs**2 + value = sqrt( value / ( thermodynamics % gamma*(thermodynamics % gamma-1.0_RP)*(Q(IRHOE,i,j,k)/Q(IRHO,i,j,k)-0.5_RP * value) ) ) + self % values ( pID+1,bufferPosition, m) = self % values (pID+1,bufferPosition, m)+ value + case("k") + value = 0.5_RP * (POW2(Q(IRHOU,i,j,k)) + POW2(Q(IRHOV,i,j,k)) + POW2(Q(IRHOW,i,j,k)))/Q(IRHO,i,j,k) + self % values ( pID+1,bufferPosition, m) = self % values (pID+1,bufferPosition, m)+ value + case("q1") + self % values ( pID+1,bufferPosition, m) = self % values (pID+1,bufferPosition, m)+ Q(IRHO,i,j,k) + case("q2") + self % values ( pID+1,bufferPosition, m) = self % values (pID+1,bufferPosition, m)+ Q(IRHOU,i,j,k) + case("q3") + self % values ( pID+1,bufferPosition, m) = self % values (pID+1,bufferPosition, m)+ Q(IRHOV,i,j,k) + case("q4") + self % values ( pID+1,bufferPosition, m) = self % values (pID+1,bufferPosition, m)+ Q(IRHOW,i,j,k) + case("q5") + self % values ( pID+1,bufferPosition, m) = self % values (pID+1,bufferPosition, m)+ Q(IRHOE,i,j,k) + end select + end associate + end associate + END DO +! Combine all result from all proc MPI into root +! ---------------------------------------------- +#ifdef _HAS_MPI_ + call mpi_barrier(MPI_COMM_WORLD, ierr) + ALLOCATE(buff(self % nUniqueAll,MPI_Process % nProcs)) + call mpi_allgather(self % values(2:self % nUniqueAll+1, bufferPosition, m), self % nUniqueAll, MPI_DOUBLE, buff, self % nUniqueAll, MPI_DOUBLE, MPI_COMM_WORLD, ierr) + self % values(2:self % nUniqueAll+1, bufferPosition, m) = SUM(buff,DIM=2) + DEALLOCATE(buff) +#endif +! Average with the number of data +! ------------------------------- + do i = 1, self % nUniqueAll + self % values(i+1, bufferPosition, m) = self % values(i+1, bufferPosition, m) / self % nMultiplyAll(i) + end do + end do + if (self % meanData) then + self % meanU = self % values (2: self % nUniqueAll+1,bufferPosition,self % iVarU) + self % meanV = self % values (2: self % nUniqueAll+1,bufferPosition,self % iVarV) + self % meanW = self % values (2: self % nUniqueAll+1,bufferPosition,self % iVarW) + end if + end if + + self % intervalCount = self % intervalCount + 1 + + if (self % intervalCount .EQ. self % interval) then + self % intervalCount = 0 + end if +#endif + end subroutine SpatialMeanNode_Update + + subroutine SpatialMeanNode_WriteToFile ( self, no_of_lines) +! +! ************************************************************* +! This subroutine writes the buffer to the file. +! ************************************************************* +! + implicit none + class(SpatialMeanNode_t) :: self + integer :: no_of_lines +! +! --------------- +! Local variables +! --------------- +! + integer :: i, j, ierr + integer :: fID + + if ( MPI_Process % isRoot ) then + DO j=1, self % nVariables + open( newunit = fID , file = trim ( self % fileName (j) ) , action = "write" , access = "stream" , status = "old", position='append' ) + do i = 1 , no_of_lines + write( fID ) self % values(:,i,j) + end do + close ( fID ) + END DO + end if + if ( no_of_lines .ne. 0 ) self % values(:,1,:) = self % values(:,no_of_lines,:) + + self % bufferLine = 0 + + self % values = 0_RP + + end subroutine SpatialMeanNode_WriteToFile + + subroutine SpatialMeanNode_LookForUniqueCoordinate(self, mesh) + use MPI_Process_Info + implicit none + class(SpatialMeanNode_t) :: self + class(HexMesh) :: mesh +! +! --------------- +! Local variables +! --------------- +! + integer :: eID, i, j, k ,l , nPotentialUnique, nUnique, nCount + integer, allocatable :: uniqueProcs(:), activeLoc(:,:) + integer :: ierr + real(kind=RP), allocatable :: yUnique(:), yFinal(:), yUniqueProc(:,:) + real(kind=RP) :: error=0.000001 + real(kind=RP) :: buff, y + logical :: unique, finish +! +! Obtain the direction at which spatial mean is not applied and store the coordinate +! ---------------------------------------------------------------------------------- + nPotentialUnique=0 + do eID=1, mesh % no_of_elements + associate(e => mesh % elements(eID)) + nPotentialUnique=nPotentialUnique+(e % Nxyz(self % dirAxis) +1) * (e % Nxyz(1) +1) !nPotentialUnique Points + end associate + end do + ALLOCATE(activeLoc(1:5, mesh % NDOF), yUnique(nPotentialUnique)) + nUnique=0 + nCount =0 + do eID=1, mesh % no_of_elements + if (.not.((any(mesh % elements(eID) % geom % x(1,:,:,:).gt.(self % pmin(1)-error))) & + .and.(any(mesh % elements(eID) % geom % x(1,:,:,:).lt.(self % pmax(1)+error))))) then + cycle + end if + if (.not.((any(mesh % elements(eID) % geom % x(2,:,:,:).gt.(self % pmin(2)-error))) & + .and.(any(mesh % elements(eID) % geom % x(2,:,:,:).lt.(self % pmax(2)+error))))) then + cycle + end if + if (.not.((any(mesh % elements(eID) % geom % x(3,:,:,:).gt.(self % pmin(3)-error))) & + .and.(any(mesh % elements(eID) % geom % x(3,:,:,:).lt.(self % pmax(3)+error))))) then + cycle + end if + + associate(e => mesh % elements(eID)) + do i=0,e % Nxyz(3); do j=0,e % Nxyz(2); do k=0,e % Nxyz(1) + if (.not.(((e % geom % x(1,k,j,i)).gt.(self % pmin(1)-error)) & + .and.((e % geom % x(1,k,j,i)).lt.(self % pmax(1)+error)))) cycle + if (.not.(((e % geom % x(2,k,j,i)).gt.(self % pmin(2)-error)) & + .and.((e % geom % x(2,k,j,i)).lt.(self % pmax(2)+error)))) cycle + if (.not.(((e % geom % x(3,k,j,i)).gt.(self % pmin(3)-error)) & + .and.((e % geom % x(3,k,j,i)).lt.(self % pmax(3)+error)))) cycle +! +! Store the location of the active node ( within the bounded x,y,z ) +! ------------------------------------------------------------------ + nCount = nCount+1 + activeLoc(1, nCount) = eID + activeLoc(2, nCount) = k + activeLoc(3, nCount) = j + activeLoc(4, nCount) = i +! +! Find the unique points +! ---------------------- + y = e % geom % x( self % dirAxis, k,j,i) + unique=.true. + do l=1,nUnique + if (abs(yUnique(l)-y).lt.error) then + unique=.false. + exit + end if + end do + + if (unique) then + nUnique = nUnique + 1 + yUnique(nUnique)=y + end if + end do ; end do ; end do + end associate + end do + + ALLOCATE(yFinal(nUnique), self % activeLoc(5,nCount)) + self % nActive = nCount + yFinal=yUnique(1:nUnique) + self % activeLoc(1:4,:) = activeLoc(1:4,1:nCount) + DEALLOCATE(yUnique, activeLoc) +! +! MPI Operation to combine operation from different MPI +! ----------------------------------------------------- +#ifdef _HAS_MPI_ + ALLOCATE(uniqueProcs(MPI_Process % nProcs)) +! +! Gather all data from all processes +! ---------------------------------- + call mpi_allgather(nUnique, 1, MPI_INT, uniqueProcs, 1, MPI_INT, MPI_COMM_WORLD, ierr) + nPotentialUnique = sum(uniqueProcs) + call mpi_barrier(MPI_COMM_WORLD, ierr) +! +! Send unique coordinate on each proc to every proc with allgather +! ---------------------------------------------------------------- + ALLOCATE(yUnique(maxval(uniqueProcs))) + ALLOCATE(yUniqueProc(maxval(uniqueProcs),MPI_Process % nProcs)) + yUnique = 0_RP + yUnique(1:nUnique)=yFinal + call mpi_allgather(yUnique, maxval(uniqueProcs), MPI_DOUBLE, yUniqueProc, maxval(uniqueProcs), MPI_DOUBLE, MPI_COMM_WORLD, ierr) + call mpi_barrier(MPI_COMM_WORLD, ierr) + + nCount = uniqueProcs(1) + DEALLOCATE(yUnique) + ALLOCATE(yUnique(nPotentialUnique)) + yUnique(1:nCount) = yUniqueProc(1:nCount,1) + + do i=2, MPI_Process % nProcs + do j=1, uniqueProcs(i) + unique=.true. + if (any((yUnique-yUniqueProc(j,i)).lt.error)) then + unique = .false. + cycle + end if + if (unique)then + nCount = nCount+1 + yUnique(nCount)=yUniqueProc(j,i) + end if + end do + end do + + DEALLOCATE(yFinal, yUniqueProc) + ALLOCATE(yFinal(1:nCount)) + yFinal=yUnique(1:nCount) + DEALLOCATE (yUnique) +! +! Sort the coordinate in ascending order +! -------------------------------------- + CALL sortDoubleArrayMinMax(nCount, yFinal, 0) + self % nUniqueAll = nCount + + call mpi_barrier(MPI_COMM_WORLD, ierr) +#else + self % nUniqueAll = nUnique +#endif + ALLOCATE(self % values (self % nUniqueAll+1, self % bufferSize, self % nVariables), self % geom (self % nUniqueAll), & + self % nMultiply ( self % nUniqueAll), self % nMultiplyAll ( self % nUniqueAll)) + self % geom= yFinal +! +! Assign location of all active node w.r.t. unique coordinate location +! -------------------------------------------------------------------- + self % nMultiply = 0 + do i=1, self % nActive + eID = self % activeLoc(1,i) + do j= 1, self % nUniqueAll + buff = mesh % elements(eID) % geom % x (self % dirAxis, self % activeLoc(2,i), self % activeLoc(3,i), self % activeLoc(4,i)) + if (abs(buff - self % geom (j)).lt.error) then + self % nMultiply(j) = self % nMultiply(j) +1 ! Local for each proc + self % activeLoc(5,i) = j + exit + end if + end do + end do +#ifdef _HAS_MPI_ +! +! Gather nMultiply to be combined into nMultiplyAll in root +! --------------------------------------------------------- + self % nMultiplyAll = 0 + ALLOCATE(activeLoc(self % nUniqueAll, MPI_Process % nProcs)) + call mpi_allgather(self % nMultiply, self % nUniqueAll, MPI_INT, activeLoc, self % nUniqueAll, MPI_INT, MPI_COMM_WORLD, ierr) + + self % nMultiplyAll = sum(activeLoc, DIM=2) + + DEALLOCATE(activeLoc) + +#else + self % nMultiplyAll = self % nMultiply + +#endif + end subroutine SpatialMeanNode_LookForUniqueCoordinate +! +!//////////////////////////////////////////////////////////////////////// +! +! This Recursive Subroutine sort Double ArraySet at nColumn from its Minimum to Maximum Value +! When call for the first time k must be set to 0 -- WARNING DO NOT USE FOR HUGE ARRAY +! + RECURSIVE SUBROUTINE sortDoubleArrayMinMax(sizeArray, ArraySet, k) + IMPLICIT NONE + INTEGER ,INTENT(IN) :: sizeArray + REAL(kind=RP) ,DIMENSION(1:sizeArray) ,INTENT(INOUT) :: ArraySet + INTEGER ,INTENT(IN) :: k +! +! --------------- +! Local variables +! --------------- +! + REAL(kind=RP) :: BufferSet + INTEGER :: i +! +! --------------- +! Perform Sorting Value +! --------------- +! + SELECT CASE(k) + CASE(0) + DO i=1,sizeArray-1 + IF (ArraySet(i).GT.ArraySet(i+1)) THEN + BufferSet = ArraySet(i) + ArraySet(i) = ArraySet(i+1) + ArraySet(i+1) = BufferSet + IF(i.GT.1) THEN + CALL sortDoubleArrayMinMax(sizeArray, ArraySet, i) + END IF + END IF + END DO + CASE DEFAULT ! For Recursive + IF (ArraySet(k-1).GT.ArraySet(k)) THEN + BufferSet = ArraySet(k-1) + ArraySet(k-1) = ArraySet(k) + ArraySet(k) = BufferSet + IF(k.GT.2) THEN + CALL sortDoubleArrayMinMax(sizeArray, ArraySet, k-1) + END IF + END IF + END SELECT + + END SUBROUTINE sortDoubleArrayMinMax + + elemental subroutine SpatialMeanNode_Destruct (self) + implicit none + class(SpatialMeanNode_t), intent(inout) :: self + + safedeallocate (self % activeLoc) + safedeallocate (self % nMultiply) + safedeallocate (self % nMultiplyAll) + safedeallocate (self % geom) + safedeallocate (self % values) + safedeallocate (self % meanU) + safedeallocate (self % meanV) + safedeallocate (self % meanW) + safedeallocate (self % fileName) + safedeallocate (self % variable) + + end subroutine SpatialMeanNode_Destruct + + elemental subroutine SpatialMeanNode_Assign (to, from) + implicit none + class(SpatialMeanNode_t), intent(inout) :: to + type(SpatialMeanNode_t) , intent(in) :: from + + + to % ID = from % ID + to % nVariables = from % nVariables + to % interval = from % interval + to % bufferSize = from % bufferSize + to % bufferLine = from % bufferLine + to % intervalCount = from % intervalCount + to % nActive = from % nActive + to % dirAxis = from % dirAxis + to % nUniqueAll = from % nUniqueAll + to % pmin = from % pmin + to % pmax = from % pmax + to % error = from % error + + safedeallocate ( to % activeLoc ) + allocate ( to % activeLoc( size(from % activeLoc,1),size(from % activeLoc,2) ) ) + to % activeLoc = from % activeLoc + + safedeallocate ( to % nMultiply ) + allocate ( to % nMultiply( size(from % nMultiply) ) ) + to % nMultiply = from % nMultiply + + safedeallocate ( to % nMultiplyAll ) + allocate ( to % nMultiplyAll( size(from % nMultiplyAll) ) ) + to % nMultiplyAll = from % nMultiplyAll + + safedeallocate ( to % geom ) + allocate ( to % geom( size(from % geom) ) ) + to % geom = from % geom + + safedeallocate ( to % values ) + allocate ( to % values( size(from % values,1),size(from % values,2),size(from % values,3) ) ) + to % values = from % values + + safedeallocate ( to % fileName ) + allocate ( to % fileName ( size(from % fileName) ) ) + to % fileName = from % fileName + + to % spatialMeanName = from % spatialMeanName + + safedeallocate ( to % variable ) + allocate ( to % variable ( size(from % variable) ) ) + to % variable = from % variable + + + + end subroutine SpatialMeanNode_Assign + +end module SpatialMeanNode \ No newline at end of file diff --git a/Solver/src/libs/monitors/SurfaceSampling.f90 b/Solver/src/libs/monitors/SurfaceSampling.f90 new file mode 100644 index 0000000000..8a916b83b1 --- /dev/null +++ b/Solver/src/libs/monitors/SurfaceSampling.f90 @@ -0,0 +1,558 @@ +#include "Includes.h" +module SurfaceSampling + use SMConstants + use HexMeshClass + use MonitorDefinitions + use PhysicsStorage + use MPI_Process_Info + use FluidData + use FileReadingUtilities, only: getRealArrayFromString, GetRealValue, getCharArrayFromString +#ifdef _HAS_MPI_ + use mpi +#endif + implicit none + + private + public SurfaceSampling_t + + +! +! ******************************** +! Surface Sampling class definition +! ******************************** +! + type SurfaceSampling_t + logical :: active=.false. + logical :: isDimensionless + integer :: ID + integer :: nVariables + integer :: marker + integer :: rank + integer :: interval + integer :: bufferSize + integer :: bufferLine + integer :: intervalCount + integer, allocatable :: nData(:) + real(kind=RP), allocatable :: values(:,:,:) + character(len=STR_LEN_MONITORS) :: SamplingName + character(len=STR_LEN_MONITORS), allocatable :: fileName(:) + character(len=STR_LEN_MONITORS), allocatable :: variable(:) + contains + procedure :: Initialization => SurfaceSampling_Initialization + procedure :: Update => SurfaceSampling_Update + procedure :: WriteToFile => SurfaceSampling_WriteToFile + procedure :: destruct => SurfaceSampling_Destruct + procedure :: copy => SurfaceSampling_Assign + generic :: assignment(=) => copy + end type SurfaceSampling_t + + contains +! +!///////////////////////////////////////////////////////////////////////// +! +! SURFACE Sampling PROCEDURES +! -------------------------- +!///////////////////////////////////////////////////////////////////////// +! + subroutine SurfaceSampling_Initialization( self , mesh , ID, solution_file , FirstCall) +! +! ***************************************************************************** +! This subroutine initializes the surface Sampling. The following +! data is obtained from the case file: +! -> Name: The Sampling name (10 characters maximum) +! -> Marker: The surface marker in which the Sampling will be computed. +! -> Variable: The variable to be Samplingized. +! ***************************************************************************** +! + use Headers + use ParamfileRegions + use MPI_Process_Info + implicit none + class(SurfaceSampling_t):: self + class(HexMesh) :: mesh + integer :: ID + character(len=*) :: solution_file + logical, intent(in) :: FirstCall +! +! --------------- +! Local variables +! --------------- +! + character(len=STR_LEN_MONITORS) :: in_label + character(len=STR_LEN_MONITORS) :: fileName + character(len=STR_LEN_MONITORS) :: paramFile + integer, allocatable :: marker + character(len=STR_LEN_MONITORS) :: markerName + character(len=STR_LEN_MONITORS) :: bufferInt + character(len=STR_LEN_MONITORS) :: formatFile + character(len=STR_LEN_MONITORS) :: variables + character(len=STR_LEN_MONITORS) :: writeInterval + integer :: pos + integer :: fID + integer :: zoneID + integer :: Nf(2) + integer :: nData + integer :: allnData(MPI_Process % nProcs), ierr + integer :: zonefID + integer :: i,j,k + real(kind=RP), allocatable :: X(:,:) + real(kind=RP), allocatable :: X1(:,:) + real(kind=RP), allocatable :: X2(:,:) + real(kind=RP), allocatable :: X3(:,:) + +! +! Get Sampling ID, assign zero to bufferLine and intervalCount +! -------------- + self % ID = ID + self % bufferLine = 0 + self % intervalCount = 0 +! +! Search for the parameters in the case file +! ------------------------------------------ + write(in_label , '(A,I0)') "#define surface sampling " , self % ID + + call get_command_argument(1, paramFile) + call readValueInRegion ( trim ( paramFile ) , "name" , self % SamplingName , in_label , "# end" ) + call readValueInRegion ( trim ( paramFile ) , "surface" , markerName , in_label , "# end" ) + call readValueInRegion ( trim ( paramFile ) , "variables" , variables , in_label , "# end" ) + call readValueInRegion ( trim ( paramFile ) , "sampling interval" , bufferInt , in_label , "# end" ) + call readValueInRegion(trim(paramFile), "write interval", writeInterval , in_label, "#end" ) +! +! Enable the Sampling +! ------------------ + self % active = .true. +! +! Get the variables and its size +! ------------------------------ + call getCharArrayFromString(variables, STR_LEN_MONITORS, self % variable) + self % nVariables = size(self % variable) +! +! Get the surface marker +! ---------------------- + self % marker = -1 + do zoneID = 1, size(mesh % zones) + if ( trim(mesh % zones(zoneID) % name) .eq. trim(markerName) ) then + self % marker = zoneID + exit + end if + end do + + if ( self % marker .eq. -1 ) then + self % active = .false. + if (MPI_Process % isRoot ) then + write(*,'(A,I0)') "Warning: Marker not specified for surface sampling ", self % ID + write(*,'(A,I0,A)') " Surface Sampling ", self % ID, " disabled." + end if + end if + + if (mesh % zones(self % marker) % no_of_faces .eq. 0) then + self % active = .false. + end if +! +! Select the variable from the available list, and compute auxiliary variables if needed +! -------------------------------------------------------------------------------------- +! + DO i=1, self % nVariables +! **************************************** + select case ( trim ( self % variable(i) ) ) +! **************************************** +! + case ("shearstress-tangent") + self % isDimensionless = .false. + + case ("shearstress-x") + self % isDimensionless = .false. + + case ("shearstress-y") + self % isDimensionless = .false. + + case ("shearstress-z") + self % isDimensionless = .false. + + case ("pressure") + self % isDimensionless = .false. + + case ("q1") + self % isDimensionless = .true. + + case ("q2") + self % isDimensionless = .true. + + case ("q3") + self % isDimensionless = .true. + + case ("q4") + self % isDimensionless = .true. + + case ("q5") + self % isDimensionless = .true. + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + case default + if (MPI_Process % isRoot ) then + if ( len_trim (self % variable(i)) .eq. 0 ) then + print*, "Variable was not specified for surface Sampling " , self % ID , "." + else + print*, 'Variable "',trim(self % variable(i)),'" surface Sampling ', self % ID, ' not implemented yet.' + print*, "Options available are:" + print*, " * shearstress-tangent" + print*, " * shearstress-x" + print*, " * shearstress-y" + print*, " * shearstress-z" + print*, " * pressure" + stop "Stopped." + + end if + end if +! +! ********** + end select +! ********** + END DO +! +! Get the number of Order and the sampling interval +! ------------------------------------------------- + if (mesh % zones(self % marker) % no_of_faces .gt.0) then + Nf(1) = mesh % faces (mesh % zones(self % marker) % faces(1))%Nf(1)+1 + Nf(2) = mesh % faces (mesh % zones(self % marker) % faces(1))%Nf(2)+1 + nData = mesh % zones(self % marker) % no_of_faces * Nf(1) * Nf(2) + else + self % active = .false. + Nf(1)=0 + Nf(2)=0 + nData=0 + end if + self % interval = GetRealValue(bufferInt) + + ALLOCATE (self % nData(MPI_Process % nProcs)) + + if ( MPI_Process % doMPIAction ) then +#ifdef _HAS_MPI_ + call mpi_allgather(nData, 1, MPI_INT, allnData, 1, MPI_INT, MPI_COMM_WORLD, ierr) + nData=sum(allnData) + self % nData = allnData + self % rank = maxloc(allnData,1)-1 + call mpi_allgather(Nf(1), 1, MPI_INT, allnData, 1, MPI_INT, MPI_COMM_WORLD, ierr) + Nf(1)=maxval(allnData) + call mpi_allgather(Nf(2), 1, MPI_INT, allnData, 1, MPI_INT, MPI_COMM_WORLD, ierr) + Nf(2)=maxval(allnData) +#endif + end if + +! +! Get the max. number of timestep in the buffer file before being written +! ----------------------------------------------------------------------- + IF (LEN(TRIM(writeInterval)) .EQ. 0) THEN + self % bufferSize = 1; + ELSE + self % bufferSize = GetRealValue(writeInterval) +! +! Failsafe to prevent too many data being written at one time +! ----------------------------------------------------------- + IF (nData .GT. 200000) THEN + self % bufferSize = 1; + END IF + END IF + + + ALLOCATE( X(nDATA,3), self % fileName(self % nVariables), X1(nData, MPI_Process % nProcs), & + X2(nData, MPI_Process % nProcs), X3(nData, MPI_Process % nProcs) ) + + if ( MPI_Process % isRoot ) then + ALLOCATE(self % values(nData+1,self % bufferSize , self % nVariables)) + end if + + k=0 + do zonefID = 1, mesh % zones(self % marker) % no_of_faces +! +! Face global ID +! -------------- + fID = mesh % zones(self % marker) % faces(zonefID) + do j = 0, Nf(2)-1; do i = 0, Nf(1)-1 + k=k+1 + X(k,1)=mesh % faces(fID) % geom % x(1,i,j) + X(k,2)=mesh % faces(fID) % geom % x(2,i,j) + X(k,3)=mesh % faces(fID) % geom % x(3,i,j) + end do; end do + + end do + + if ( MPI_Process % doMPIAction ) then +#ifdef _HAS_MPI_ + call mpi_allgather(X(:,1), nData, MPI_DOUBLE, X1, nData, MPI_DOUBLE, MPI_COMM_WORLD, ierr) + call mpi_allgather(X(:,2), nData, MPI_DOUBLE, X2, nData, MPI_DOUBLE, MPI_COMM_WORLD, ierr) + call mpi_allgather(X(:,3), nData, MPI_DOUBLE, X3, nData, MPI_DOUBLE, MPI_COMM_WORLD, ierr) + if ( MPI_Process % isRoot ) then + X(:,1)=X1(:,1) + X(:,2)=X2(:,1) + X(:,3)=X3(:,1) + do i=1, MPI_Process % nProcs-1 + X(sum(self % nData(1:i))+1:sum(self % nData(1:i+1)),1)=X1(1:self % nData(i+1),i+1) + X(sum(self % nData(1:i))+1:sum(self % nData(1:i+1)),2)=X2(1:self % nData(i+1),i+1) + X(sum(self % nData(1:i))+1:sum(self % nData(1:i+1)),3)=X3(1:self % nData(i+1),i+1) + end do + end if +#endif + end if + + + DO i=1, self % nVariables +! +! Prepare the file in which the Sampling is exported +! ------------------------------------------------- + write( self % fileName(i) , '(A,A,A,A,A,A,I0,A)') trim(solution_file) , "_" , trim(markerName) , "_" , trim(self % variable(i)) & + , "_surface_" , self % ID,".sampling" + +#if defined(NAVIERSTOKES) && (!(INCNS)) +! +! Create file +! ----------- + if (FirstCall) then + if (MPI_Process % isRoot ) then + open ( newunit = fID , file = trim(self % fileName(i)) , action = "write" , access = "stream" , status = "replace", position='append' ) + +! +! Write the file headers +! ---------------------- + write( fID) self % ID + write( fID) Nf(1) + write( fID) Nf(2) + write( fID) nData + write( fID ) self % interval + write( fID ) refValues % rho + write( fID ) refValues % V + write( fID ) refValues % p + write( fID ) refValues % T + write( fID ) refValues % mu + write( fID ) refValues % AoATheta + write( fID ) refValues % AoAPhi + write( fID ) X(:,1) + write( fID ) X(:,2) + write( fID ) X(:,3) + close ( fID ) + end if + end if +#endif + END DO + +! +! Write Information into log +! -------------------------- + write( formatFile , '(A,A,A,A,I0,A)') trim(solution_file) , "_" , trim(markerName) , "_'variable'_surface_", self % ID, ".sampling" + + if ( .not. MPI_Process % isRoot ) return +! +! Write Information +! ----------------------------------------------- + + write(STD_OUT,'(/)') + call SubSection_Header("Surface Samplings") + write(STD_OUT,'(30X,A,A27,I4)') "->" , "Surface Sampling ID: " , self % ID + write(STD_OUT,'(30X,A,A27,A27)') "->" , "Surface Name: " , markerName + write(STD_OUT,'(30X,A,A27,A128)') "->" , "Variables: " , variables + write(STD_OUT,'(30X,A,A27,I4)') "->" , "Samplings Interval: ", self % interval + write(STD_OUT,'(30X,A,A27,A128)') "->" , "Filename: ", formatFile + write(STD_OUT,'(30X,A,A27,A25)') "->" ,"Note: ","Extracted with dimension" + + deallocate(X,X1,X2,X3) + + end subroutine SurfaceSampling_Initialization + + subroutine SurfaceSampling_Update ( self, mesh, bufferPosition, t ) +! +! ******************************************************************* +! This subroutine updates the Sampling value computing it from +! the mesh. It is stored in the "bufferPosition" position of the +! buffer. +! ******************************************************************* +! + use SamplingOperator + implicit none + class ( SurfaceSampling_t ) :: self + class ( HexMesh ) :: mesh + integer :: bufferPosition + real(kind=RP) :: t + integer :: i, j, k, recv_req(MPI_Process % nProcs-1), send_req, ierr + real(kind=RP) :: F(NDIM) + real(kind=RP), allocatable :: data_out(:) + + + if (self % intervalCount .EQ. 0 ) then + self % bufferLine = self % bufferLine + 1 + + if ((self % nData(MPI_Process % rank +1) .gt.0 ).or. MPI_Process % isRoot) then + DO i=1, self % nVariables + select case ( trim ( self % variable(i) ) ) + + case ("shearstress-tangent") + Call VectorSurfaceSampling(mesh, self % marker, SHEAR_STRESS_TANGENT, self % SamplingName, data_out) + + case ("shearstress-x") + Call VectorSurfaceSampling(mesh, self % marker, SHEAR_STRESS_X, self % SamplingName, data_out) + + case ("shearstress-y") + Call VectorSurfaceSampling(mesh, self % marker, SHEAR_STRESS_Y, self % SamplingName, data_out) + + case ("shearstress-z") + Call VectorSurfaceSampling(mesh, self % marker, SHEAR_STRESS_Z, self % SamplingName, data_out) + + case ("pressure") + Call VectorSurfaceSampling(mesh, self % marker, PRESSURE_SURF, self % SamplingName, data_out) + + case ("q1") + Call VectorSurfaceSampling(mesh, self % marker, Q1, self % SamplingName, data_out) + + case ("q2") + Call VectorSurfaceSampling(mesh, self % marker, Q2, self % SamplingName, data_out) + + case ("q3") + Call VectorSurfaceSampling(mesh, self % marker, Q3, self % SamplingName, data_out) + + case ("q4") + Call VectorSurfaceSampling(mesh, self % marker, Q4, self % SamplingName, data_out) + + case ("q5") + Call VectorSurfaceSampling(mesh, self % marker, Q5, self % SamplingName, data_out) + + end select + + if ( MPI_Process % doMPIAction ) then +#ifdef _HAS_MPI_ + if ( MPI_Process % isRoot ) then + + self % values(:,bufferPosition,i)=0_RP + self % values(1,bufferPosition,i)=t + if (self % nData(1) .gt.0 ) then + self % values(2:self % nData(1) +1,bufferPosition,i)=data_out + end if + k=0 + DO j=1, MPI_Process % nProcs -1 + if (self % nData(j+1) .gt.0 ) then + k=k+1 + call mpi_irecv(self % values(sum(self % nData(1:j))+1:sum(self % nData(1:j))+1+self % nData(j+1),bufferPosition,i), & + self % nData(j+1), MPI_DOUBLE, j, MPI_ANY_TAG, MPI_COMM_WORLD, recv_req(k), ierr) + call mpi_wait(recv_req(k), MPI_STATUS_IGNORE, ierr) + end if + END DO + else + + if (self % nData(MPI_Process % rank +1) .gt.0 ) then + call mpi_isend(data_out, self % nData(MPI_Process % rank +1), MPI_DOUBLE, 0, & + DEFAULT_TAG, MPI_COMM_WORLD, send_req, ierr) + call mpi_wait(send_req, MPI_STATUS_IGNORE, ierr) + send_req=0 + end if + end if +#endif + else + self % values(1,bufferPosition,i)=t + self % values(2:size(data_out,1)+1,bufferPosition,i)=data_out + end if + + END DO + + end if + + + end if + + self % intervalCount = self % intervalCount + 1 + + if (self % intervalCount .EQ. self % interval) then + self % intervalCount = 0 + end if + + end subroutine SurfaceSampling_Update + + subroutine SurfaceSampling_WriteToFile ( self , no_of_lines) +! +! ************************************************************* +! This subroutine writes the buffer to the file. +! ************************************************************* +! + implicit none + class(SurfaceSampling_t) :: self + integer :: no_of_lines +! +! --------------- +! Local variables +! --------------- +! + integer :: i,k + integer :: fID + + if ( MPI_Process % isRoot ) then + + DO k=1, self % nVariables + + open( newunit = fID , file = trim ( self % fileName(k) ) , action = "write" , access = "stream" , status = "old", position='append' ) + + do i = 1 , no_of_lines + write( fID ) self % values(:,i,k) + + end do + + close ( fID ) + + if ( no_of_lines .ne. 0 ) self % values(:,1,k) = self % values(:,no_of_lines,k) + + END DO + + end if + + self % bufferLine = 0 + + end subroutine SurfaceSampling_WriteToFile +! + elemental subroutine SurfaceSampling_Destruct (self) + implicit none + class(SurfaceSampling_t), intent(inout) :: self + + safedeallocate (self % values) + safedeallocate (self % nData) + safedeallocate (self % fileName) + safedeallocate (self % variable) + + end subroutine SurfaceSampling_Destruct +! +!////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// +! + elemental subroutine SurfaceSampling_Assign(to, from) + implicit none + class(SurfaceSampling_t), intent(inout) :: to + type(SurfaceSampling_t), intent(in) :: from + + if ( from % active) then + to % active = from % active + to % isDimensionless = from % isDimensionless + to % ID = from % ID + to % nVariables = from % nVariables + to % marker = from % marker + to % interval = from % interval + to % bufferSize = from % bufferSize + to % bufferLine = from % bufferLine + to % intervalCount = from % intervalCount + + safedeallocate(to % values) + allocate ( to % values( size(from % values,1),size(from % values,2),size(from % values,3) ) ) + to % values = from % values + + to % SamplingName = from % SamplingName + + safedeallocate(to % fileName) + allocate ( to % fileName( size(from % fileName) ) ) + to % fileName = from % fileName + + safedeallocate(to % variable) + allocate ( to % variable( size(from % variable) ) ) + to % variable = from % variable + + else + to % active = .FALSE. + end if + + end subroutine SurfaceSampling_Assign +end module SurfaceSampling + diff --git a/Solver/src/libs/sources/SpongeClass.f90 b/Solver/src/libs/sources/SpongeClass.f90 index 9c0076a072..24908bb504 100644 --- a/Solver/src/libs/sources/SpongeClass.f90 +++ b/Solver/src/libs/sources/SpongeClass.f90 @@ -28,6 +28,10 @@ Module SpongeClass ! real(kind=RP) :: amplitude ! amplitude of the source term real(kind=RP) :: delta ! temporal filter width real(kind=RP) :: rampWidth ! length of the ramp width + real(kind=RP), dimension(:), allocatable :: fringeStart ! start of Fringe Region + real(kind=RP), dimension(:), allocatable :: fringeEnd ! end of Fringe Region + real(kind=RP), dimension(:), allocatable :: deltaRise ! rampRise of Fringe Region + real(kind=RP), dimension(:), allocatable :: deltaFall ! rampFall of Fringe Region real(kind=RP), dimension(:,:),allocatable :: x0 ! position of start of the sponge (for cylindrical in the center) real(kind=RP), dimension(:), allocatable :: radious ! radious of the ramp zone in cylindrical/cirular sponges real(kind=RP), dimension(:,:) ,allocatable :: axis ! axis vector of the sponge. In cylindrical axis of the cylinder, in cartesian, the aligned vector @@ -57,9 +61,10 @@ Module SpongeClass ! integer, parameter :: KEYWORD_LENGTH = 132 character(len=KEYWORD_LENGTH), parameter :: SPONGE_CYLINDRICAL_NAME = "cylindrical" character(len=KEYWORD_LENGTH), parameter :: SPONGE_CARTESIAN_NAME = "cartesian" + character(len=KEYWORD_LENGTH), parameter :: SPONGE_FRINGE_NAME = "fringe" enum, bind(C) - enumerator :: SPONGE_CYLINDRICAL = 1, SPONGE_CARTESIAN + enumerator :: SPONGE_CYLINDRICAL = 1, SPONGE_CARTESIAN, FRINGE end enum contains @@ -100,6 +105,7 @@ Subroutine spongeConstruct(self, mesh, controlVariables) allocate(self % radious(self % numOfSponges)) allocate(self % axis(self % numOfSponges,NDIM)) allocate(self % x0(self % numOfSponges,NDIM)) + allocate(self % fringeStart(self % numOfSponges),self % fringeEnd(self % numOfSponges),self % deltaRise(self % numOfSponges),self % deltaFall(self % numOfSponges)) do i = 1, self% numOfSponges write(tmp, '("sponge shape ",I0)') i @@ -119,23 +125,52 @@ Subroutine spongeConstruct(self, mesh, controlVariables) end if self % radious(i) = controlVariables % getValueOrDefault(tmp,0.0_RP) endif - if (useNumberedKeys) then - write(tmp, '("sponge axis ",I0)') i - else - write(tmp, '("sponge axis")') - end if - axis = controlVariables % stringValueForKey(tmp, requestedLength = STRING_CONSTANT_LENGTH) - self % axis(i,:) = getRealArrayFromString(trim(axis)) - !normalize axis - self % axis(i,:) = self % axis(i,:)/sqrt(sum(self % axis(i,:)**2)) - - if (useNumberedKeys) then - write(tmp, '("sponge start position ",I0)') i - else - write(tmp, '("sponge start position")') - end if - coordinates = controlVariables % stringValueForKey(tmp, requestedLength = STRING_CONSTANT_LENGTH) - self % x0(i,:) = getRealArrayFromString(trim(coordinates)) + if((self % shapeType(i) == "cylindrical").OR.(self % shapeType(i) == "cartesian")) then + if (useNumberedKeys) then + write(tmp, '("sponge axis ",I0)') i + else + write(tmp, '("sponge axis")') + end if + axis = controlVariables % stringValueForKey(tmp, requestedLength = STRING_CONSTANT_LENGTH) + self % axis(i,:) = getRealArrayFromString(trim(axis)) + !normalize axis + self % axis(i,:) = self % axis(i,:)/sqrt(sum(self % axis(i,:)**2)) + + if (useNumberedKeys) then + write(tmp, '("sponge start position ",I0)') i + else + write(tmp, '("sponge start position")') + end if + coordinates = controlVariables % stringValueForKey(tmp, requestedLength = STRING_CONSTANT_LENGTH) + self % x0(i,:) = getRealArrayFromString(trim(coordinates)) + end if + + if(self % shapeType(i) == "fringe") then + if (useNumberedKeys) then + write(tmp, '("fringe start ",I0)') i + else + write(tmp, '("fringe start")') + end if + self % fringeStart(i) = controlVariables % getValueOrDefault(tmp,0.0_RP) + if (useNumberedKeys) then + write(tmp, '("fringe end ",I0)') i + else + write(tmp, '("fringe end")') + end if + self % fringeEnd(i) = controlVariables % getValueOrDefault(tmp,0.0_RP) + if (useNumberedKeys) then + write(tmp, '("delta rise ",I0)') i + else + write(tmp, '("delta rise")') + end if + self % deltaRise(i) = controlVariables % getValueOrDefault(tmp,0.0_RP) + if (useNumberedKeys) then + write(tmp, '("delta fall ",I0)') i + else + write(tmp, '("delta fall")') + end if + self % deltaFall(i) = controlVariables % getValueOrDefault(tmp,0.0_RP) + endif end do ! Get the file name of the solution @@ -180,8 +215,16 @@ Subroutine spongeConstruct(self, mesh, controlVariables) if(self % shapeType(i) == "cylindrical") then write(STD_OUT,'(30X,A,A28,F10.2)') "->", "Ramp radious: ", self % radious(i) endif - write(STD_OUT,'(30X,A,A28,3(F10.2))') "->", "Axis: ", self % axis(i,:) - write(STD_OUT,'(30X,A,A28,3(F10.2))') "->", "Ramp start: ", self % x0(i,:) + if((self % shapeType(i) == "cylindrical").OR.(self % shapeType(i) == "cartesian")) then + write(STD_OUT,'(30X,A,A28,3(F10.2))') "->", "Axis: ", self % axis(i,:) + write(STD_OUT,'(30X,A,A28,3(F10.2))') "->", "Ramp start: ", self % x0(i,:) + end if + if(self % shapeType(i) == "fringe") then + write(STD_OUT,'(30X,A,A28,3(F10.2))') "->", "Fringe Start: ", self % fringeStart(i) + write(STD_OUT,'(30X,A,A28,3(F10.2))') "->", "Fringe End : ", self % fringeEnd(i) + write(STD_OUT,'(30X,A,A28,3(F10.2))') "->", "delta Rise : ", self % deltaRise(i) + write(STD_OUT,'(30X,A,A28,3(F10.2))') "->", "delta Fall : ", self % deltaFall(i) + end if end do if (self % readBaseFLowFlag) write(STD_OUT,'(30X,A,A28,A)') "->", "Initial base flow file: ", self % initialFileName @@ -200,6 +243,10 @@ Subroutine spongeDestruct(self) safedeallocate(self % radious) safedeallocate(self % axis) safedeallocate(self % shapeType) + safedeallocate(self % fringeStart) + safedeallocate(self % fringeEnd) + safedeallocate(self % deltaRise) + safedeallocate(self % deltaFall) End Subroutine spongeDestruct ! @@ -215,6 +262,7 @@ Subroutine creatRamp(self, mesh) !local variables real(kind=RP), dimension(:,:,:), allocatable :: xStar, sigma real(kind=RP), dimension(NDIM) :: r_vector + real(kind=RP) :: xData(2), Sf(2) logical, dimension(:), allocatable :: hasSponge integer :: i, j, k, eID, counter, spongeEID, ierr integer, dimension(NDIM) :: Nxyz @@ -237,6 +285,8 @@ Subroutine creatRamp(self, mesh) whichSponge = SPONGE_CYLINDRICAL case (SPONGE_CARTESIAN_NAME) whichSponge = SPONGE_CARTESIAN + case (SPONGE_FRINGE_NAME) + whichSponge = FRINGE case default print*, "Sponge name not recognized." errorMessage(STD_OUT) @@ -259,6 +309,11 @@ Subroutine creatRamp(self, mesh) case (SPONGE_CARTESIAN) ! in this case xstar is the distance to the plane xStar(i,j,k) = sum(r_vector(:)*self % axis(sponge_number,:)) + case (FRINGE) + ! in this case xstar is the distance to the plane + if ((e % geom % x(1,i,j,k).ge.self % fringeStart(sponge_number)).and.(e % geom % x(1,i,j,k).le.self % fringeEnd(sponge_number))) then + xStar(i,j,k) = 1.0_RP + end if end select end do ; end do ; end do if (any(xStar .ge. 0.0_RP) .AND. .not.hasSponge(eID) ) then @@ -298,6 +353,8 @@ Subroutine creatRamp(self, mesh) whichSponge = SPONGE_CYLINDRICAL case (SPONGE_CARTESIAN_NAME) whichSponge = SPONGE_CARTESIAN + case (SPONGE_FRINGE_NAME) + whichSponge = FRINGE end select do spongeEID = 1, self % nElements @@ -319,18 +376,50 @@ Subroutine creatRamp(self, mesh) xStar(i,j,k) = sqrt(abs(sum(r_vector*r_vector) - POW2(self % radious(sponge_number)))) / self % rampWidth case (SPONGE_CARTESIAN) xStar(i,j,k) = sum(r_vector(:)*self % axis(sponge_number,:))/(self % rampWidth) + case (FRINGE) + xData(1)=(e % geom % x(1,i,j,k)-self % fringeStart(sponge_number))/self % deltaRise(sponge_number) + xData(2)=((e % geom % x(1,i,j,k)-self % fringeEnd(sponge_number))/self % deltaFall(sponge_number))+1 + + if (xData(1).LE.0.0_RP) then + Sf(1)=0.0_RP + else if (xData(1).GE.1.0_RP) then + Sf(1)=1.0_RP + else + !if (abs(xData(1) - 1.0_RP) > epsilon(xData(1)) .and. abs(xData(1)) > epsilon(xData(1))) then + Sf(1)=1.0_RP/(1.0_RP+exp((1.0_RP/max((xData(1)-1.0_RP),0.01_RP))+(1.0_RP/max(xData(1),0.01_RP)))) + ! else + ! Sf(1)=1.0_RP + ! end if + end if + + if (xData(2).LE.0.0_RP) then + Sf(2)=0.0_RP + else if (xData(2).GE.1.0_RP) then + Sf(2)=1.0_RP + else + !if (abs(xData(2) - 1.0_RP) > epsilon(xData(2)) .and. abs(xData(2)) > epsilon(xData(2))) then + Sf(2) = 1.0_RP / (1.0_RP + exp((1.0_RP / max((xData(2) - 1.0_RP),0.01_RP)) + (1.0_RP / max(xData(2),0.01_RP)))) + ! else + ! Sf(2) = 1.0_RP + ! end if + end if + + e % storage % intensitySponge(i,j,k)= self % amplitude*(Sf(1)-Sf(2))+e % storage % intensitySponge(i,j,k) end select end do ; end do ; end do - - ! limit xStar to [0,1] since after 1 should be constant at the amplitude value - xStar = max(0.0_RP,xStar) - xStar = min(1.0_RP,xStar) - ! Sponge Ramping Function, taken from Beck, A., and Munz, C.-D., Direct Aeroacoustic Simulations Based on High Order Discontinuous Galerkin Schemes, Springer, Cham, 2018, Vol. 579 - sigma = 6.0_RP*xStar**5.0_RP - 15.0_RP*xStar**4.0_RP + 10.0_RP*xStar**3.0_RP - ! limit sigms <=1 since after 1 should be constant at the amplitude value - sigma = MIN(1.0_RP,sigma) - ! pre computed intensity, including amplitude and ramp damping - e % storage % intensitySponge(:,:,:) = max(e % storage % intensitySponge(:,:,:), self % amplitude * sigma(:,:,:)) + if((self % shapeType(sponge_number) == "cylindrical").OR.(self % shapeType(sponge_number) == "cartesian")) then + if (any(xStar .ge. 0.0_RP)) then + ! limit xStar to [0,1] since after 1 should be constant at the amplitude value + xStar = max(0.0_RP,xStar) + xStar = min(1.0_RP,xStar) + ! Sponge Ramping Function, taken from Beck, A., and Munz, C.-D., Direct Aeroacoustic Simulations Based on High Order Discontinuous Galerkin Schemes, Springer, Cham, 2018, Vol. 579 + sigma = 6.0_RP*xStar**5.0_RP - 15.0_RP*xStar**4.0_RP + 10.0_RP*xStar**3.0_RP + ! limit sigms <=1 since after 1 should be constant at the amplitude value + sigma = MIN(1.0_RP,sigma) + ! pre computed intensity, including amplitude and ramp damping + e % storage % intensitySponge(:,:,:) = max(e % storage % intensitySponge(:,:,:), self % amplitude * sigma(:,:,:)) + end if + end if end associate end do end do diff --git a/Solver/src/libs/timeintegrator/TimeIntegrator.f90 b/Solver/src/libs/timeintegrator/TimeIntegrator.f90 index 87098441fa..fddb064064 100644 --- a/Solver/src/libs/timeintegrator/TimeIntegrator.f90 +++ b/Solver/src/libs/timeintegrator/TimeIntegrator.f90 @@ -23,6 +23,7 @@ MODULE TimeIntegratorClass use MPI_Process_Info use TimeIntegratorDefinitions use MonitorsClass + use Samplings use ParticlesClass use Utilities , only: ToLower, AlmostEqual use FileReadingUtilities , only: getFileName @@ -149,9 +150,7 @@ SUBROUTINE constructTimeIntegrator(self,controlVariables, sem, initial_iter, ini self % outputInterval = controlVariables % integerValueForKey("output interval") self % tolerance = controlVariables % doublePrecisionValueForKey("convergence tolerance") self % RKStep => TakeRK3Step - - -! call sem % mesh % MLRK % construct(sem % mesh, 1) ! default 1 level + if (controlVariables % containsKey(TIME_INTEGRATION_KEY)) then self % integration_method = controlVariables % stringValueForKey(TIME_INTEGRATION_KEY, LINE_LENGTH) else @@ -359,7 +358,7 @@ END SUBROUTINE destructTimeIntegrator ! ! //////////////////////////////////////////////////////////////////////////////////////// ! - SUBROUTINE Integrate( self, sem, controlVariables, monitors, ComputeTimeDerivative, ComputeTimeDerivativeIsolated) + SUBROUTINE Integrate( self, sem, controlVariables, monitors, samplings, ComputeTimeDerivative, ComputeTimeDerivativeIsolated) USE FASMultigridClass IMPLICIT NONE ! @@ -371,6 +370,7 @@ SUBROUTINE Integrate( self, sem, controlVariables, monitors, ComputeTimeDerivati TYPE(DGSem) :: sem TYPE(FTValueDictionary) :: controlVariables class(Monitor_t) :: monitors + class(Sampling_t) :: samplings procedure(ComputeTimeDerivative_f) :: ComputeTimeDerivative procedure(ComputeTimeDerivative_f) :: ComputeTimeDerivativeIsolated @@ -388,6 +388,7 @@ SUBROUTINE Integrate( self, sem, controlVariables, monitors, ComputeTimeDerivati sem % numberOfTimeSteps = self % initial_iter if (.not. self % Compute_dt) monitors % dt_restriction = DT_FIXED + if (.not. self % Compute_dt) samplings % dt_restriction = DT_FIXED ! Measure solver time ! ------------------- @@ -427,7 +428,7 @@ SUBROUTINE Integrate( self, sem, controlVariables, monitors, ComputeTimeDerivati ! Lower the residual to 0.1 * truncation error threshold ! -> See Kompenhans et al. "Adaptation strategies for high order discontinuous Galerkin methods based on Tau-estimation." Journal of Computational Physics 306 (2016): 216-236. ! ------------------------------------------------------ - call IntegrateInTime( self, sem, controlVariables, monitors, ComputeTimeDerivative, ComputeTimeDerivativeIsolated, self % pAdaptator % reqTE*0.1_RP) + call IntegrateInTime( self, sem, controlVariables, monitors, samplings, ComputeTimeDerivative, ComputeTimeDerivativeIsolated, self % pAdaptator % reqTE*0.1_RP) end if call self % pAdaptator % pAdapt(sem,sem % numberOfTimeSteps, self % time, ComputeTimeDerivative, ComputeTimeDerivativeIsolated, controlVariables) @@ -438,7 +439,7 @@ SUBROUTINE Integrate( self, sem, controlVariables, monitors, ComputeTimeDerivati ! Finish time integration ! ----------------------- - call IntegrateInTime( self, sem, controlVariables, monitors, ComputeTimeDerivative, ComputeTimeDerivativeIsolated ) + call IntegrateInTime( self, sem, controlVariables, monitors, samplings, ComputeTimeDerivative, ComputeTimeDerivativeIsolated ) ! Measure solver time ! ------------------- @@ -453,7 +454,7 @@ END SUBROUTINE Integrate ! -> If "tolerance" is provided, the value in controlVariables is ignored. ! This is only relevant for STEADY_STATE computations. ! ------------------------------------------------------------------------ - subroutine IntegrateInTime( self, sem, controlVariables, monitors, ComputeTimeDerivative, ComputeTimeDerivativeIsolated, tolerance, CTD_linear, CTD_nonlinear) + subroutine IntegrateInTime( self, sem, controlVariables, monitors, samplings, ComputeTimeDerivative, ComputeTimeDerivativeIsolated, tolerance, CTD_linear, CTD_nonlinear) use BDFTimeIntegrator use FASMultigridClass @@ -486,6 +487,7 @@ subroutine IntegrateInTime( self, sem, controlVariables, monitors, ComputeTimeDe TYPE(DGSem) :: sem TYPE(FTValueDictionary), intent(in) :: controlVariables class(Monitor_t) :: monitors + class(Sampling_t) :: samplings procedure(ComputeTimeDerivative_f) :: ComputeTimeDerivative procedure(ComputeTimeDerivative_f) :: ComputeTimeDerivativeIsolated real(kind=RP), optional, intent(in) :: tolerance !< ? tolerance to integrate down to @@ -589,12 +591,14 @@ subroutine IntegrateInTime( self, sem, controlVariables, monitors, ComputeTimeDe maxResidual = ComputeMaxResiduals(sem % mesh) sem % maxResidual = maxval(maxResidual) call Monitors % UpdateValues( sem % mesh, t, sem % numberOfTimeSteps, maxResidual, .false., dt ) + call Samplings % UpdateValues( sem % mesh, t) call self % Display(sem % mesh, monitors, sem % numberOfTimeSteps) if (self % pAdaptator % adaptation_mode == ADAPT_DYNAMIC_TIME .and. & self % pAdaptator % nextAdaptationTime == self % time) then call self % pAdaptator % pAdapt(sem,sem % numberOfTimeSteps,t, ComputeTimeDerivative, ComputeTimeDerivativeIsolated, controlVariables) self % pAdaptator % nextAdaptationTime = self % pAdaptator % nextAdaptationTime + self % pAdaptator % time_interval + call samplings % UpdateInterp(sem % mesh) end if call monitors % WriteToFile(sem % mesh) @@ -761,6 +765,10 @@ subroutine IntegrateInTime( self, sem, controlVariables, monitors, ComputeTimeDe ! --------------- call Monitors % UpdateValues( sem % mesh, t, k+1, maxResidual, self% autosave% Autosave(k+1), dt ) ! +! Update samplings +! ---------------- + call Samplings % UpdateValues( sem % mesh, t ) +! ! Exit if the target is reached ! ----------------------------- IF (self % integratorType == STEADY_STATE) THEN @@ -796,10 +804,7 @@ subroutine IntegrateInTime( self, sem, controlVariables, monitors, ComputeTimeDe call sem % particles % inject( sem % mesh ) endif endif - endif - - #endif ! ! Print monitors @@ -810,6 +815,7 @@ subroutine IntegrateInTime( self, sem, controlVariables, monitors, ComputeTimeDe ! -------------- IF( self % pAdaptator % hasToAdapt(k+1) ) then call self % pAdaptator % pAdapt(sem,k,t, ComputeTimeDerivative, ComputeTimeDerivativeIsolated, controlVariables) + call samplings % UpdateInterp(sem % mesh) end if call self % TauEstimator % estimate(sem, k+1, t, ComputeTimeDerivative, ComputeTimeDerivativeIsolated) ! @@ -840,14 +846,15 @@ subroutine IntegrateInTime( self, sem, controlVariables, monitors, ComputeTimeDe ! Flush monitors ! -------------- call monitors % WriteToFile(sem % mesh) - + call samplings % WriteToFile(sem % mesh) sem % numberOfTimeSteps = k + 1 END DO ! -! Flush the remaining information in the monitors -! ----------------------------------------------- +! Flush the remaining information in the monitors and samplings +! ------------------------------------------------------------- if ( k .ne. 0 ) then call Monitors % writeToFile(sem % mesh, force = .true. ) + call Samplings % writeToFile(sem % mesh, force = .true. ) #if defined(NAVIERSTOKES) && (!(SPALARTALMARAS)) call sem % fwh % writeToFile( force = .TRUE. ) #endif diff --git a/Solver/src/libs/timeintegrator/pAdaptationClass.f90 b/Solver/src/libs/timeintegrator/pAdaptationClass.f90 index 746803c17a..262f7badd9 100644 --- a/Solver/src/libs/timeintegrator/pAdaptationClass.f90 +++ b/Solver/src/libs/timeintegrator/pAdaptationClass.f90 @@ -363,7 +363,7 @@ subroutine OverEnriching_Initialize(this,oID) ! ----------------------------------------------- write(STD_OUT,'(/)') call SubSection_Header("Initialize Adaptation based on Overenriching Box") - write(STD_OUT,'(30X,A,A27,A20)') "->" , "ID: " , this % ID + write(STD_OUT,'(30X,A,A27,I4)') "->" , "ID: " , this % ID write(STD_OUT,'(30X,A,A27,A20)') "->" , "x span: " , x_span write(STD_OUT,'(30X,A,A27,A20)') "->" , "y span: " , y_span write(STD_OUT,'(30X,A,A27,A20)') "->" , "z span: " , z_span diff --git a/Solver/src/libs/timeintegrator/pAdaptationClassRL.f90 b/Solver/src/libs/timeintegrator/pAdaptationClassRL.f90 index 244a8b416f..337041cef5 100644 --- a/Solver/src/libs/timeintegrator/pAdaptationClassRL.f90 +++ b/Solver/src/libs/timeintegrator/pAdaptationClassRL.f90 @@ -341,12 +341,12 @@ subroutine pAdaptation_Construct(this, controlVariables, t0, mesh) case ("rho") this % acoustic_variable = 8 case default - WRITE(STD_OUT,*) 'Not recognized acoustic variable. Using pressure by default.' this % acoustic_variable = 7 + if ( MPI_Process % isRoot ) WRITE(STD_OUT,*) 'Undefined acoustic variable. Using pressure by default.' end select else - WRITE(STD_OUT,*) 'Undefined acoustic variable. Using pressure by default.' this % acoustic_variable = 7 + if ( MPI_Process % isRoot ) WRITE(STD_OUT,*) 'Undefined acoustic variable. Using pressure by default.' end if call mesh % DefineAcousticElements(this % observer, this % acoustic_sources, this % acoustic_distance, surfacesMesh % zones)