diff --git a/src/forces.f90 b/src/forces.f90
index e861dbb383a9bc25098445e13a243d76c0316bff..8d4bdb7a044a796df4b0330a11a2f993fc3566c3 100644
--- a/src/forces.f90
+++ b/src/forces.f90
@@ -640,9 +640,9 @@ contains
!momentum flux
call crossProduct(angularVelocity,[xm,ym,zm]-position,rotationalComponent)
- uxmid = half*(ux1(i,j,k)+ux1(i+1,j,k)) ! - linearVelocity(1) - rotationalComponent(1)
- uymid = half*(uy1(i,j,k)+uy1(i+1,j,k)) !- linearVelocity(2) - rotationalComponent(2)
- uzmid = half*(uz1(i,j,k)+uz1(i+1,j,k)) !- linearVelocity(3) - rotationalComponent(3)
+ uxmid = half*(ux1(i,j,k)+ux1(i+1,j,k)) - linearVelocity(1) - rotationalComponent(1)
+ uymid = half*(uy1(i,j,k)+uy1(i+1,j,k)) - linearVelocity(2) - rotationalComponent(2)
+ uzmid = half*(uz1(i,j,k)+uz1(i+1,j,k)) - linearVelocity(3) - rotationalComponent(3)
fcvx = fcvx -uxmid*uymid*dx*dz
fcvy = fcvy -uymid*uymid*dx*dz
@@ -699,9 +699,9 @@ contains
!momentum flux
call crossProduct(angularVelocity,[xm,ym,zm]-position,rotationalComponent)
- uxmid = half*(ux1(i,j,k)+ux1(i+1,j,k))! - linearVelocity(1) - rotationalComponent(1)
- uymid = half*(uy1(i,j,k)+uy1(i+1,j,k))! - linearVelocity(2) - rotationalComponent(2)
- uzmid = half*(uz1(i,j,k)+uz1(i+1,j,k))! - linearVelocity(3) - rotationalComponent(3)
+ uxmid = half*(ux1(i,j,k)+ux1(i+1,j,k)) - linearVelocity(1) - rotationalComponent(1)
+ uymid = half*(uy1(i,j,k)+uy1(i+1,j,k)) - linearVelocity(2) - rotationalComponent(2)
+ uzmid = half*(uz1(i,j,k)+uz1(i+1,j,k)) - linearVelocity(3) - rotationalComponent(3)
fcvx = fcvx +uxmid*uymid*dx*dz
fcvy = fcvy +uymid*uymid*dx*dz
@@ -758,9 +758,9 @@ contains
! write(*,*) 'Calculating force at left x boundary', [xm,ym,zm]
!momentum flux
call crossProduct(angularVelocity,[xm,ym,zm]-position,rotationalComponent)
- uxmid = half*(ux2(i,j,k)+ux2(i,j+1,k))! - linearVelocity(1) - rotationalComponent(1)
- uymid = half*(uy2(i,j,k)+uy2(i,j+1,k))! - linearVelocity(2) - rotationalComponent(2)
- uzmid = half*(uz2(i,j,k)+uz2(i,j+1,k))! - linearVelocity(3) - rotationalComponent(3)
+ uxmid = half*(ux2(i,j,k)+ux2(i,j+1,k)) - linearVelocity(1) - rotationalComponent(1)
+ uymid = half*(uy2(i,j,k)+uy2(i,j+1,k)) - linearVelocity(2) - rotationalComponent(2)
+ uzmid = half*(uz2(i,j,k)+uz2(i,j+1,k)) - linearVelocity(3) - rotationalComponent(3)
fcvx = fcvx -uxmid*uxmid*del_y(j)*dz
@@ -816,9 +816,9 @@ contains
!momentum flux
call crossProduct(angularVelocity,[xm,ym,zm]-position,rotationalComponent)
- uxmid = half*(ux2(i,j,k)+ux2(i,j+1,k))! - linearVelocity(1) - rotationalComponent(1)
- uymid = half*(uy2(i,j,k)+uy2(i,j+1,k))! - linearVelocity(2) - rotationalComponent(2)
- uzmid = half*(uz2(i,j,k)+uz2(i,j+1,k))! - linearVelocity(3) - rotationalComponent(3)
+ uxmid = half*(ux2(i,j,k)+ux2(i,j+1,k)) - linearVelocity(1) - rotationalComponent(1)
+ uymid = half*(uy2(i,j,k)+uy2(i,j+1,k)) - linearVelocity(2) - rotationalComponent(2)
+ uzmid = half*(uz2(i,j,k)+uz2(i,j+1,k)) - linearVelocity(3) - rotationalComponent(3)
fcvx = fcvx + uxmid*uxmid*del_y(j)*dz
@@ -880,9 +880,9 @@ contains
!momentum flux
call crossProduct(angularVelocity,[xm,ym,zm]-position,rotationalComponent)
- uxmid = half*(ux1(i,j,k)+ux1(i+1,j,k))! - linearVelocity(1) - rotationalComponent(1)
- uymid = half*(uy1(i,j,k)+uy1(i+1,j,k))! - linearVelocity(2) - rotationalComponent(2)
- uzmid = half*(uz1(i,j,k)+uz1(i+1,j,k))! - linearVelocity(3) - rotationalComponent(3)
+ uxmid = half*(ux1(i,j,k)+ux1(i+1,j,k)) - linearVelocity(1) - rotationalComponent(1)
+ uymid = half*(uy1(i,j,k)+uy1(i+1,j,k)) - linearVelocity(2) - rotationalComponent(2)
+ uzmid = half*(uz1(i,j,k)+uz1(i+1,j,k)) - linearVelocity(3) - rotationalComponent(3)
fcvx= fcvx +uxmid*uzmid*dx*dy
fcvy= fcvy +uymid*uzmid*dx*dy
@@ -941,9 +941,9 @@ contains
call crossProduct(angularVelocity,[xm,ym,zm]-position,rotationalComponent)
! write(*,*) 'Calculating force at right z boundary', [xm,ym,zm]
- uxmid = half*(ux1(i,j,k)+ux1(i+1,j,k))! - linearVelocity(1) - rotationalComponent(1)
- uymid = half*(uy1(i,j,k)+uy1(i+1,j,k))! - linearVelocity(2) - rotationalComponent(2)
- uzmid = half*(uz1(i,j,k)+uz1(i+1,j,k))! - linearVelocity(3) - rotationalComponent(3)
+ uxmid = half*(ux1(i,j,k)+ux1(i+1,j,k)) - linearVelocity(1) - rotationalComponent(1)
+ uymid = half*(uy1(i,j,k)+uy1(i+1,j,k)) - linearVelocity(2) - rotationalComponent(2)
+ uzmid = half*(uz1(i,j,k)+uz1(i+1,j,k)) - linearVelocity(3) - rotationalComponent(3)
fcvx= fcvx -uxmid*uzmid*dx*dy
fcvy= fcvy -uymid*uzmid*dx*dy
@@ -1048,25 +1048,27 @@ contains
! endif
enddo
- do k = 1, xsize(3)
- do j = 1, xsize(2)
- do i = 1, xsize(1)
- ux11(i,j,k)=ux01(i,j,k)
- uy11(i,j,k)=uy01(i,j,k)
- uz11(i,j,k)=uz01(i,j,k)
- ux01(i,j,k)=ux1(i,j,k)
- uy01(i,j,k)=uy1(i,j,k)
- uz01(i,j,k)=uz1(i,j,k)
- enddo
- enddo
- enddo
+ if (torques_flag.eq.0) then !only update velocity fields if torque calculation not called.
+ do k = 1, xsize(3)
+ do j = 1, xsize(2)
+ do i = 1, xsize(1)
+ ux11(i,j,k)=ux01(i,j,k)
+ uy11(i,j,k)=uy01(i,j,k)
+ uz11(i,j,k)=uz01(i,j,k)
+ ux01(i,j,k)=ux1(i,j,k)
+ uy01(i,j,k)=uy1(i,j,k)
+ uz01(i,j,k)=uz1(i,j,k)
+ enddo
+ enddo
+ enddo
+ endif
return
end subroutine force
- subroutine torque(ux1,uy1,uz1,ep1,dra1,dra2,dra3,record_var)
+ subroutine torque_calc(ux1,uy1,uz1,ep1,dra1,dra2,dra3,record_var)
USE param
USE variables
@@ -1133,6 +1135,7 @@ contains
real(mytype), dimension(nz) :: drag1, drag2, drag11, drag22
real(mytype), dimension(nz) :: drag3, drag4, drag33, drag44
real(mytype) :: mom1, mom2, mom3, tp1, tp2, tp3
+ real(mytype) :: radial(3),angular_velocity_result(3)
! write(*,*) 'Inside FORCE'
@@ -1257,16 +1260,17 @@ contains
! relative to the volume, and, therefore, this has a sense
! of a "source".
! fac = (1.5*ux1(i,j,k)-2.0*ux01(i,j,k)+0.5*ux11(i,j,k))*epcv1(i,j,k)
- tsumx = tsumx+crossProduct((fac1-coriolis(1)-centrifugal(1)),[xm,ym,zm]-position)*dx*del_y(j+(xstart(2)-1))*dz/dt !tsumx+fac*dx*dy/dt
+ call crossProduct(([fac1,fac2,fac3]),[xm,ym,zm]-position, angular_velocity_result)
+ tsumx = tsumx+angular_velocity_result(1)*dx*del_y(j+(xstart(2)-1))*dz/dt !tsumx+fac*dx*dy/dt
! tsumx = tsumx+fac1*dx*del_y(j+xstart(2)-1)*dz/dt
!sumx(k) = sumx(k)+dudt1*dx*dy
! fac = (1.5*uy1(i,j,k)-2.0*uy01(i,j,k)+0.5*uy11(i,j,k))*epcv1(i,j,k)
- tsumy = tsumy+(fac2-coriolis(2)-centrifugal(2))*dx*del_y(j+(xstart(2)-1))*dz/dt !tsumy+fac*dx*dy/dt
+ tsumy = tsumy+angular_velocity_result(2)*dx*del_y(j+(xstart(2)-1))*dz/dt !tsumy+fac*dx*dy/dt
! tsumy = tsumy+fac2*dx*del_y(j+xstart(2)-1)*dz/dt
!sumy(k) = sumy(k)+dudt1*dx*dy
- tsumz = tsumz+(fac3-coriolis(3)-centrifugal(3))*dx*del_y(j+(xstart(2)-1))*dz/dt
+ tsumz = tsumz+angular_velocity_result(3)*dx*del_y(j+(xstart(2)-1))*dz/dt
! tsumz = tsumz+fac3*dx*del_y(j+xstart(2)-1)*dz/dt
enddo
enddo
@@ -1802,5 +1806,5 @@ contains
return
- end subroutine torque
+ end subroutine torque_calc
end module forces
diff --git a/src/module_param.f90 b/src/module_param.f90
index 9d16bbe878ea1fbff3cb582848cf6a5e16aa7016..96afe048d215dad0dd6cd7770bd7b00e9700a027 100644
--- a/src/module_param.f90
+++ b/src/module_param.f90
@@ -622,6 +622,6 @@ module ibm_param
real(mytype) :: position_1(3),linearVelocity_1(3),orientation_1(4),angularVelocity_1(4)
real(mytype) :: chord,thickness,omega, tconv2_sign
integer :: inana ! Analytical BC as Input
- integer :: imove, nozdrift, force_csv, bodies_fixed, cube_flag
+ integer :: imove, nozdrift, force_csv, bodies_fixed, cube_flag, torques_flag
end module ibm_param
!############################################################################
diff --git a/src/xcompact3d.f90 b/src/xcompact3d.f90
index b5718b0b468d5d5bec991e7a69bb65c0adef6c4e..9716d7bc702285b2d8e6e29d5d41651ae9bf0387 100644
--- a/src/xcompact3d.f90
+++ b/src/xcompact3d.f90
@@ -17,9 +17,9 @@ program xcompact3d
use ibm, only : body
use genepsi, only : genepsi3d
use ellipsoid_utils, only: lin_step, ang_step, QuaternionNorm
- use forces, only : force, init_forces, iforces,update_forces, xld,xrd,yld,yud,zld,zrd
+ use forces, only : force, init_forces, iforces,update_forces, xld,xrd,yld,yud,zld,zrd,torque_calc
implicit none
- real(mytype) :: dummy,drag,lift,lat,grav_effy,grav_effx,grav_effz
+ real(mytype) :: dummy,drag,lift,lat,grav_effy,grav_effx,grav_effz,xtorq,ytorq,ztorq
integer :: iounit,ierr,i
real, dimension(100) :: x
@@ -139,8 +139,13 @@ program xcompact3d
close(20)
endif
- torque(:)=zero
-
+
+ if (torques_flag.eq.1) then
+ call torque_calc(ux1,uy1,uz1,ep1,xtorq,ytorq,ztorq,1)
+ torque = [xtorq,ytorq,ztorq]
+ else
+ torque(:) = zero
+ endif
! if (nrank==0) then
! if (bodies_fixed==0) then