Skip to content

Omarino dev #268

New issue

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

Sign up for GitHub

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

Already on GitHub? Sign in to your account

Merged
oscarmarino merged 6 commits into from
Jul 31, 2025
Merged

Omarino dev #268

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
ca6936b
Fix AL getting of elementsActuated, to include minor radious
oscarmarino Jul 14, 2025
27e672b
Precalculate epsilon in AL, use it as constant per turbine in the whole
oscarmarino Jul 29, 2025
2968d26
Fix AL getting of elementsActuated, to include minor radious
oscarmarino Jul 14, 2025
53ee674
Add sponge to CAA solver and change it to be gpu compatible, removing
oscarmarino Jul 30, 2025
88a74c6
Merge branch 'master' into omarinoDev
oscarmarino Jul 30, 2025
d323d6b
Update AL tests
oscarmarino Jul 30, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions Solver/src/AcousticSolver/SpatialDiscretization.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
122 changes: 53 additions & 69 deletions Solver/src/libs/sources/ActuatorLine.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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,:)=' '
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
! -------------------
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
!
Expand Down
Loading
Loading