diff --git a/Solver/src/AcousticSolver/SpatialDiscretization.f90 b/Solver/src/AcousticSolver/SpatialDiscretization.f90 index 8e23497670..115c300c95 100644 --- a/Solver/src/AcousticSolver/SpatialDiscretization.f90 +++ b/Solver/src/AcousticSolver/SpatialDiscretization.f90 @@ -208,7 +208,7 @@ END SUBROUTINE ComputeTimeDerivativeIsolated subroutine TimeDerivative_ComputeQDot( mesh , particles, t) ! use ActuatorLine, only: farm - ! use SpongeClass, only: sponge + use SpongeClass, only: sponge, addSourceSponge implicit none type(HexMesh) :: mesh type(Particles_t) :: particles @@ -332,7 +332,7 @@ subroutine TimeDerivative_ComputeQDot( mesh , particles, t) end do !$omp end do ! for the sponge, loops are in the internal subroutine as values are precalculated - ! call sponge % addSource(mesh) + call addSourceSponge(sponge,mesh) ! ! *********************** ! Now add the source term diff --git a/Solver/src/libs/sources/ActuatorLine.f90 b/Solver/src/libs/sources/ActuatorLine.f90 index 384cf2b7bb..789fdebf7c 100644 --- a/Solver/src/libs/sources/ActuatorLine.f90 +++ b/Solver/src/libs/sources/ActuatorLine.f90 @@ -45,7 +45,7 @@ module ActuatorLine real(KIND=RP), allocatable :: local_velocity_temp(:) real(KIND=RP), allocatable :: local_angle_temp(:) real(KIND=RP), allocatable :: local_root_bending(:) ! N.m - real(KIND=RP), allocatable :: gauss_epsil_delta(:) ! HO element size, for calculate force Gaussian + ! real(KIND=RP), allocatable :: gauss_epsil_delta(:) ! HO element size, for calculate force Gaussian real(KIND=RP) :: tip_c1,tip_c2 ! tip force corrections real(KIND=RP), allocatable :: local_gaussian_sum(:) ! necessary for Gaussian weighted average real(KIND=RP), allocatable :: local_Re(:) ! local Re based on local conditions and the chord of the airfoil at the blade section @@ -67,6 +67,7 @@ module ActuatorLine real(KIND=RP) :: Cp ! turbine power coef. real(KIND=RP) :: Ct ! turbine thrust coef. real(KIND=RP), allocatable :: average_conditions(:,:) ! time and blade average local variables at the blade section + real(KIND=RP) :: gauss_epsil ! force Gaussian width end type type Farm_t @@ -212,7 +213,7 @@ subroutine ConstructFarm(self, controlVariables, t0, mesh) self%turbine_t(k)%blade_t(j)%point_xyz_loc(num_blade_sections,3),self%turbine_t(k)%blade_t(j)%local_torque(num_blade_sections), & self%turbine_t(k)%blade_t(j)%local_thrust(num_blade_sections),self%turbine_t(k)%blade_t(j)%local_root_bending(num_blade_sections), & self%turbine_t(k)%blade_t(j)%local_rotor_force(num_blade_sections),self%turbine_t(k)%blade_t(j)%local_gaussian_sum(num_blade_sections), & - self%turbine_t(k)%blade_t(j)%local_Re(num_blade_sections), self%turbine_t(k)%blade_t(j)%gauss_epsil_delta(num_blade_sections) ) + self%turbine_t(k)%blade_t(j)%local_Re(num_blade_sections) ) do i=1, num_blade_sections self%turbine_t(k)%blade_t(j)%airfoil_files(i,:)=' ' @@ -290,7 +291,7 @@ subroutine ConstructFarm(self, controlVariables, t0, mesh) case (2) write(STD_OUT,'(30X,A)') 'Epsilon calculated based on element size and polynomial order' write(STD_OUT,'(30X,A,A28,F10.3)') "->", 'Constant for Epsilon: ',self%gauss_epsil - if (self%calculate_with_projection) write(STD_OUT,'(30X,A)') 'Warining, epsilon calculated using properties of element 1' + write(STD_OUT,'(30X,A)') 'Warining, epsilon calculated using properties of element 1' case default write(STD_OUT,'(30X,A,A28,ES10.3)') "->", 'Fixed Epsilon value: ',self%gauss_epsil end select @@ -413,7 +414,7 @@ subroutine ConstructFarm(self, controlVariables, t0, mesh) element_loop:do eID = 1, mesh%no_of_elements do k=1, self%num_turbines tolerance = self%tolerance_factor*self%turbine_t(k)%radius - r_square = POW2(minval(mesh%elements(eID)%geom%x(2,:,:,:))-self%turbine_t(k)%hub_cood_y) + POW2(minval(mesh%elements(eID)%geom%x(3,:,:,:))-self%turbine_t(k)%hub_cood_z) + r_square = minval(POW2(mesh%elements(eID)%geom%x(2,:,:,:)-self%turbine_t(k)%hub_cood_y)) + minval(POW2(mesh%elements(eID)%geom%x(3,:,:,:)-self%turbine_t(k)%hub_cood_z)) if( r_square <= POW2(self%turbine_t(k)%radius+tolerance) & .and. minval(mesh%elements(eID)%geom%x(1,:,:,:)) < self%turbine_t(k)%hub_cood_x+tolerance & .and. maxval(mesh%elements(eID)%geom%x(1,:,:,:)) >self%turbine_t(k)%hub_cood_x-tolerance) then @@ -428,7 +429,7 @@ subroutine ConstructFarm(self, controlVariables, t0, mesh) element_loop2:do eID = 1, mesh%no_of_elements do k=1, self%num_turbines tolerance = self%tolerance_factor*self%turbine_t(k)%radius - r_square = POW2(minval(mesh%elements(eID)%geom%x(2,:,:,:))-self%turbine_t(k)%hub_cood_y) + POW2(minval(mesh%elements(eID)%geom%x(3,:,:,:))-self%turbine_t(k)%hub_cood_z) + r_square = minval(POW2(mesh%elements(eID)%geom%x(2,:,:,:)-self%turbine_t(k)%hub_cood_y)) + minval(POW2(mesh%elements(eID)%geom%x(3,:,:,:)-self%turbine_t(k)%hub_cood_z)) if( r_square <= POW2(self%turbine_t(k)%radius+tolerance) & .and. minval(mesh%elements(eID)%geom%x(1,:,:,:)) < self%turbine_t(k)%hub_cood_x+tolerance & .and. maxval(mesh%elements(eID)%geom%x(1,:,:,:)) >self%turbine_t(k)%hub_cood_x-tolerance) then @@ -440,35 +441,44 @@ subroutine ConstructFarm(self, controlVariables, t0, mesh) end do end do element_loop2 - ! precalculate epsilon for projection mode - if (self % calculate_with_projection) then - if (MPI_Process % doMPIAction) then - if (nelem .gt. 0) then - delta_temp = (mesh % elements(elementsActuated(1)) % geom % Volume / product(mesh % elements(elementsActuated(1)) % Nxyz + 1)) ** (1.0_RP / 3.0_RP) - delta_count = 1 - else - delta_temp = 0.0_RP - delta_count = 0 - end if + select case (self % epsilon_type) + case (0) + ! EPSILON - option 1 (from file) + do k = 1, self%num_turbines + self % turbine_t(k) % gauss_epsil = self % gauss_epsil + enddo + case (2) + ! EPSILON - option 3 (k is from file) + ! eps = k*delta; k is in gauss_epsil of farm + ! precalculate delta, for now only using element 1 + if (MPI_Process % doMPIAction) then + if (nelem .gt. 0) then + delta_temp = (mesh % elements(elementsActuated(1)) % geom % Volume / product(mesh % elements(elementsActuated(1)) % Nxyz + 1)) ** (1.0_RP / 3.0_RP) + delta_count = 1 + else + delta_temp = 0.0_RP + delta_count = 0 + end if #ifdef _HAS_MPI_ - call mpi_allreduce(delta_temp, delta, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD, ierr) - call mpi_allreduce(delta_count, delta_paritions, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD, ierr) + call mpi_allreduce(delta_temp, delta, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD, ierr) + call mpi_allreduce(delta_count, delta_paritions, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD, ierr) #endif - delta = delta / real(delta_paritions,kind=RP) - else - delta = (mesh % elements(elementsActuated(1)) % geom % Volume / product(mesh % elements(elementsActuated(1)) % Nxyz + 1)) ** (1.0_RP / 3.0_RP) - end if -!$omp do schedule(runtime)private(i,j,k) - do k = 1, self%num_turbines - do j = 1, self%turbine_t(k)%num_blades - do i = 1, self%turbine_t(k)%num_blade_sections - self % turbine_t(k) % blade_t(j) % gauss_epsil_delta(i) = delta - end do - enddo - enddo -!$omp end do - end if + delta = delta / real(delta_paritions,kind=RP) + else + delta = (mesh % elements(elementsActuated(1)) % geom % Volume / product(mesh % elements(elementsActuated(1)) % Nxyz + 1)) ** (1.0_RP / 3.0_RP) + end if + do k = 1, self%num_turbines + ! e = k*delta + self % turbine_t(k) % gauss_epsil = self % gauss_epsil * delta + enddo + + case default + ! EPSILON - option 2 using Cd not implementd + do k = 1, self%num_turbines + self % turbine_t(k) % gauss_epsil = self % gauss_epsil + enddo + end select ! ! Create output files ! ------------------- @@ -556,7 +566,7 @@ subroutine UpdateFarm(self,time, mesh) !local variables integer :: ii, jj, i, j, k, kk real(kind=RP) :: dt, interp, delta_temp - logical :: found + logical :: found, allfound integer :: eID, ierr real(kind=RP), dimension(NDIM) :: x, xi real(kind=RP), dimension(NCONS) :: Q, Qtemp @@ -613,7 +623,7 @@ subroutine UpdateFarm(self,time, mesh) ! use the local Q based on the position of the actuator line point ! ---------------------------------------------------------------- ! -!$omp do schedule(runtime)private(ii,jj,kk,eID,Q,Qtemp,delta_temp,xi,found) +!$omp do schedule(runtime)private(ii,jj,kk,eID,Q,Qtemp,delta_temp,xi,found,allfound) do kk = 1, self%num_turbines do jj = 1, self%turbine_t(kk)%num_blades @@ -638,22 +648,21 @@ subroutine UpdateFarm(self,time, mesh) if (found) then ! interpolate state values in the element Qtemp = interpolateQ(mesh,eID, xi) - delta_temp = (mesh % elements(eID) % geom % Volume / product(mesh % elements(eID) % Nxyz + 1)) ** (1.0_RP / 3.0_RP) else Qtemp = 0.0_RP - delta_temp = 0.0_RP end if if ( (MPI_Process % doMPIAction) ) then #ifdef _HAS_MPI_ call mpi_allreduce(Qtemp, Q, NCONS, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD, ierr) - call mpi_allreduce(delta_temp, self%turbine_t(kk)%blade_t(jj)%gauss_epsil_delta(ii), 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD, ierr) + ! call mpi_allreduce(delta_temp, self%turbine_t(kk)%blade_t(jj)%gauss_epsil_delta(ii), 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD, ierr) + ! call mpi_allreduce(found, allfound, NCONS, MPI_LOGICAL, MPI_LOR, MPI_COMM_WORLD, ierr) + ! found = allfound #endif else Q = Qtemp - self % turbine_t(kk) % blade_t(jj) % gauss_epsil_delta(ii) = delta_temp - end if if (all(Q .eq. 0.0_RP)) then + ! if (.not. found) then print*, "Actuator line point not found in mesh, x: ", x call exit(99) end if @@ -720,8 +729,8 @@ subroutine ForcesFarm(self, mesh, time, Level) do jj = 1, self % turbine_t(kk) % num_blades do ii = 1, self % turbine_t(kk) % num_blade_sections - interp = GaussianInterpolation(self%epsilon_type, mesh % elements(eID) % geom % x(:,i,j,k), self%turbine_t(kk)%blade_t(jj)%point_xyz_loc(ii,:), & - self%turbine_t(kk)%blade_t(jj)%chord(ii), self%gauss_epsil,self%turbine_t(kk)%blade_t(jj)%gauss_epsil_delta(ii)) + interp = GaussianInterpolation( mesh % elements(eID) % geom % x(:,i,j,k), self%turbine_t(kk)%blade_t(jj)%point_xyz_loc(ii,:), & + self%turbine_t(kk)%blade_t(jj)%chord(ii), self%turbine_t(kk)%gauss_epsil) call FarmGetLocalForces(self, ii, jj, kk, mesh%elements(eID)%storage%Q(:,i,j,k), interp, local_angle, local_velocity, local_Re, local_thrust, local_rotor_force) ! minus account action-reaction effect, is the force on the fliud @@ -775,8 +784,8 @@ subroutine ForcesFarm(self, mesh, time, Level) do jj = 1, self % turbine_t(kk) % num_blades do ii = 1, self % turbine_t(kk) % num_blade_sections - interp = GaussianInterpolation(self%epsilon_type, mesh % elements(eID) % geom % x(:,i,j,k), self%turbine_t(kk)%blade_t(jj)%point_xyz_loc(ii,:), & - self%turbine_t(kk)%blade_t(jj)%chord(ii), self%gauss_epsil,self%turbine_t(kk)%blade_t(jj)%gauss_epsil_delta(ii)) + interp = GaussianInterpolation(mesh % elements(eID) % geom % x(:,i,j,k), self%turbine_t(kk)%blade_t(jj)%point_xyz_loc(ii,:), & + self%turbine_t(kk)%blade_t(jj)%chord(ii), self%turbine_t(kk)%gauss_epsil) ! minus account action-reaction effect, is the force on the fluid actuator_source(1) = actuator_source(1) - self%turbine_t(kk)%blade_t(jj)%local_thrust(ii) * interp @@ -1156,40 +1165,15 @@ End Subroutine FarmUpdateBladeForces ! !/////////////////////////////////////////////////////////////////////////////////////// ! - Function GaussianInterpolation(epsilon_type, x, x_point, chord, gauss_epsil, delta, Cd) + Function GaussianInterpolation(x, x_point, chord, gauss_epsil) implicit none - integer, intent(in) :: epsilon_type real(kind=RP), intent(in) :: x(NDIM) real(kind=RP), intent(in) :: x_point(NDIM) real(kind=RP), intent(in) :: chord real(kind=RP), intent(in) :: gauss_epsil - real(kind=RP), intent(in) :: delta - real(kind=RP), intent(in), optional :: Cd real(kind=RP) :: GaussianInterpolation - !local variables - real(kind=RP) :: epsil - - select case (epsilon_type) - case (0) -! EPSILON - option 1 (from file) - epsil = gauss_epsil - case (1) -! EPSILON - option 2 - if (present(Cd)) then - epsil = max(chord/4.0_RP,chord*Cd/2.0_RP) - else - epsil = gauss_epsil - end if - case (2) -! EPSILON - option 3 (k is from file) -! eps = k*delta; k is in gauss_epsil, gauss_epsil_delta is obtained in UpdateFarm - epsil = gauss_epsil * delta - case default - epsil = gauss_epsil - end select - - GaussianInterpolation = exp( -(POW2(x(1) - x_point(1)) + POW2(x(2) - x_point(2)) + POW2(x(3) - x_point(3))) / POW2(epsil) ) / ( POW3(epsil) * pi**(3.0_RP/2.0_RP) ) + GaussianInterpolation = exp( -(POW2(x(1) - x_point(1)) + POW2(x(2) - x_point(2)) + POW2(x(3) - x_point(3))) / POW2(gauss_epsil) ) / ( POW3(gauss_epsil) * pi**(3.0_RP/2.0_RP) ) End Function GaussianInterpolation ! diff --git a/Solver/src/libs/sources/SpongeClass.f90 b/Solver/src/libs/sources/SpongeClass.f90 index 24908bb504..fc2025a577 100644 --- a/Solver/src/libs/sources/SpongeClass.f90 +++ b/Solver/src/libs/sources/SpongeClass.f90 @@ -17,7 +17,7 @@ Module SpongeClass ! Implicit None private - public sponge, addSourceSponge + public sponge, addSourceSponge, ConstructSponge, DestructSponge, UpdateBaseFlowSponge, WriteBaseFlowSponge, creatRampSponge !definition of sponge class type sponge_t @@ -28,10 +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 :: 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 @@ -42,17 +42,6 @@ Module SpongeClass ! logical :: useMovingAverage ! to use Qbase as a moving average logical :: isActive = .false. - contains - - procedure :: construct => spongeConstruct - procedure :: destruct => spongeDestruct - procedure :: creatRamp - ! procedure :: addSource - procedure :: initializeBaseFlow - procedure :: updateBaseFlow - procedure :: readBaseFlow - procedure :: writeBaseFlow - end type sponge_t type(sponge_t) :: sponge @@ -61,7 +50,7 @@ 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" + character(len=KEYWORD_LENGTH), parameter :: SPONGE_FRINGE_NAME = "fringe" enum, bind(C) enumerator :: SPONGE_CYLINDRICAL = 1, SPONGE_CARTESIAN, FRINGE @@ -73,7 +62,7 @@ Module SpongeClass ! ! SPONGE PROCEDURES -------------------------- !///////////////////////////////////////////////////////////////////////// - Subroutine spongeConstruct(self, mesh, controlVariables) + Subroutine ConstructSponge(self, mesh, controlVariables) use FileReadingUtilities, only: getRealArrayFromString use FTValueDictionaryClass use MPI_Process_Info @@ -81,8 +70,8 @@ Subroutine spongeConstruct(self, mesh, controlVariables) use mainKeywordsModule use FileReadingUtilities, only: getFileName Implicit None - class(sponge_t) :: self - type(HexMesh), intent(inout) :: mesh + type(sponge_t) :: self + type(HexMesh), intent(inout) :: mesh type(FTValueDictionary), intent(in) :: controlVariables !local variables @@ -105,7 +94,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)) + 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 @@ -126,45 +115,45 @@ Subroutine spongeConstruct(self, mesh, controlVariables) self % radious(i) = controlVariables % getValueOrDefault(tmp,0.0_RP) endif 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, '("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 + 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 + 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 + if (useNumberedKeys) then write(tmp, '("delta fall ",I0)') i else write(tmp, '("delta fall")') @@ -199,8 +188,8 @@ Subroutine spongeConstruct(self, mesh, controlVariables) self % isActive = .true. ! create arrays and pre calculate values - call self % creatRamp(mesh) - call self % initializeBaseFlow(mesh) + call creatRampSponge(self,mesh) + call initializeBaseFlow(self,mesh) if ( .not. MPI_Process % isRoot ) return call Subsection_Header("Sponge") @@ -215,23 +204,23 @@ 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 - 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 + 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 - End Subroutine spongeConstruct + End Subroutine ConstructSponge ! - Subroutine spongeDestruct(self) + Subroutine DestructSponge(self) Implicit None class(sponge_t), intent(inout) :: self @@ -243,26 +232,26 @@ Subroutine spongeDestruct(self) safedeallocate(self % radious) safedeallocate(self % axis) safedeallocate(self % shapeType) - safedeallocate(self % fringeStart) - safedeallocate(self % fringeEnd) - safedeallocate(self % deltaRise) - safedeallocate(self % deltaFall) + safedeallocate(self % fringeStart) + safedeallocate(self % fringeEnd) + safedeallocate(self % deltaRise) + safedeallocate(self % deltaFall) - End Subroutine spongeDestruct + End Subroutine DestructSponge ! - Subroutine creatRamp(self, mesh) + Subroutine creatRampSponge(self, mesh) use MPI_Process_Info #ifdef _HAS_MPI_ use mpi #endif Implicit None - class(sponge_t) :: self + type(sponge_t) :: self type(HexMesh), intent(inout) :: mesh !local variables real(kind=RP), dimension(:,:,:), allocatable :: xStar, sigma real(kind=RP), dimension(NDIM) :: r_vector - real(kind=RP) :: xData(2), Sf(2) + real(kind=RP) :: xData(2), Sf(2) logical, dimension(:), allocatable :: hasSponge integer :: i, j, k, eID, counter, spongeEID, ierr integer, dimension(NDIM) :: Nxyz @@ -285,7 +274,7 @@ Subroutine creatRamp(self, mesh) whichSponge = SPONGE_CYLINDRICAL case (SPONGE_CARTESIAN_NAME) whichSponge = SPONGE_CARTESIAN - case (SPONGE_FRINGE_NAME) + case (SPONGE_FRINGE_NAME) whichSponge = FRINGE case default print*, "Sponge name not recognized." @@ -309,11 +298,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) + 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 + 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 @@ -353,7 +342,7 @@ Subroutine creatRamp(self, mesh) whichSponge = SPONGE_CYLINDRICAL case (SPONGE_CARTESIAN_NAME) whichSponge = SPONGE_CARTESIAN - case (SPONGE_FRINGE_NAME) + case (SPONGE_FRINGE_NAME) whichSponge = FRINGE end select @@ -376,50 +365,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) + 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 - 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 + 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 @@ -429,11 +418,11 @@ Subroutine creatRamp(self, mesh) if (allocated(sigma)) deallocate(sigma) if (allocated(hasSponge)) deallocate(hasSponge) - End Subroutine creatRamp + End Subroutine creatRampSponge ! Subroutine addSourceSponge(self,mesh) Implicit None - class(sponge_t) :: self + type(sponge_t) :: self type(HexMesh), intent(inout) :: mesh !local variables @@ -448,8 +437,6 @@ Subroutine addSourceSponge(self,mesh) eID = self % elementIndexMap(spongeEID) associate(e => mesh % elements(eID)) do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1) - ! e % storage % S_NS(:,i,j,k) = e % storage % S_NS(:,i,j,k) - e % storage % intensitySponge(i,j,k) * & - ! (e % storage % Q(:,i,j,k) - e % storage % QbaseSponge(:,i,j,k)) e % storage % S_NS(:,i,j,k) = - e % storage % intensitySponge(i,j,k) * (e % storage % Q(:,i,j,k) - e % storage % QbaseSponge(:,i,j,k)) end do ; end do ; end do end associate @@ -461,14 +448,14 @@ End Subroutine addSourceSponge Subroutine initializeBaseFlow(self,mesh) use ElementClass Implicit None - class(sponge_t) :: self + type(sponge_t) :: self type(HexMesh), intent(inout) :: mesh !local variables integer :: eID if (self % readBaseFLowFlag) then - call self % readBaseFlow(mesh) + call readBaseFlow(self,mesh) else !$omp do schedule(runtime) private(eID) do eID = 1, mesh % no_of_elements @@ -482,9 +469,9 @@ Subroutine initializeBaseFlow(self,mesh) End Subroutine initializeBaseFlow ! advance in time base flow as single euler step, taken from Beck 2018 - Subroutine updateBaseFlow(self, mesh, dt) + Subroutine UpdateBaseFlowSponge(self, mesh, dt) Implicit None - class(sponge_t) :: self + type(sponge_t) :: self type(HexMesh), intent(inout) :: mesh real(kind=RP), intent(in) :: dt @@ -507,16 +494,16 @@ Subroutine updateBaseFlow(self, mesh, dt) end do !$omp end do - End Subroutine updateBaseFlow + End Subroutine UpdateBaseFlowSponge ! Subroutine readBaseFlow(self,mesh) Implicit None - class(sponge_t) :: self + type(sponge_t) :: self type(HexMesh), intent(inout) :: mesh !local variables INTEGER :: fID, eID, fileType, no_of_elements, flag - integer :: padding, pos + integer(kind=AddrInt) :: pos, padding integer :: Nxp1, Nyp1, Nzp1, no_of_eqs, array_rank real(kind=RP), allocatable :: Q(:,:,:,:) integer, dimension(NDIM) :: Nxyz @@ -564,7 +551,7 @@ Subroutine readBaseFlow(self,mesh) Nxyz = e % Nxyz if (allocated(Q)) deallocate(Q) allocate(Q(NCONS, 0:Nxyz(1), 0:Nxyz(2), 0:Nxyz(3))) - pos = POS_INIT_DATA + (e % globID-1)*5*SIZEOF_INT + padding*e % offsetIO*SIZEOF_RP + pos = POS_INIT_DATA + (e % globID-1)*5_AddrInt*SIZEOF_INT + 1_AddrInt*padding*e % offsetIO * SIZEOF_RP read(fID, pos=pos) array_rank read(fID) no_of_eqs, Nxp1, Nyp1, Nzp1 if ( ((Nxp1-1) .ne. e % Nxyz(1)) & @@ -600,7 +587,7 @@ Subroutine readBaseFlow(self,mesh) End Subroutine readBaseFlow ! - Subroutine writeBaseFlow(self,mesh,iter,time,last) + Subroutine WriteBaseFlowSponge(self,mesh,iter,time,last) use FluidData, only: thermodynamics, refValues, dimensionless Implicit None class(sponge_t) :: self @@ -650,19 +637,22 @@ Subroutine writeBaseFlow(self,mesh,iter,time,last) refs(V_REF) = refValues % V refs(T_REF) = refValues % T refs(MACH_REF) = dimensionless % Mach - ! refs(RE_REF) = dimensionless % Re -#endif -#if defined(INCNS) +#elif defined(INCNS) refs(GAMMA_REF) = 0.0_RP refs(RGAS_REF) = 0.0_RP refs(RHO_REF) = refValues % rho refs(V_REF) = refValues % V refs(T_REF) = 0.0_RP refs(MACH_REF) = 0.0_RP - !refs(RE_REF) = dimensionless % Re !throws an erro in debug -#endif -#if defined(Multiphase) +#elif defined(Multiphase) refs = 0.0_RP +#elif defined(ACOUSTIC) + refs(GAMMA_REF) = thermodynamics % gamma + refs(RGAS_REF) = thermodynamics % R + refs(RHO_REF) = refValues % rho + refs(V_REF) = refValues % V + refs(T_REF) = refValues % T + refs(MACH_REF) = dimensionless % Mach #endif ! Create new file ! --------------- @@ -692,7 +682,7 @@ Subroutine writeBaseFlow(self,mesh,iter,time,last) ! -------------- call SealSolutionFile(trim(name)) - End Subroutine writeBaseFlow + End Subroutine WriteBaseFlowSponge ! #endif End Module SpongeClass diff --git a/Solver/src/libs/timeintegrator/TimeIntegrator.f90 b/Solver/src/libs/timeintegrator/TimeIntegrator.f90 index 8d48cd73d4..61fe680193 100644 --- a/Solver/src/libs/timeintegrator/TimeIntegrator.f90 +++ b/Solver/src/libs/timeintegrator/TimeIntegrator.f90 @@ -510,8 +510,10 @@ subroutine IntegrateInTime( self, sem, controlVariables, monitors, samplings, Co use WallFunctionDefinitions, only: useAverageV use WallFunctionConnectivity, only: Initialize_WallConnection, WallUpdateMeanV, useWallFunc #endif +#if defined(FLOW) + use SpongeClass, only: sponge, ConstructSponge, DestructSponge, UpdateBaseFlowSponge, WriteBaseFlowSponge +#endif #if defined(NAVIERSTOKES) || defined(INCNS) || defined(MULTIPHASE) - use SpongeClass, only: sponge use ActuatorLine, only: farm, ConstructFarm, DestructFarm, UpdateFarm, WriteFarmForces #endif @@ -593,7 +595,9 @@ subroutine IntegrateInTime( self, sem, controlVariables, monitors, samplings, Co call ConstructFarm(farm, controlVariables, t, sem % mesh) call UpdateFarm(farm, t, sem % mesh) end if - call sponge % construct(sem % mesh,controlVariables) +#endif +#if defined(FLOW) + call ConstructSponge(sponge,sem % mesh,controlVariables) #endif ! ! ---------------------------------- @@ -798,7 +802,9 @@ subroutine IntegrateInTime( self, sem, controlVariables, monitors, samplings, Co #if defined(NAVIERSTOKES) || defined(INCNS) || defined(MULTIPHASE) if(ActuatorLineFlag) call WriteFarmForces(farm, t, k) - call sponge % updateBaseFlow(sem % mesh,dt) +#endif +#if defined(FLOW) + call UpdateBaseFlowSponge(sponge,sem % mesh,dt) #endif ! ! Compute the new time @@ -918,7 +924,9 @@ subroutine IntegrateInTime( self, sem, controlVariables, monitors, samplings, Co #endif #if defined(NAVIERSTOKES) || defined(INCNS) || defined(MULTIPHASE) if(ActuatorLineFlag) call WriteFarmForces(farm, t, k, last=.true.) - call sponge % writeBaseFlow(sem % mesh, k, t, last=.true.) +#endif +#if defined(FLOW) + call WriteBaseFlowSponge(sponge, sem % mesh, k, t, last=.true.) #endif end if @@ -949,8 +957,10 @@ subroutine IntegrateInTime( self, sem, controlVariables, monitors, samplings, Co if (useTrip) call randomTrip % destruct #endif +#if defined(FLOW) + call DestructSponge(sponge) +#endif #if defined(NAVIERSTOKES) || defined(INCNS) || defined(MULTIPHASE) - call sponge % destruct() if(ActuatorLineFlag) call DestructFarm(farm) #endif if (saveOrders) call sem % mesh % ExportOrders(SolutionFileName) @@ -1005,8 +1015,8 @@ end subroutine TimeIntegrator_Display !///////////////////////////////////////////////////////////////////////////////////////////////// ! SUBROUTINE SaveRestart(sem,k,t,RestFileName, saveGradients, saveSensor, saveLES, saveSource) -#if defined(NAVIERSTOKES) || defined(INCNS) || defined(MULTIPHASE) - use SpongeClass, only: sponge +#if defined(FLOW) + use SpongeClass, only: sponge, WriteBaseFlowSponge #endif IMPLICIT NONE ! @@ -1031,8 +1041,8 @@ SUBROUTINE SaveRestart(sem,k,t,RestFileName, saveGradients, saveSensor, saveLES, WRITE(FinalName,'(2A,I10.10,A)') TRIM(RestFileName),'_',k,'.hsol' if ( MPI_Process % isRoot ) write(STD_OUT,'(A,A,A,ES10.3,A)') '*** Writing file "',trim(FinalName),'", with t = ',t,'.' call sem % mesh % SaveSolution(k,t,trim(finalName),saveGradients,saveSensor, saveLES, saveSource) -#if defined(NAVIERSTOKES) || defined(INCNS) || defined(MULTIPHASE) - call sponge % writeBaseFlow(sem % mesh, k, t) +#if defined(FLOW) + call WriteBaseFlowSponge(sponge, sem % mesh, k, t) #endif END SUBROUTINE SaveRestart ! diff --git a/Solver/src/libs/timeintegrator/pAdaptationClassRL.f90 b/Solver/src/libs/timeintegrator/pAdaptationClassRL.f90 index 5a4e4e5f33..a62c546d63 100644 --- a/Solver/src/libs/timeintegrator/pAdaptationClassRL.f90 +++ b/Solver/src/libs/timeintegrator/pAdaptationClassRL.f90 @@ -466,8 +466,8 @@ end subroutine pAdaptation_Destruct ! ------------------------------------------------------------------------ subroutine pAdaptation_pAdapt(this, sem, itera, t, computeTimeDerivative, ComputeTimeDerivativeIsolated, controlVariables, adaptiveTimeStep) use AnisFASMultigridClass -#if defined(NAVIERSTOKES) || defined(INCNS) || defined(MULTIPHASE) - use SpongeClass, only: sponge +#if defined(FLOW) + use SpongeClass, only: sponge, creatRampSponge #endif implicit none !-arguments---------------------------- @@ -667,8 +667,8 @@ subroutine pAdaptation_pAdapt(this, sem, itera, t, computeTimeDerivative, Comput ! Reconstruct sponge ! ---------------------------------- ! -#if defined(NAVIERSTOKES) || defined(INCNS) || defined(MULTIPHASE) - call sponge % creatRamp(sem % mesh) +#if defined(FLOW) + call creatRampSponge(sponge, sem % mesh) #endif ! Reconstruct probes @@ -1152,4 +1152,4 @@ subroutine pAdaptation_pAdaptRL_SelectElemPolorders (this, e, NNew, tolerance, N end subroutine pAdaptation_pAdaptRL_SelectElemPolorders -end module pAdaptationClassRL \ No newline at end of file +end module pAdaptationClassRL diff --git a/Solver/test/IncompressibleNS/ActuatorLineAverage/SETUP/ProblemFile.f90 b/Solver/test/IncompressibleNS/ActuatorLineAverage/SETUP/ProblemFile.f90 deleted file mode 100644 index a7dbce6e2e..0000000000 --- a/Solver/test/IncompressibleNS/ActuatorLineAverage/SETUP/ProblemFile.f90 +++ /dev/null @@ -1,638 +0,0 @@ -! -!//////////////////////////////////////////////////////////////////////// -! -! The Problem File contains user defined procedures -! that are used to "personalize" i.e. define a specific -! problem to be solved. These procedures include initial conditions, -! exact solutions (e.g. for tests), etc. and allow modifications -! without having to modify the main code. -! -! The procedures, *even if empty* that must be defined are -! -! UserDefinedSetUp -! UserDefinedInitialCondition(mesh) -! UserDefinedPeriodicOperation(mesh) -! UserDefinedFinalize(mesh) -! UserDefinedTermination -! -!//////////////////////////////////////////////////////////////////////// -! -#include "Includes.h" -module ProblemFileFunctions - implicit none - - abstract interface - subroutine UserDefinedStartup_f - end subroutine UserDefinedStartup_f - - SUBROUTINE UserDefinedFinalSetup_f(mesh & -#ifdef FLOW - , thermodynamics_ & - , dimensionless_ & - , refValues_ & -#endif -#ifdef CAHNHILLIARD - , multiphase_ & -#endif - ) - USE HexMeshClass - use FluidData - IMPLICIT NONE - CLASS(HexMesh) :: mesh -#ifdef FLOW - type(Thermodynamics_t), intent(in) :: thermodynamics_ - type(Dimensionless_t), intent(in) :: dimensionless_ - type(RefValues_t), intent(in) :: refValues_ -#endif -#ifdef CAHNHILLIARD - type(Multiphase_t), intent(in) :: multiphase_ -#endif - END SUBROUTINE UserDefinedFinalSetup_f - - subroutine UserDefinedInitialCondition_f(mesh & -#ifdef FLOW - , thermodynamics_ & - , dimensionless_ & - , refValues_ & -#endif -#ifdef CAHNHILLIARD - , multiphase_ & -#endif - ) - use smconstants - use physicsstorage - use hexmeshclass - use fluiddata - implicit none - class(hexmesh) :: mesh -#ifdef FLOW - type(Thermodynamics_t), intent(in) :: thermodynamics_ - type(Dimensionless_t), intent(in) :: dimensionless_ - type(RefValues_t), intent(in) :: refValues_ -#endif -#ifdef CAHNHILLIARD - type(Multiphase_t), intent(in) :: multiphase_ -#endif - end subroutine UserDefinedInitialCondition_f -#ifdef FLOW - subroutine UserDefinedState_f(x, t, nHat, Q, thermodynamics_, dimensionless_, refValues_) - use SMConstants - use PhysicsStorage - use FluidData - implicit none - real(kind=RP) :: x(NDIM) - real(kind=RP) :: t - real(kind=RP) :: nHat(NDIM) - real(kind=RP) :: Q(NCONS) - type(Thermodynamics_t), intent(in) :: thermodynamics_ - type(Dimensionless_t), intent(in) :: dimensionless_ - type(RefValues_t), intent(in) :: refValues_ - end subroutine UserDefinedState_f - - subroutine UserDefinedGradVars_f(x, t, nHat, Q, U, thermodynamics_, dimensionless_, refValues_) - use SMConstants - use PhysicsStorage - use FluidData - implicit none - real(kind=RP), intent(in) :: x(NDIM) - real(kind=RP), intent(in) :: t - real(kind=RP), intent(in) :: nHat(NDIM) - real(kind=RP), intent(in) :: Q(NCONS) - real(kind=RP), intent(inout) :: U(NGRAD) - type(Thermodynamics_t), intent(in) :: thermodynamics_ - type(Dimensionless_t), intent(in) :: dimensionless_ - type(RefValues_t), intent(in) :: refValues_ - end subroutine UserDefinedGradVars_f - - - subroutine UserDefinedNeumann_f(x, t, nHat, Q, U_x, U_y, U_z, flux, thermodynamics_, dimensionless_, refValues_) - use SMConstants - use PhysicsStorage - use FluidData - implicit none - real(kind=RP), intent(in) :: x(NDIM) - real(kind=RP), intent(in) :: t - real(kind=RP), intent(in) :: nHat(NDIM) - real(kind=RP), intent(in) :: Q(NCONS) - real(kind=RP), intent(in) :: U_x(NGRAD) - real(kind=RP), intent(in) :: U_y(NGRAD) - real(kind=RP), intent(in) :: U_z(NGRAD) - real(kind=RP), intent(inout) :: flux(NCONS) - type(Thermodynamics_t), intent(in) :: thermodynamics_ - type(Dimensionless_t), intent(in) :: dimensionless_ - type(RefValues_t), intent(in) :: refValues_ - end subroutine UserDefinedNeumann_f - -#endif -! -!//////////////////////////////////////////////////////////////////////// -! - SUBROUTINE UserDefinedPeriodicOperation_f(mesh, time, dt, Monitors) - use SMConstants - USE HexMeshClass - use MonitorsClass - IMPLICIT NONE - CLASS(HexMesh) :: mesh - REAL(KIND=RP) :: time - REAL(KIND=RP) :: dt - type(Monitor_t), intent(in) :: monitors - END SUBROUTINE UserDefinedPeriodicOperation_f -! -!//////////////////////////////////////////////////////////////////////// -! -#ifdef FLOW - subroutine UserDefinedSourceTermNS_f(x, Q, time, S, thermodynamics_, dimensionless_, refValues_ & -#ifdef CAHNHILLIARD -,multiphase_ & -#endif -) - use SMConstants - USE HexMeshClass - use FluidData - use PhysicsStorage - IMPLICIT NONE - real(kind=RP), intent(in) :: x(NDIM) - real(kind=RP), intent(in) :: Q(NCONS) - real(kind=RP), intent(in) :: time - real(kind=RP), intent(inout) :: S(NCONS) - type(Thermodynamics_t), intent(in) :: thermodynamics_ - type(Dimensionless_t), intent(in) :: dimensionless_ - type(RefValues_t), intent(in) :: refValues_ -#ifdef CAHNHILLIARD - type(Multiphase_t), intent(in) :: multiphase_ -#endif - end subroutine UserDefinedSourceTermNS_f -#endif -! -!//////////////////////////////////////////////////////////////////////// -! - SUBROUTINE UserDefinedFinalize_f(mesh, time, iter, maxResidual & -#ifdef FLOW - , thermodynamics_ & - , dimensionless_ & - , refValues_ & -#endif -#ifdef CAHNHILLIARD - , multiphase_ & -#endif - , monitors, & - elapsedTime, & - CPUTime ) - use SMConstants - USE HexMeshClass - use FluidData - use MonitorsClass - IMPLICIT NONE - CLASS(HexMesh) :: mesh - REAL(KIND=RP) :: time - integer :: iter - real(kind=RP) :: maxResidual -#ifdef FLOW - type(Thermodynamics_t), intent(in) :: thermodynamics_ - type(Dimensionless_t), intent(in) :: dimensionless_ - type(RefValues_t), intent(in) :: refValues_ -#endif -#ifdef CAHNHILLIARD - type(Multiphase_t), intent(in) :: multiphase_ -#endif - type(Monitor_t), intent(in) :: monitors - real(kind=RP), intent(in) :: elapsedTime - real(kind=RP), intent(in) :: CPUTime - END SUBROUTINE UserDefinedFinalize_f - - SUBROUTINE UserDefinedTermination_f - implicit none - END SUBROUTINE UserDefinedTermination_f - end interface - -end module ProblemFileFunctions - - SUBROUTINE UserDefinedStartup -! -! -------------------------------- -! Called before any other routines -! -------------------------------- -! - IMPLICIT NONE - END SUBROUTINE UserDefinedStartup -! -!//////////////////////////////////////////////////////////////////////// -! - SUBROUTINE UserDefinedFinalSetup(mesh & -#ifdef FLOW - , thermodynamics_ & - , dimensionless_ & - , refValues_ & -#endif -#ifdef CAHNHILLIARD - , multiphase_ & -#endif - ) -! -! ---------------------------------------------------------------------- -! Called after the mesh is read in to allow mesh related initializations -! or memory allocations. -! ---------------------------------------------------------------------- -! - USE HexMeshClass - use PhysicsStorage - use FluidData - IMPLICIT NONE - CLASS(HexMesh) :: mesh -#ifdef FLOW - type(Thermodynamics_t), intent(in) :: thermodynamics_ - type(Dimensionless_t), intent(in) :: dimensionless_ - type(RefValues_t), intent(in) :: refValues_ -#endif -#ifdef CAHNHILLIARD - type(Multiphase_t), intent(in) :: multiphase_ -#endif - END SUBROUTINE UserDefinedFinalSetup -! -!//////////////////////////////////////////////////////////////////////// -! - subroutine UserDefinedInitialCondition(mesh & -#ifdef FLOW - , thermodynamics_ & - , dimensionless_ & - , refValues_ & -#endif -#ifdef CAHNHILLIARD - , multiphase_ & -#endif - ) -! -! ------------------------------------------------ -! called to set the initial condition for the flow -! - by default it sets an uniform initial -! condition. -! ------------------------------------------------ -! - use smconstants - use physicsstorage - use hexmeshclass - use fluiddata - implicit none - class(hexmesh) :: mesh -#ifdef FLOW - type(Thermodynamics_t), intent(in) :: thermodynamics_ - type(Dimensionless_t), intent(in) :: dimensionless_ - type(RefValues_t), intent(in) :: refValues_ -#endif -#ifdef CAHNHILLIARD - type(Multiphase_t), intent(in) :: multiphase_ -#endif -! -! --------------- -! local variables -! --------------- -! - integer :: eid, i, j, k - real(kind=RP) :: qq, u, v, w, p -#if defined(NAVIERSTOKES) - real(kind=RP) :: Q(NCONS), phi, theta -#endif - -! -! --------------------------------------- -! Navier-Stokes default initial condition -! --------------------------------------- -! -#if defined(NAVIERSTOKES) - associate ( gammaM2 => dimensionless_ % gammaM2, & - gamma => thermodynamics_ % gamma ) - theta = refvalues_ % AOAtheta*(pi/180.0_RP) - phi = refvalues_ % AOAphi*(pi/180.0_RP) - - do eID = 1, mesh % no_of_elements - associate( Nx => mesh % elements(eID) % Nxyz(1), & - ny => mesh % elemeNts(eID) % nxyz(2), & - Nz => mesh % elements(eID) % Nxyz(3) ) - do k = 0, Nz; do j = 0, Ny; do i = 0, Nx - qq = 1.0_RP - u = qq*cos(theta)*cos(phi) - v = qq*sin(theta)*cos(phi) - w = qq*sin(phi) - - q(1) = 1.0_RP - p = 1.0_RP/(gammaM2) - q(2) = q(1)*u - q(3) = q(1)*v - q(4) = q(1)*w - q(5) = p/(gamma - 1._RP) + 0.5_RP*q(1)*(u**2 + v**2 + w**2) - - mesh % elements(eID) % storage % q(:,i,j,k) = q - end do; end do; end do - end associate - end do - - end associate -#endif -! -! ------------------------------------------------------ -! Incompressible Navier-Stokes default initial condition -! ------------------------------------------------------ -! -#if defined(INCNS) - do eID = 1, mesh % no_of_elements - associate( Nx => mesh % elements(eID) % Nxyz(1), & - ny => mesh % elemeNts(eID) % nxyz(2), & - Nz => mesh % elements(eID) % Nxyz(3) ) - do k = 0, Nz; do j = 0, Ny; do i = 0, Nx - mesh % elements(eID) % storage % q(:,i,j,k) = [1.0_RP, 1.0_RP,0.0_RP,0.0_RP,0.0_RP] - end do; end do; end do - end associate - end do -#endif - -! -! --------------------------------------- -! Cahn-Hilliard default initial condition -! --------------------------------------- -! -#ifdef CAHNHILLIARD - call random_seed() - - do eid = 1, mesh % no_of_elements - associate( Nx => mesh % elements(eid) % Nxyz(1), & - Ny => mesh % elements(eid) % Nxyz(2), & - Nz => mesh % elements(eid) % Nxyz(3) ) - associate(e => mesh % elements(eID) % storage) - call random_number(e % c) - e % c = 2.0_RP * (e % c - 0.5_RP) - end associate - end associate - end do -#endif - - end subroutine UserDefinedInitialCondition -#ifdef FLOW - subroutine UserDefinedState1(x, t, nHat, Q, thermodynamics_, dimensionless_, refValues_) - use SMConstants - use PhysicsStorage - use FluidData - implicit none - real(kind=RP), intent(in) :: x(NDIM) - real(kind=RP), intent(in) :: t - real(kind=RP), intent(in) :: nHat(NDIM) - real(kind=RP), intent(inout) :: Q(NCONS) - type(Thermodynamics_t), intent(in) :: thermodynamics_ - type(Dimensionless_t), intent(in) :: dimensionless_ - type(RefValues_t), intent(in) :: refValues_ - end subroutine UserDefinedState1 - - subroutine UserDefinedGradVars1(x, t, nHat, Q, U, thermodynamics_, dimensionless_, refValues_) - use SMConstants - use PhysicsStorage - use FluidData - implicit none - real(kind=RP), intent(in) :: x(NDIM) - real(kind=RP), intent(in) :: t - real(kind=RP), intent(in) :: nHat(NDIM) - real(kind=RP), intent(in) :: Q(NCONS) - real(kind=RP), intent(inout) :: U(NGRAD) - type(Thermodynamics_t), intent(in) :: thermodynamics_ - type(Dimensionless_t), intent(in) :: dimensionless_ - type(RefValues_t), intent(in) :: refValues_ - end subroutine UserDefinedGradVars1 - - subroutine UserDefinedNeumann1(x, t, nHat, Q, U_x, U_y, U_z, flux, thermodynamics_, dimensionless_, refValues_) - use SMConstants - use PhysicsStorage - use FluidData - implicit none - real(kind=RP), intent(in) :: x(NDIM) - real(kind=RP), intent(in) :: t - real(kind=RP), intent(in) :: nHat(NDIM) - real(kind=RP), intent(in) :: Q(NCONS) - real(kind=RP), intent(in) :: U_x(NGRAD) - real(kind=RP), intent(in) :: U_y(NGRAD) - real(kind=RP), intent(in) :: U_z(NGRAD) - real(kind=RP), intent(inout) :: flux(NCONS) - type(Thermodynamics_t), intent(in) :: thermodynamics_ - type(Dimensionless_t), intent(in) :: dimensionless_ - type(RefValues_t), intent(in) :: refValues_ - end subroutine UserDefinedNeumann1 -#endif -! -!//////////////////////////////////////////////////////////////////////// -! - SUBROUTINE UserDefinedPeriodicOperation(mesh, time, dt, Monitors) -! -! ---------------------------------------------------------- -! Called before every time-step to allow periodic operations -! to be performed -! ---------------------------------------------------------- -! - use SMConstants - USE HexMeshClass - use MonitorsClass - IMPLICIT NONE - CLASS(HexMesh) :: mesh - REAL(KIND=RP) :: time - REAL(KIND=RP) :: dt - type(Monitor_t), intent(in) :: monitors - - END SUBROUTINE UserDefinedPeriodicOperation -! -!//////////////////////////////////////////////////////////////////////// -! -#ifdef FLOW - subroutine UserDefinedSourceTermNS(x, Q, time, S, thermodynamics_, dimensionless_, refValues_ & -#ifdef CAHNHILLIARD -, multiphase_ & -#endif -) -! -! -------------------------------------------- -! Called to apply source terms to the equation -! -------------------------------------------- -! - use SMConstants - USE HexMeshClass - use PhysicsStorage - use FluidData - IMPLICIT NONE - real(kind=RP), intent(in) :: x(NDIM) - real(kind=RP), intent(in) :: Q(NCONS) - real(kind=RP), intent(in) :: time - real(kind=RP), intent(inout) :: S(NCONS) - type(Thermodynamics_t), intent(in) :: thermodynamics_ - type(Dimensionless_t), intent(in) :: dimensionless_ - type(RefValues_t), intent(in) :: refValues_ -#ifdef CAHNHILLIARD - type(Multiphase_t), intent(in) :: multiphase_ -#endif -! -! --------------- -! Local variables -! --------------- -! - integer :: i, j, k, eID -! -! Usage example -! ------------- -! S(:) = x(1) + x(2) + x(3) + time - - end subroutine UserDefinedSourceTermNS -#endif -! -!//////////////////////////////////////////////////////////////////////// -! - SUBROUTINE UserDefinedFinalize(mesh, time, iter, maxResidual & -#ifdef FLOW - , thermodynamics_ & - , dimensionless_ & - , refValues_ & -#endif -#ifdef CAHNHILLIARD - , multiphase_ & -#endif - , monitors, & - elapsedTime, & - CPUTime ) -! -! -------------------------------------------------------- -! Called after the solution computed to allow, for example -! error tests to be performed -! -------------------------------------------------------- -! - use SMConstants - use FTAssertions - USE HexMeshClass - use PhysicsStorage - use FluidData - use MonitorsClass - IMPLICIT NONE - CLASS(HexMesh) :: mesh - REAL(KIND=RP) :: time - integer :: iter - real(kind=RP) :: maxResidual -#ifdef FLOW - type(Thermodynamics_t), intent(in) :: thermodynamics_ - type(Dimensionless_t), intent(in) :: dimensionless_ - type(RefValues_t), intent(in) :: refValues_ -#endif -#ifdef CAHNHILLIARD - type(Multiphase_t), intent(in) :: multiphase_ -#endif - type(Monitor_t), intent(in) :: monitors - real(kind=RP), intent(in) :: elapsedTime - real(kind=RP), intent(in) :: CPUTime -! -! --------------- -! Local variables -! --------------- -! - CHARACTER(LEN=29) :: testName = "Actuator Line" - REAL(KIND=RP) :: maxError - REAL(KIND=RP), ALLOCATABLE :: QExpected(:,:,:,:) - INTEGER :: eID - INTEGER :: i, j, k, N - TYPE(FTAssertionsManager), POINTER :: sharedManager - LOGICAL :: success -! -! -------------------------------------------------- -! Expected Solutions: Volume integral of source term -! -------------------------------------------------- -! -#if defined(INCNS) - - real(kind=RP), parameter :: residuals(5) = [ 7.9573578444244626E-03_RP,& - 6.2238059594071782E-01_RP,& - 2.4243699814964531E-01_RP,& - 2.5626372779868756E-01_RP,& - 9.3606484796252865E+00_RP] - - real(kind=RP), parameter :: source(5) = [ 0.0000000000000000E+00_RP,& - -6.2315845438003122E-03_RP,& - 8.3992715844761168E-06_RP,& - -4.7487359880172579E-06_RP,& - 0.0000000000000000E+00_RP] - -! - CALL initializeSharedAssertionsManager - sharedManager => sharedAssertionsManager() - - CALL FTAssertEqual(expectedValue = residuals(1)+1.0_RP, & - actualValue = monitors % residuals % values(1,1)+1.0_RP, & - tol = 1.d-11, & - msg = "Continuity residual") - - CALL FTAssertEqual(expectedValue = residuals(2)+1.0_RP, & - actualValue = monitors % residuals % values(2,1)+1.0_RP, & - tol = 1.d-11, & - msg = "X-Momentum residual") - - CALL FTAssertEqual(expectedValue = residuals(3)+1.0_RP, & - actualValue = monitors % residuals % values(3,1)+1.0_RP, & - tol = 1.d-11, & - msg = "Y-Momentum residual") - - CALL FTAssertEqual(expectedValue = residuals(4)+1.0_RP, & - actualValue = monitors % residuals % values(4,1)+1.0_RP, & - tol = 1.d-11, & - msg = "Z-Momentum residual") - - CALL FTAssertEqual(expectedValue = residuals(5)+1.0_RP, & - actualValue = monitors % residuals % values(5,1)+1.0_RP, & - tol = 1.d-10, & - msg = "Energy residual") - - CALL FTAssertEqual(expectedValue = source(1)+1.0_RP, & - actualValue = monitors % volumeMonitors(1) % values(1,1)+1.0_RP, & - tol = 1.0e-11_RP, & - msg = "Continuity source") - - CALL FTAssertEqual(expectedValue = source(2)+1.0_RP, & - actualValue = monitors % volumeMonitors(1) % values(2,1)+1.0_RP, & - tol = 1.0e-11_RP, & - msg = "X-Momentum source") - - CALL FTAssertEqual(expectedValue = source(3)+1.0_RP, & - actualValue = monitors % volumeMonitors(1) % values(3,1)+1.0_RP, & - tol = 1.0e-11_RP, & - msg = "Y-Momentum source") - - CALL FTAssertEqual(expectedValue = source(4)+1.0_RP, & - actualValue = monitors % volumeMonitors(1) % values(4,1)+1.0_RP, & - tol = 1.0e-11_RP, & - msg = "Z-Momentum source") - - CALL FTAssertEqual(expectedValue = source(5)+1.0_RP, & - actualValue = monitors % volumeMonitors(1) % values(5,1)+1.0_RP, & - tol = 1.0e-11_RP, & - msg = "Pressure source") - - CALL sharedManager % summarizeAssertions(title = testName,iUnit = 6) - - IF ( sharedManager % numberOfAssertionFailures() == 0 ) THEN - WRITE(6,*) testName, " ... Passed" - WRITE(6,*) "This test case has no expected solution yet, only checks the residual after 50 iterations." - ELSE - WRITE(6,*) testName, " ... Failed" - WRITE(6,*) "NOTE: Failure is expected when the max eigenvalue procedure is changed." - WRITE(6,*) " If that is done, re-compute the expected values and modify this procedure" - error stop 99 - END IF - WRITE(6,*) - - CALL finalizeSharedAssertionsManager - CALL detachSharedAssertionsManager -#endif - - - END SUBROUTINE UserDefinedFinalize -! -!//////////////////////////////////////////////////////////////////////// -! - SUBROUTINE UserDefinedTermination -! -! ----------------------------------------------- -! Called at the the end of the main driver after -! everything else is done. -! ----------------------------------------------- -! - IMPLICIT NONE - END SUBROUTINE UserDefinedTermination - diff --git a/Solver/test/IncompressibleNS/ActuatorLineInterpolation/SETUP/ProblemFile.f90 b/Solver/test/IncompressibleNS/ActuatorLineInterpolation/SETUP/ProblemFile.f90 index 9c75424145..247f5ec179 100644 --- a/Solver/test/IncompressibleNS/ActuatorLineInterpolation/SETUP/ProblemFile.f90 +++ b/Solver/test/IncompressibleNS/ActuatorLineInterpolation/SETUP/ProblemFile.f90 @@ -538,16 +538,16 @@ SUBROUTINE UserDefinedFinalize(mesh, time, iter, maxResidual & ! #if defined(INCNS) - real(kind=RP), parameter :: residuals(5) = [ 7.9496787733853297E-03_RP,& - 6.2110979387193821E-01_RP,& - 2.4170243767877711E-01_RP,& - 2.5557185261713861E-01_RP,& - 9.3493024139426666E+00_RP] + real(kind=RP), parameter :: residuals(5) = [ 7.9548917706895768E-03_RP,& + 6.1688371769286465E-01_RP,& + 2.4146703351376803E-01_RP,& + 2.5555065598213400E-01_RP,& + 9.3633084274341147E+00_RP] real(kind=RP), parameter :: source(5) = [ 0.0000000000000000E+00_RP,& - -6.2188688155108999E-03_RP,& - 7.8674469596670224E-06_RP,& - -4.3675311615581007E-06_RP,& + -6.2957515511053745E-03_RP,& + -1.1658960895005827E-06_RP,& + 5.9379529484098832E-07_RP,& 0.0000000000000000E+00_RP] ! diff --git a/Solver/test/Multiphase/ActuatorLineAverage/SETUP/ProblemFile.f90 b/Solver/test/Multiphase/ActuatorLineAverage/SETUP/ProblemFile.f90 deleted file mode 100644 index c865b97e72..0000000000 --- a/Solver/test/Multiphase/ActuatorLineAverage/SETUP/ProblemFile.f90 +++ /dev/null @@ -1,691 +0,0 @@ -! -!//////////////////////////////////////////////////////////////////////// -! -! The Problem File contains user defined procedures -! that are used to "personalize" i.e. define a specific -! problem to be solved. These procedures include initial conditions, -! exact solutions (e.g. for tests), etc. and allow modifications -! without having to modify the main code. -! -! The procedures, *even if empty* that must be defined are -! -! UserDefinedSetUp -! UserDefinedInitialCondition(mesh) -! UserDefinedPeriodicOperation(mesh) -! UserDefinedFinalize(mesh) -! UserDefinedTermination -! -!//////////////////////////////////////////////////////////////////////// -! -#include "Includes.h" -module ProblemFileFunctions - implicit none - - abstract interface - subroutine UserDefinedStartup_f - end subroutine UserDefinedStartup_f - - SUBROUTINE UserDefinedFinalSetup_f(mesh & -#ifdef FLOW - , thermodynamics_ & - , dimensionless_ & - , refValues_ & -#endif -#ifdef CAHNHILLIARD - , multiphase_ & -#endif - ) - USE HexMeshClass - use FluidData - IMPLICIT NONE - CLASS(HexMesh) :: mesh -#ifdef FLOW - type(Thermodynamics_t), intent(in) :: thermodynamics_ - type(Dimensionless_t), intent(in) :: dimensionless_ - type(RefValues_t), intent(in) :: refValues_ -#endif -#ifdef CAHNHILLIARD - type(Multiphase_t), intent(in) :: multiphase_ -#endif - END SUBROUTINE UserDefinedFinalSetup_f - - subroutine UserDefinedInitialCondition_f(mesh & -#ifdef FLOW - , thermodynamics_ & - , dimensionless_ & - , refValues_ & -#endif -#ifdef CAHNHILLIARD - , multiphase_ & -#endif - ) - use smconstants - use physicsstorage - use hexmeshclass - use fluiddata - implicit none - class(hexmesh) :: mesh -#ifdef FLOW - type(Thermodynamics_t), intent(in) :: thermodynamics_ - type(Dimensionless_t), intent(in) :: dimensionless_ - type(RefValues_t), intent(in) :: refValues_ -#endif -#ifdef CAHNHILLIARD - type(Multiphase_t), intent(in) :: multiphase_ -#endif - end subroutine UserDefinedInitialCondition_f -#ifdef FLOW - subroutine UserDefinedState_f(x, t, nHat, Q, thermodynamics_, dimensionless_, refValues_) - use SMConstants - use PhysicsStorage - use FluidData - implicit none - real(kind=RP) :: x(NDIM) - real(kind=RP) :: t - real(kind=RP) :: nHat(NDIM) - real(kind=RP) :: Q(NCONS) - type(Thermodynamics_t), intent(in) :: thermodynamics_ - type(Dimensionless_t), intent(in) :: dimensionless_ - type(RefValues_t), intent(in) :: refValues_ - end subroutine UserDefinedState_f - - subroutine UserDefinedGradVars_f(x, t, nHat, Q, U, thermodynamics_, dimensionless_, refValues_) - use SMConstants - use PhysicsStorage - use FluidData - implicit none - real(kind=RP), intent(in) :: x(NDIM) - real(kind=RP), intent(in) :: t - real(kind=RP), intent(in) :: nHat(NDIM) - real(kind=RP), intent(in) :: Q(NCONS) - real(kind=RP), intent(inout) :: U(NGRAD) - type(Thermodynamics_t), intent(in) :: thermodynamics_ - type(Dimensionless_t), intent(in) :: dimensionless_ - type(RefValues_t), intent(in) :: refValues_ - end subroutine UserDefinedGradVars_f - - - subroutine UserDefinedNeumann_f(x, t, nHat, Q, U_x, U_y, U_z, flux, thermodynamics_, dimensionless_, refValues_) - use SMConstants - use PhysicsStorage - use FluidData - implicit none - real(kind=RP), intent(in) :: x(NDIM) - real(kind=RP), intent(in) :: t - real(kind=RP), intent(in) :: nHat(NDIM) - real(kind=RP), intent(in) :: Q(NCONS) - real(kind=RP), intent(in) :: U_x(NGRAD) - real(kind=RP), intent(in) :: U_y(NGRAD) - real(kind=RP), intent(in) :: U_z(NGRAD) - real(kind=RP), intent(inout) :: flux(NCONS) - type(Thermodynamics_t), intent(in) :: thermodynamics_ - type(Dimensionless_t), intent(in) :: dimensionless_ - type(RefValues_t), intent(in) :: refValues_ - end subroutine UserDefinedNeumann_f - -#endif -! -!//////////////////////////////////////////////////////////////////////// -! - SUBROUTINE UserDefinedPeriodicOperation_f(mesh, time, dt, Monitors) - use SMConstants - USE HexMeshClass - use MonitorsClass - IMPLICIT NONE - CLASS(HexMesh) :: mesh - REAL(KIND=RP) :: time - REAL(KIND=RP) :: dt - type(Monitor_t), intent(in) :: monitors - END SUBROUTINE UserDefinedPeriodicOperation_f -! -!//////////////////////////////////////////////////////////////////////// -! -#ifdef FLOW - subroutine UserDefinedSourceTermNS_f(x, Q, time, S, thermodynamics_, dimensionless_, refValues_ & -#ifdef CAHNHILLIARD -,multiphase_ & -#endif -) - use SMConstants - USE HexMeshClass - use FluidData - use PhysicsStorage - IMPLICIT NONE - real(kind=RP), intent(in) :: x(NDIM) - real(kind=RP), intent(in) :: Q(NCONS) - real(kind=RP), intent(in) :: time - real(kind=RP), intent(inout) :: S(NCONS) - type(Thermodynamics_t), intent(in) :: thermodynamics_ - type(Dimensionless_t), intent(in) :: dimensionless_ - type(RefValues_t), intent(in) :: refValues_ -#ifdef CAHNHILLIARD - type(Multiphase_t), intent(in) :: multiphase_ -#endif - end subroutine UserDefinedSourceTermNS_f -#endif -! -!//////////////////////////////////////////////////////////////////////// -! - SUBROUTINE UserDefinedFinalize_f(mesh, time, iter, maxResidual & -#ifdef FLOW - , thermodynamics_ & - , dimensionless_ & - , refValues_ & -#endif -#ifdef CAHNHILLIARD - , multiphase_ & -#endif - , monitors, & - elapsedTime, & - CPUTime ) - use SMConstants - USE HexMeshClass - use FluidData - use MonitorsClass - IMPLICIT NONE - CLASS(HexMesh) :: mesh - REAL(KIND=RP) :: time - integer :: iter - real(kind=RP) :: maxResidual -#ifdef FLOW - type(Thermodynamics_t), intent(in) :: thermodynamics_ - type(Dimensionless_t), intent(in) :: dimensionless_ - type(RefValues_t), intent(in) :: refValues_ -#endif -#ifdef CAHNHILLIARD - type(Multiphase_t), intent(in) :: multiphase_ -#endif - type(Monitor_t), intent(in) :: monitors - real(kind=RP), intent(in) :: elapsedTime - real(kind=RP), intent(in) :: CPUTime - END SUBROUTINE UserDefinedFinalize_f - - SUBROUTINE UserDefinedTermination_f - implicit none - END SUBROUTINE UserDefinedTermination_f - end interface - -end module ProblemFileFunctions - - SUBROUTINE UserDefinedStartup -! -! -------------------------------- -! Called before any other routines -! -------------------------------- -! - IMPLICIT NONE - END SUBROUTINE UserDefinedStartup -! -!//////////////////////////////////////////////////////////////////////// -! - SUBROUTINE UserDefinedFinalSetup(mesh & -#ifdef FLOW - , thermodynamics_ & - , dimensionless_ & - , refValues_ & -#endif -#ifdef CAHNHILLIARD - , multiphase_ & -#endif - ) -! -! ---------------------------------------------------------------------- -! Called after the mesh is read in to allow mesh related initializations -! or memory allocations. -! ---------------------------------------------------------------------- -! - USE HexMeshClass - use PhysicsStorage - use FluidData - IMPLICIT NONE - CLASS(HexMesh) :: mesh -#ifdef FLOW - type(Thermodynamics_t), intent(in) :: thermodynamics_ - type(Dimensionless_t), intent(in) :: dimensionless_ - type(RefValues_t), intent(in) :: refValues_ -#endif -#ifdef CAHNHILLIARD - type(Multiphase_t), intent(in) :: multiphase_ -#endif - END SUBROUTINE UserDefinedFinalSetup -! -!//////////////////////////////////////////////////////////////////////// -! - subroutine UserDefinedInitialCondition(mesh & -#ifdef FLOW - , thermodynamics_ & - , dimensionless_ & - , refValues_ & -#endif -#ifdef CAHNHILLIARD - , multiphase_ & -#endif - ) -! -! ------------------------------------------------ -! called to set the initial condition for the flow -! - by default it sets an uniform initial -! condition. -! ------------------------------------------------ -! - use smconstants - use physicsstorage - use hexmeshclass - use fluiddata - implicit none - class(hexmesh) :: mesh -#ifdef FLOW - type(Thermodynamics_t), intent(in) :: thermodynamics_ - type(Dimensionless_t), intent(in) :: dimensionless_ - type(RefValues_t), intent(in) :: refValues_ -#endif -#ifdef CAHNHILLIARD - type(Multiphase_t), intent(in) :: multiphase_ -#endif -! -! --------------- -! local variables -! --------------- -! - integer :: eid, i, j, k - integer :: Nx,Ny,Nz - real(kind=RP) :: qq, u, v, w, p, c -#if defined(NAVIERSTOKES) - real(kind=RP) :: Q(NCONS), phi, theta -#endif - -! -! --------------------------------------- -! Navier-Stokes default initial condition -! --------------------------------------- -! -#if defined(NAVIERSTOKES) - associate ( gammaM2 => dimensionless_ % gammaM2, & - gamma => thermodynamics_ % gamma ) - theta = refvalues_ % AOAtheta*(pi/180.0_RP) - phi = refvalues_ % AOAphi*(pi/180.0_RP) - - do eID = 1, mesh % no_of_elements - associate( Nx => mesh % elements(eID) % Nxyz(1), & - ny => mesh % elemeNts(eID) % nxyz(2), & - Nz => mesh % elements(eID) % Nxyz(3) ) - do k = 0, Nz; do j = 0, Ny; do i = 0, Nx - qq = 1.0_RP - u = qq*cos(theta)*cos(phi) - v = qq*sin(theta)*cos(phi) - w = qq*sin(phi) - - q(1) = 1.0_RP - p = 1.0_RP/(gammaM2) - q(2) = q(1)*u - q(3) = q(1)*v - q(4) = q(1)*w - q(5) = p/(gamma - 1._RP) + 0.5_RP*q(1)*(u**2 + v**2 + w**2) - - mesh % elements(eID) % storage % q(:,i,j,k) = q - end do; end do; end do - end associate - end do - - end associate -#endif -! -! ------------------------------------------------------ -! Incompressible Navier-Stokes default initial condition -! ------------------------------------------------------ -! -#if defined(INCNS) - do eID = 1, mesh % no_of_elements - associate( Nx => mesh % elements(eID) % Nxyz(1), & - ny => mesh % elemeNts(eID) % nxyz(2), & - Nz => mesh % elements(eID) % Nxyz(3) ) - do k = 0, Nz; do j = 0, Ny; do i = 0, Nx - mesh % elements(eID) % storage % q(:,i,j,k) = [1.0_RP, 1.0_RP,0.0_RP,0.0_RP,0.0_RP] - end do; end do; end do - end associate - end do -#endif - -! -! --------------------------------------- -! Cahn-Hilliard default initial condition -! --------------------------------------- -! -#ifdef CAHNHILLIARD - call random_seed() - - do eid = 1, mesh % no_of_elements - associate( Nx => mesh % elements(eid) % Nxyz(1), & - Ny => mesh % elements(eid) % Nxyz(2), & - Nz => mesh % elements(eid) % Nxyz(3) ) - associate(e => mesh % elements(eID) % storage) - call random_number(e % c) - e % c = 2.0_RP * (e % c - 0.5_RP) - end associate - end associate - end do -#endif - -#if defined(MULTIPHASE) - - - DO eID = 1, SIZE(mesh % elements) - Nx = mesh % elements(eID) % Nxyz(1) - Ny = mesh % elements(eID) % Nxyz(2) - Nz = mesh % elements(eID) % Nxyz(3) - - DO k = 0, Nz - DO j = 0, Ny - DO i = 0, Nx - - ! c = 1.0 - 0.5*(1.0+tanh(2.0*(( mesh % elements(eID) % geom % x(IX,i,j,k) + 0.0))/multiphase_ % eps)) - ! single fluid - c = 0.0_RP - !rho= dimensionless_ % rho(1)*c + dimensionless_ % rho(2)*(1-c) - - mesh % elements(eID) % storage % Q(1,i,j,k) = c - mesh % elements(eID) % storage % Q(2,i,j,k) = 1.0_RP - mesh % elements(eID) % storage % Q(3,i,j,k) = 0.0_RP - mesh % elements(eID) % storage % Q(4,i,j,k) = 0.0_RP - mesh % elements(eID) % storage % Q(5,i,j,k) = 0.0_RP - - END DO - END DO - END DO - - END DO - -#endif - - end subroutine UserDefinedInitialCondition -#ifdef FLOW - subroutine UserDefinedState1(x, t, nHat, Q, thermodynamics_, dimensionless_, refValues_) - use SMConstants - use PhysicsStorage - use FluidData - implicit none - real(kind=RP), intent(in) :: x(NDIM) - real(kind=RP), intent(in) :: t - real(kind=RP), intent(in) :: nHat(NDIM) - real(kind=RP), intent(inout) :: Q(NCONS) - type(Thermodynamics_t), intent(in) :: thermodynamics_ - type(Dimensionless_t), intent(in) :: dimensionless_ - type(RefValues_t), intent(in) :: refValues_ - - real(kind=RP) :: c, u, v, w, p, rho, sqrtRho -#ifdef MULTIPHASE - - ! single fluid - c = 0.0_RP - - rho = dimensionless_ % rho(1)*c + dimensionless_ % rho(2)*(1.0_RP - c) - - sqrtRho = sqrt(rho) - - u = 1.0_RP * c - v = 0.0_RP - w = 0.0_RP - - ! w = Vmax1 * c + Vmax2 * (1.0_RP - c) - - p = 0.0_RP - - Q = [c,sqrtRho*u, sqrtRho*v, sqrtRho*w, p] - -#endif - end subroutine UserDefinedState1 - - subroutine UserDefinedGradVars1(x, t, nHat, Q, U, thermodynamics_, dimensionless_, refValues_) - use SMConstants - use PhysicsStorage - use FluidData - implicit none - real(kind=RP), intent(in) :: x(NDIM) - real(kind=RP), intent(in) :: t - real(kind=RP), intent(in) :: nHat(NDIM) - real(kind=RP), intent(in) :: Q(NCONS) - real(kind=RP), intent(inout) :: U(NGRAD) - type(Thermodynamics_t), intent(in) :: thermodynamics_ - type(Dimensionless_t), intent(in) :: dimensionless_ - type(RefValues_t), intent(in) :: refValues_ - end subroutine UserDefinedGradVars1 - - subroutine UserDefinedNeumann1(x, t, nHat, Q, U_x, U_y, U_z, flux, thermodynamics_, dimensionless_, refValues_) - use SMConstants - use PhysicsStorage - use FluidData - implicit none - real(kind=RP), intent(in) :: x(NDIM) - real(kind=RP), intent(in) :: t - real(kind=RP), intent(in) :: nHat(NDIM) - real(kind=RP), intent(in) :: Q(NCONS) - real(kind=RP), intent(in) :: U_x(NGRAD) - real(kind=RP), intent(in) :: U_y(NGRAD) - real(kind=RP), intent(in) :: U_z(NGRAD) - real(kind=RP), intent(inout) :: flux(NCONS) - type(Thermodynamics_t), intent(in) :: thermodynamics_ - type(Dimensionless_t), intent(in) :: dimensionless_ - type(RefValues_t), intent(in) :: refValues_ - end subroutine UserDefinedNeumann1 -#endif -! -!//////////////////////////////////////////////////////////////////////// -! - SUBROUTINE UserDefinedPeriodicOperation(mesh, time, dt, Monitors) -! -! ---------------------------------------------------------- -! Called before every time-step to allow periodic operations -! to be performed -! ---------------------------------------------------------- -! - use SMConstants - USE HexMeshClass - use MonitorsClass - IMPLICIT NONE - CLASS(HexMesh) :: mesh - REAL(KIND=RP) :: time - REAL(KIND=RP) :: dt - type(Monitor_t), intent(in) :: monitors - - END SUBROUTINE UserDefinedPeriodicOperation -! -!//////////////////////////////////////////////////////////////////////// -! -#ifdef FLOW - subroutine UserDefinedSourceTermNS(x, Q, time, S, thermodynamics_, dimensionless_, refValues_ & -#ifdef CAHNHILLIARD -, multiphase_ & -#endif -) -! -! -------------------------------------------- -! Called to apply source terms to the equation -! -------------------------------------------- -! - use SMConstants - USE HexMeshClass - use PhysicsStorage - use FluidData - IMPLICIT NONE - real(kind=RP), intent(in) :: x(NDIM) - real(kind=RP), intent(in) :: Q(NCONS) - real(kind=RP), intent(in) :: time - real(kind=RP), intent(inout) :: S(NCONS) - type(Thermodynamics_t), intent(in) :: thermodynamics_ - type(Dimensionless_t), intent(in) :: dimensionless_ - type(RefValues_t), intent(in) :: refValues_ -#ifdef CAHNHILLIARD - type(Multiphase_t), intent(in) :: multiphase_ -#endif -! -! --------------- -! Local variables -! --------------- -! - integer :: i, j, k, eID -! -! Usage example -! ------------- -! S(:) = x(1) + x(2) + x(3) + time - - end subroutine UserDefinedSourceTermNS -#endif -! -!//////////////////////////////////////////////////////////////////////// -! - SUBROUTINE UserDefinedFinalize(mesh, time, iter, maxResidual & -#ifdef FLOW - , thermodynamics_ & - , dimensionless_ & - , refValues_ & -#endif -#ifdef CAHNHILLIARD - , multiphase_ & -#endif - , monitors, & - elapsedTime, & - CPUTime ) -! -! -------------------------------------------------------- -! Called after the solution computed to allow, for example -! error tests to be performed -! -------------------------------------------------------- -! - use SMConstants - use FTAssertions - USE HexMeshClass - use PhysicsStorage - use FluidData - use MonitorsClass - IMPLICIT NONE - CLASS(HexMesh) :: mesh - REAL(KIND=RP) :: time - integer :: iter - real(kind=RP) :: maxResidual -#ifdef FLOW - type(Thermodynamics_t), intent(in) :: thermodynamics_ - type(Dimensionless_t), intent(in) :: dimensionless_ - type(RefValues_t), intent(in) :: refValues_ -#endif -#ifdef CAHNHILLIARD - type(Multiphase_t), intent(in) :: multiphase_ -#endif - type(Monitor_t), intent(in) :: monitors - real(kind=RP), intent(in) :: elapsedTime - real(kind=RP), intent(in) :: CPUTime -! -! --------------- -! Local variables -! --------------- -! - CHARACTER(LEN=29) :: testName = "Actuator Line" - REAL(KIND=RP) :: maxError - REAL(KIND=RP), ALLOCATABLE :: QExpected(:,:,:,:) - INTEGER :: eID - INTEGER :: i, j, k, N - TYPE(FTAssertionsManager), POINTER :: sharedManager - LOGICAL :: success -! -! -------------------------------------------------- -! Expected Solutions: Volume integral of source term -! -------------------------------------------------- -! -#if defined(MULTIPHASE) - - real(kind=RP), parameter :: residuals(5) = [ 0.0000000000000000E+00_RP,& - 6.6814255992417429E+01_RP,& - 2.7045594375426185E-01_RP,& - 2.9162963953534043E-01_RP,& - 2.2437631252183919E+03_RP] - - real(kind=RP), parameter :: source(5) = [ 0.0000000000000000E+00_RP,& - -5.9895760348747610E-03_RP,& - 1.9436694053209671E-05_RP,& - -1.9856398207925524E-05_RP,& - 0.0000000000000000E+00_RP] -! - CALL initializeSharedAssertionsManager - sharedManager => sharedAssertionsManager() - - CALL FTAssertEqual(expectedValue = residuals(1)+1.0_RP, & - actualValue = monitors % residuals % values(1,1)+1.0_RP, & - tol = 1.d-11, & - msg = "Continuity residual") - - CALL FTAssertEqual(expectedValue = residuals(2)+1.0_RP, & - actualValue = monitors % residuals % values(2,1)+1.0_RP, & - tol = 1.d-11, & - msg = "X-Momentum residual") - - CALL FTAssertEqual(expectedValue = residuals(3)+1.0_RP, & - actualValue = monitors % residuals % values(3,1)+1.0_RP, & - tol = 1.d-11, & - msg = "Y-Momentum residual") - - CALL FTAssertEqual(expectedValue = residuals(4)+1.0_RP, & - actualValue = monitors % residuals % values(4,1)+1.0_RP, & - tol = 1.d-11, & - msg = "Z-Momentum residual") - - CALL FTAssertEqual(expectedValue = residuals(5)+1.0_RP, & - actualValue = monitors % residuals % values(5,1)+1.0_RP, & - tol = 1.d-10, & - msg = "Energy residual") - - CALL FTAssertEqual(expectedValue = source(1)+1.0_RP, & - actualValue = monitors % volumeMonitors(1) % values(1,1)+1.0_RP, & - tol = 1.0e-11_RP, & - msg = "Continuity source") - - CALL FTAssertEqual(expectedValue = source(2)+1.0_RP, & - actualValue = monitors % volumeMonitors(1) % values(2,1)+1.0_RP, & - tol = 1.0e-11_RP, & - msg = "X-Momentum source") - - CALL FTAssertEqual(expectedValue = source(3)+1.0_RP, & - actualValue = monitors % volumeMonitors(1) % values(3,1)+1.0_RP, & - tol = 1.0e-11_RP, & - msg = "Y-Momentum source") - - CALL FTAssertEqual(expectedValue = source(4)+1.0_RP, & - actualValue = monitors % volumeMonitors(1) % values(4,1)+1.0_RP, & - tol = 1.0e-11_RP, & - msg = "Z-Momentum source") - - CALL FTAssertEqual(expectedValue = source(5)+1.0_RP, & - actualValue = monitors % volumeMonitors(1) % values(5,1)+1.0_RP, & - tol = 1.0e-11_RP, & - msg = "Pressure source") - - CALL sharedManager % summarizeAssertions(title = testName,iUnit = 6) - - IF ( sharedManager % numberOfAssertionFailures() == 0 ) THEN - WRITE(6,*) testName, " ... Passed" - WRITE(6,*) "This test case has no expected solution yet, only checks the residual after 50 iterations." - ELSE - WRITE(6,*) testName, " ... Failed" - WRITE(6,*) "NOTE: Failure is expected when the max eigenvalue procedure is changed." - WRITE(6,*) " If that is done, re-compute the expected values and modify this procedure" - error stop 99 - END IF - WRITE(6,*) - - CALL finalizeSharedAssertionsManager - CALL detachSharedAssertionsManager -#endif - - - END SUBROUTINE UserDefinedFinalize -! -!//////////////////////////////////////////////////////////////////////// -! - SUBROUTINE UserDefinedTermination -! -! ----------------------------------------------- -! Called at the the end of the main driver after -! everything else is done. -! ----------------------------------------------- -! - IMPLICIT NONE - END SUBROUTINE UserDefinedTermination - diff --git a/Solver/test/Multiphase/ActuatorLineInterpolation/SETUP/ProblemFile.f90 b/Solver/test/Multiphase/ActuatorLineInterpolation/SETUP/ProblemFile.f90 index c448ef0b2f..0d95d89076 100644 --- a/Solver/test/Multiphase/ActuatorLineInterpolation/SETUP/ProblemFile.f90 +++ b/Solver/test/Multiphase/ActuatorLineInterpolation/SETUP/ProblemFile.f90 @@ -594,14 +594,14 @@ SUBROUTINE UserDefinedFinalize(mesh, time, iter, maxResidual & real(kind=RP), parameter :: residuals(5) = [ 0.0000000000000000E+00_RP,& 6.6848885010291738E+01_RP,& - 0.2700343177172555E+00_RP,& - 0.2914828786903740E+00_RP,& + 2.6691099501763155E-01_RP,& + 2.9148212224592934E-01_RP,& 2.2442342918587788E+03_RP] real(kind=RP), parameter :: source(5) = [ 0.0000000000000000E+00_RP,& - -5.98253977235252E-03_RP,& - 0.0000192653654527E+00_RP,& - -1.9887137428830E-05_RP,& + -6.2203300268581108E-03_RP,& + 9.7418137119553618E-08_RP,& + -1.2709186043159373E-07_RP,& 0.0000000000000000E+00_RP] ! CALL initializeSharedAssertionsManager diff --git a/Solver/test/NavierStokes/ActuatorLineInterpolation/SETUP/ProblemFile.f90 b/Solver/test/NavierStokes/ActuatorLineInterpolation/SETUP/ProblemFile.f90 index 639ae95176..322e35180b 100644 --- a/Solver/test/NavierStokes/ActuatorLineInterpolation/SETUP/ProblemFile.f90 +++ b/Solver/test/NavierStokes/ActuatorLineInterpolation/SETUP/ProblemFile.f90 @@ -538,16 +538,16 @@ SUBROUTINE UserDefinedFinalize(mesh, time, iter, maxResidual & ! #if defined(NAVIERSTOKES) - real(kind=RP), parameter :: residuals(5) = [ 8.3181093712812824E-03_RP,& - 6.0521912591069082E-01_RP,& - 2.3307309187131547E-01_RP,& - 2.5823932578223946E-01_RP,& - 2.3829928894542491E+01_RP] + real(kind=RP), parameter :: residuals(5) = [ 8.3335278353326506E-03_RP,& + 6.0509613944336882E-01_RP,& + 2.3289244686369220E-01_RP,& + 2.5854985032385236E-01_RP,& + 2.3874448934744198E+01_RP] real(kind=RP), parameter :: source(5) = [ 0.0000000000000000E+00_RP,& - -6.2212865110828552E-03_RP,& - 7.0018589188196696E-06_RP,& - -3.6995113997700409E-06_RP,& + -6.2981548797418864E-03_RP,& + -2.0524924862822906E-06_RP,& + 1.2742567523920291E-06_RP,& 0.0000000000000000E+00_RP] ! CALL initializeSharedAssertionsManager diff --git a/Solver/test/NavierStokes/ActuatorLineProjection/SETUP/ProblemFile.f90 b/Solver/test/NavierStokes/ActuatorLineProjection/SETUP/ProblemFile.f90 index 37481edebf..23fe409128 100644 --- a/Solver/test/NavierStokes/ActuatorLineProjection/SETUP/ProblemFile.f90 +++ b/Solver/test/NavierStokes/ActuatorLineProjection/SETUP/ProblemFile.f90 @@ -538,16 +538,16 @@ SUBROUTINE UserDefinedFinalize(mesh, time, iter, maxResidual & ! #if defined(NAVIERSTOKES) - real(kind=RP), parameter :: residuals(5) = [ 8.3294147062644638E-03_RP,& - 6.0430810336897400E-01_RP,& - 2.3415808821048351E-01_RP,& - 2.5918707902588395E-01_RP,& - 2.3862335992216614E+01_RP] + real(kind=RP), parameter :: residuals(5) = [ 8.3448348265495256E-03_RP,& + 6.0418180386378295E-01_RP,& + 2.3396790810953455E-01_RP,& + 2.5949374476580855E-01_RP,& + 2.3906860845643735E+01_RP] real(kind=RP), parameter :: source(5) = [ 0.0000000000000000E+00_RP,& - -6.2505986327932193E-03_RP,& - 8.0633445199070679E-06_RP,& - -4.7800389962274271E-06_RP,& + -6.3277328483268723E-03_RP,& + -1.1490645475318143E-06_RP,& + 2.9094005738307613E-07_RP,& 0.0000000000000000E+00_RP] !