Skip to content

Commit ba26917

Browse files
oscarmarinooscarmarino
authored
Merge pull request #268 from loganoz/omarinoDev
Omarino dev
2 parents b241433 + d323d6b commit ba26917

11 files changed

Lines changed: 232 additions & 1577 deletions

File tree

‎Solver/src/AcousticSolver/SpatialDiscretization.f90‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -208,7 +208,7 @@ END SUBROUTINE ComputeTimeDerivativeIsolated
208208

209209
subroutine TimeDerivative_ComputeQDot( mesh , particles, t)
210210
! use ActuatorLine, only: farm
211-
! use SpongeClass, only: sponge
211+
use SpongeClass, only: sponge, addSourceSponge
212212
implicit none
213213
type(HexMesh) :: mesh
214214
type(Particles_t) :: particles
@@ -332,7 +332,7 @@ subroutine TimeDerivative_ComputeQDot( mesh , particles, t)
332332
end do
333333
!$omp end do
334334
! for the sponge, loops are in the internal subroutine as values are precalculated
335-
! call sponge % addSource(mesh)
335+
call addSourceSponge(sponge,mesh)
336336
!
337337
! ***********************
338338
! Now add the source term

‎Solver/src/libs/sources/ActuatorLine.f90‎

Lines changed: 53 additions & 69 deletions
Original file line numberDiff line numberDiff line change
@@ -45,7 +45,7 @@ module ActuatorLine
4545
real(KIND=RP), allocatable :: local_velocity_temp(:)
4646
real(KIND=RP), allocatable :: local_angle_temp(:)
4747
real(KIND=RP), allocatable :: local_root_bending(:) ! N.m
48-
real(KIND=RP), allocatable :: gauss_epsil_delta(:) ! HO element size, for calculate force Gaussian
48+
! real(KIND=RP), allocatable :: gauss_epsil_delta(:) ! HO element size, for calculate force Gaussian
4949
real(KIND=RP) :: tip_c1,tip_c2 ! tip force corrections
5050
real(KIND=RP), allocatable :: local_gaussian_sum(:) ! necessary for Gaussian weighted average
5151
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
6767
real(KIND=RP) :: Cp ! turbine power coef.
6868
real(KIND=RP) :: Ct ! turbine thrust coef.
6969
real(KIND=RP), allocatable :: average_conditions(:,:) ! time and blade average local variables at the blade section
70+
real(KIND=RP) :: gauss_epsil ! force Gaussian width
7071
end type
7172

7273
type Farm_t
@@ -212,7 +213,7 @@ subroutine ConstructFarm(self, controlVariables, t0, mesh)
212213
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), &
213214
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), &
214215
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), &
215-
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) )
216+
self%turbine_t(k)%blade_t(j)%local_Re(num_blade_sections) )
216217

217218
do i=1, num_blade_sections
218219
self%turbine_t(k)%blade_t(j)%airfoil_files(i,:)=' '
@@ -290,7 +291,7 @@ subroutine ConstructFarm(self, controlVariables, t0, mesh)
290291
case (2)
291292
write(STD_OUT,'(30X,A)') 'Epsilon calculated based on element size and polynomial order'
292293
write(STD_OUT,'(30X,A,A28,F10.3)') "->", 'Constant for Epsilon: ',self%gauss_epsil
293-
if (self%calculate_with_projection) write(STD_OUT,'(30X,A)') 'Warining, epsilon calculated using properties of element 1'
294+
write(STD_OUT,'(30X,A)') 'Warining, epsilon calculated using properties of element 1'
294295
case default
295296
write(STD_OUT,'(30X,A,A28,ES10.3)') "->", 'Fixed Epsilon value: ',self%gauss_epsil
296297
end select
@@ -413,7 +414,7 @@ subroutine ConstructFarm(self, controlVariables, t0, mesh)
413414
element_loop:do eID = 1, mesh%no_of_elements
414415
do k=1, self%num_turbines
415416
tolerance = self%tolerance_factor*self%turbine_t(k)%radius
416-
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)
417+
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))
417418
if( r_square <= POW2(self%turbine_t(k)%radius+tolerance) &
418419
.and. minval(mesh%elements(eID)%geom%x(1,:,:,:)) < self%turbine_t(k)%hub_cood_x+tolerance &
419420
.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)
428429
element_loop2:do eID = 1, mesh%no_of_elements
429430
do k=1, self%num_turbines
430431
tolerance = self%tolerance_factor*self%turbine_t(k)%radius
431-
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)
432+
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))
432433
if( r_square <= POW2(self%turbine_t(k)%radius+tolerance) &
433434
.and. minval(mesh%elements(eID)%geom%x(1,:,:,:)) < self%turbine_t(k)%hub_cood_x+tolerance &
434435
.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)
440441
end do
441442
end do element_loop2
442443

443-
! precalculate epsilon for projection mode
444-
if (self % calculate_with_projection) then
445-
if (MPI_Process % doMPIAction) then
446-
if (nelem .gt. 0) then
447-
delta_temp = (mesh % elements(elementsActuated(1)) % geom % Volume / product(mesh % elements(elementsActuated(1)) % Nxyz + 1)) ** (1.0_RP / 3.0_RP)
448-
delta_count = 1
449-
else
450-
delta_temp = 0.0_RP
451-
delta_count = 0
452-
end if
444+
select case (self % epsilon_type)
445+
case (0)
446+
! EPSILON - option 1 (from file)
447+
do k = 1, self%num_turbines
448+
self % turbine_t(k) % gauss_epsil = self % gauss_epsil
449+
enddo
450+
case (2)
451+
! EPSILON - option 3 (k is from file)
452+
! eps = k*delta; k is in gauss_epsil of farm
453+
! precalculate delta, for now only using element 1
454+
if (MPI_Process % doMPIAction) then
455+
if (nelem .gt. 0) then
456+
delta_temp = (mesh % elements(elementsActuated(1)) % geom % Volume / product(mesh % elements(elementsActuated(1)) % Nxyz + 1)) ** (1.0_RP / 3.0_RP)
457+
delta_count = 1
458+
else
459+
delta_temp = 0.0_RP
460+
delta_count = 0
461+
end if
453462
#ifdef _HAS_MPI_
454-
call mpi_allreduce(delta_temp, delta, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD, ierr)
455-
call mpi_allreduce(delta_count, delta_paritions, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD, ierr)
463+
call mpi_allreduce(delta_temp, delta, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD, ierr)
464+
call mpi_allreduce(delta_count, delta_paritions, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD, ierr)
456465
#endif
457-
delta = delta / real(delta_paritions,kind=RP)
458-
else
459-
delta = (mesh % elements(elementsActuated(1)) % geom % Volume / product(mesh % elements(elementsActuated(1)) % Nxyz + 1)) ** (1.0_RP / 3.0_RP)
460-
end if
461-
!$omp do schedule(runtime)private(i,j,k)
462-
do k = 1, self%num_turbines
463-
do j = 1, self%turbine_t(k)%num_blades
464-
do i = 1, self%turbine_t(k)%num_blade_sections
465-
self % turbine_t(k) % blade_t(j) % gauss_epsil_delta(i) = delta
466-
end do
467-
enddo
468-
enddo
469-
!$omp end do
470-
end if
466+
delta = delta / real(delta_paritions,kind=RP)
467+
else
468+
delta = (mesh % elements(elementsActuated(1)) % geom % Volume / product(mesh % elements(elementsActuated(1)) % Nxyz + 1)) ** (1.0_RP / 3.0_RP)
469+
end if
470+
do k = 1, self%num_turbines
471+
! e = k*delta
472+
self % turbine_t(k) % gauss_epsil = self % gauss_epsil * delta
473+
enddo
474+
475+
case default
476+
! EPSILON - option 2 using Cd not implementd
477+
do k = 1, self%num_turbines
478+
self % turbine_t(k) % gauss_epsil = self % gauss_epsil
479+
enddo
471480

481+
end select
472482
!
473483
! Create output files
474484
! -------------------
@@ -556,7 +566,7 @@ subroutine UpdateFarm(self,time, mesh)
556566
!local variables
557567
integer :: ii, jj, i, j, k, kk
558568
real(kind=RP) :: dt, interp, delta_temp
559-
logical :: found
569+
logical :: found, allfound
560570
integer :: eID, ierr
561571
real(kind=RP), dimension(NDIM) :: x, xi
562572
real(kind=RP), dimension(NCONS) :: Q, Qtemp
@@ -613,7 +623,7 @@ subroutine UpdateFarm(self,time, mesh)
613623
! use the local Q based on the position of the actuator line point
614624
! ----------------------------------------------------------------
615625
!
616-
!$omp do schedule(runtime)private(ii,jj,kk,eID,Q,Qtemp,delta_temp,xi,found)
626+
!$omp do schedule(runtime)private(ii,jj,kk,eID,Q,Qtemp,delta_temp,xi,found,allfound)
617627
do kk = 1, self%num_turbines
618628
do jj = 1, self%turbine_t(kk)%num_blades
619629

@@ -638,22 +648,21 @@ subroutine UpdateFarm(self,time, mesh)
638648
if (found) then
639649
! interpolate state values in the element
640650
Qtemp = interpolateQ(mesh,eID, xi)
641-
delta_temp = (mesh % elements(eID) % geom % Volume / product(mesh % elements(eID) % Nxyz + 1)) ** (1.0_RP / 3.0_RP)
642651
else
643652
Qtemp = 0.0_RP
644-
delta_temp = 0.0_RP
645653
end if
646654
if ( (MPI_Process % doMPIAction) ) then
647655
#ifdef _HAS_MPI_
648656
call mpi_allreduce(Qtemp, Q, NCONS, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD, ierr)
649-
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)
657+
! 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)
658+
! call mpi_allreduce(found, allfound, NCONS, MPI_LOGICAL, MPI_LOR, MPI_COMM_WORLD, ierr)
659+
! found = allfound
650660
#endif
651661
else
652662
Q = Qtemp
653-
self % turbine_t(kk) % blade_t(jj) % gauss_epsil_delta(ii) = delta_temp
654-
655663
end if
656664
if (all(Q .eq. 0.0_RP)) then
665+
! if (.not. found) then
657666
print*, "Actuator line point not found in mesh, x: ", x
658667
call exit(99)
659668
end if
@@ -720,8 +729,8 @@ subroutine ForcesFarm(self, mesh, time, Level)
720729

721730
do jj = 1, self % turbine_t(kk) % num_blades
722731
do ii = 1, self % turbine_t(kk) % num_blade_sections
723-
interp = GaussianInterpolation(self%epsilon_type, mesh % elements(eID) % geom % x(:,i,j,k), self%turbine_t(kk)%blade_t(jj)%point_xyz_loc(ii,:), &
724-
self%turbine_t(kk)%blade_t(jj)%chord(ii), self%gauss_epsil,self%turbine_t(kk)%blade_t(jj)%gauss_epsil_delta(ii))
732+
interp = GaussianInterpolation( mesh % elements(eID) % geom % x(:,i,j,k), self%turbine_t(kk)%blade_t(jj)%point_xyz_loc(ii,:), &
733+
self%turbine_t(kk)%blade_t(jj)%chord(ii), self%turbine_t(kk)%gauss_epsil)
725734
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)
726735

727736
! minus account action-reaction effect, is the force on the fliud
@@ -775,8 +784,8 @@ subroutine ForcesFarm(self, mesh, time, Level)
775784

776785
do jj = 1, self % turbine_t(kk) % num_blades
777786
do ii = 1, self % turbine_t(kk) % num_blade_sections
778-
interp = GaussianInterpolation(self%epsilon_type, mesh % elements(eID) % geom % x(:,i,j,k), self%turbine_t(kk)%blade_t(jj)%point_xyz_loc(ii,:), &
779-
self%turbine_t(kk)%blade_t(jj)%chord(ii), self%gauss_epsil,self%turbine_t(kk)%blade_t(jj)%gauss_epsil_delta(ii))
787+
interp = GaussianInterpolation(mesh % elements(eID) % geom % x(:,i,j,k), self%turbine_t(kk)%blade_t(jj)%point_xyz_loc(ii,:), &
788+
self%turbine_t(kk)%blade_t(jj)%chord(ii), self%turbine_t(kk)%gauss_epsil)
780789

781790
! minus account action-reaction effect, is the force on the fluid
782791
actuator_source(1) = actuator_source(1) - self%turbine_t(kk)%blade_t(jj)%local_thrust(ii) * interp
@@ -1156,40 +1165,15 @@ End Subroutine FarmUpdateBladeForces
11561165
!
11571166
!///////////////////////////////////////////////////////////////////////////////////////
11581167
!
1159-
Function GaussianInterpolation(epsilon_type, x, x_point, chord, gauss_epsil, delta, Cd)
1168+
Function GaussianInterpolation(x, x_point, chord, gauss_epsil)
11601169
implicit none
1161-
integer, intent(in) :: epsilon_type
11621170
real(kind=RP), intent(in) :: x(NDIM)
11631171
real(kind=RP), intent(in) :: x_point(NDIM)
11641172
real(kind=RP), intent(in) :: chord
11651173
real(kind=RP), intent(in) :: gauss_epsil
1166-
real(kind=RP), intent(in) :: delta
1167-
real(kind=RP), intent(in), optional :: Cd
11681174
real(kind=RP) :: GaussianInterpolation
11691175

1170-
!local variables
1171-
real(kind=RP) :: epsil
1172-
1173-
select case (epsilon_type)
1174-
case (0)
1175-
! EPSILON - option 1 (from file)
1176-
epsil = gauss_epsil
1177-
case (1)
1178-
! EPSILON - option 2
1179-
if (present(Cd)) then
1180-
epsil = max(chord/4.0_RP,chord*Cd/2.0_RP)
1181-
else
1182-
epsil = gauss_epsil
1183-
end if
1184-
case (2)
1185-
! EPSILON - option 3 (k is from file)
1186-
! eps = k*delta; k is in gauss_epsil, gauss_epsil_delta is obtained in UpdateFarm
1187-
epsil = gauss_epsil * delta
1188-
case default
1189-
epsil = gauss_epsil
1190-
end select
1191-
1192-
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) )
1176+
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) )
11931177

11941178
End Function GaussianInterpolation
11951179
!

0 commit comments

Comments
 (0)