diff --git a/src/BC-Ellipsoid.f90 b/src/BC-Ellipsoid.f90
index 7da85bea95ce3d2dfdadf15207228795f295a0d2..9df6e2b5b75e8d4d481de2a83191f686951b0898 100644
--- a/src/BC-Ellipsoid.f90
+++ b/src/BC-Ellipsoid.f90
@@ -22,7 +22,7 @@ subroutine geomcomplex_ellip(epsi,nxi,nxf,ny,nyi,nyf,nzi,nzf,dx,yp,remp)
use param, only : one, two, ten
use ibm_param
use dbg_schemes, only: sqrt_prec
- use ellipsoid_utils, only: EllipsoidalRadius, NormalizeQuaternion
+ use ellipsoid_utils, only: NormalizeQuaternion, is_inside_ellipsoid
implicit none
@@ -36,6 +36,7 @@ subroutine geomcomplex_ellip(epsi,nxi,nxf,ny,nyi,nyf,nzi,nzf,dx,yp,remp)
real(mytype) :: zeromach
real(mytype) :: cexx,ceyy,cezz,dist_axi
real(mytype) :: point(3)
+ logical :: is_inside
zeromach=one
do while ((one + zeromach / two) .gt. one)
@@ -54,16 +55,16 @@ subroutine geomcomplex_ellip(epsi,nxi,nxf,ny,nyi,nyf,nzi,nzf,dx,yp,remp)
! Update center of moving ellipsoid
- if (t.ne.0.) then
- cexx=cex+lvx*(t-ifirst*dt)
- ceyy=cey+lvy*(t-ifirst*dt)
- cezz=cez+lvz*(t-ifirst*dt)
- else
- cexx=cex
- ceyy=cey
- cezz=cez
- endif
- position=[cexx,ceyy,cezz]
+ ! if (t.ne.0.) then
+ ! cexx=cex+lvx*(t-ifirst*dt)
+ ! ceyy=cey+lvy*(t-ifirst*dt)
+ ! cezz=cez+lvz*(t-ifirst*dt)
+ ! else
+ ! cexx=cex
+ ! ceyy=cey
+ ! cezz=cez
+ ! endif
+ ! position=[cexx,ceyy,cezz]
! write(*,*) position
! ce=[cexx, ceyy, cezz]
!
@@ -78,11 +79,12 @@ subroutine geomcomplex_ellip(epsi,nxi,nxf,ny,nyi,nyf,nzi,nzf,dx,yp,remp)
do i=nxi,nxf
xm=real(i-1,mytype)*dx
point=[xm, ym, zm]
- call EllipsoidalRadius(point, position, orientation, shape, r)
+ ! call EllipsoidalRadius(point, position, orientation, shape, r)
+ call is_inside_ellipsoid(point, position, orientation, shape, ra, zeromach, is_inside)
! r=sqrt_prec((xm-cexx)**two+(ym-ceyy)**two+(zm-cezz)**two)
! r=sqrt_prec((xm-cexx)**two+(ym-ceyy)**two)
- if (r-ra.gt.zeromach) then
+ if (.not.is_inside) then
! write(*,*) i, j, k
cycle
endif
@@ -230,7 +232,7 @@ subroutine init_ellip (ux1,uy1,uz1,phi1)
USE MPI
USE ibm_param
use dbg_schemes, only: exp_prec
- use ellipsoid_utils, only: NormalizeQuaternion
+ use ellipsoid_utils, only: NormalizeQuaternion,ellipInertiaCalculate
implicit none
@@ -248,13 +250,18 @@ subroutine init_ellip (ux1,uy1,uz1,phi1)
call NormalizeQuaternion(orientation)
position=[cex,cey,cez]
linearVelocity=[lvx,lvy,lvz]
- angularVelocity=[avx,avy,avz]
-
- write(*,*) 'set shape = ', shape
- write(*,*) 'set orientation = ', orientation
- write(*,*) 'set position = ', position
- write(*,*) 'set linear velocity = ', linearVelocity
- write(*,*) 'set angular velocity = ', angularVelocity
+ angularVelocity=[zero,avx,avy,avz]
+ call ellipInertiaCalculate(shape,rho_s,inertia)
+
+ if (nrank==0) then
+ write(*,*) 'set shape = ', shape
+ write(*,*) 'set orientation = ', orientation
+ write(*,*) 'set position = ', position
+ write(*,*) 'set linear velocity = ', linearVelocity
+ write(*,*) 'set angular velocity = ', angularVelocity
+ write(*,*) 'set moment of inertia = ', inertia
+ write(*,*) 'density of solid = ', rho_s
+ end if
if (iscalar==1) then
diff --git a/src/ellip_utils.f90 b/src/ellip_utils.f90
index 32ed16ee217d9c78cf767bc540a99b6b42262fff..1e5ab49b00ece7ae3121f7e2e1223ab5ec001271 100644
--- a/src/ellip_utils.f90
+++ b/src/ellip_utils.f90
@@ -242,6 +242,18 @@ contains
pointVelocity = crossed + linearVelocity
end subroutine CalculatePointVelocity
+ subroutine is_inside_ellipsoid(point, centre, orientation, shape, ra, zeromach, is_inside)
+ real(mytype), intent(in) :: point(3), centre(3), orientation(4), shape(3), ra, zeromach
+ logical,intent(out) :: is_inside
+ real(mytype) :: r
+
+ call EllipsoidalRadius(point,centre,orientation,shape,r)
+
+ is_inside = ((r-ra).lt.zeromach)
+
+ end subroutine is_inside_ellipsoid
+
+
subroutine navierFieldGen(center, linearVelocity, angularVelocity, ep1, ep1_x, ep1_y, ep1_z)
use param
use decomp_2d
@@ -379,6 +391,7 @@ contains
real(mytype),intent(out) :: ang_accel(4)
real(mytype) :: inertia_inv(3,3),omega_v(3),torque_v(3),test(3),crossed(3)
+ ! write(*,*) 'inverting ', inertia
call invert_3x3_matrix(inertia,inertia_inv)
omega_v=omega(2:4)
torque_v=torque_b(2:4)
@@ -390,9 +403,9 @@ contains
- subroutine ang_half_step(q, omega_q, torque, q_new, o_new)
+ subroutine ang_half_step(q, omega_q, torque_vec, q_new, o_new)
use param
- real(mytype),intent(in) :: q(4),omega_q(4),torque(3)
+ real(mytype),intent(in) :: q(4),omega_q(4),torque_vec(3)
real(mytype),intent(out) :: q_new(4),o_new(4)
real(mytype) :: inertia(3,3)
real(mytype) :: omega_b(4),torque_q(4),ang_accel_b(4),torque_b(4)
@@ -401,7 +414,7 @@ contains
call lab_to_body(omega_q,q,omega_b)
torque_q(1)=zero
- torque_q(2:4)=torque
+ torque_q(2:4)=torque_vec(:)
call lab_to_body(torque_q,q,torque_b)
@@ -420,9 +433,9 @@ contains
end subroutine ang_half_step
- subroutine ang_full_step(q,omega_q,q_half,omega_n_half,torque,q_full,omega_full)
+ subroutine ang_full_step(q,omega_q,q_half,omega_n_half,torque_vec,q_full,omega_full)
use param
- real(mytype),intent(in) :: q(4),omega_q(4),q_half(4),omega_n_half(4),torque(3)
+ real(mytype),intent(in) :: q(4),omega_q(4),q_half(4),omega_n_half(4),torque_vec(3)
real(mytype),intent(out) :: q_full(4),omega_full(4)
real(mytype) :: inertia(3,3)
real(mytype) :: omega_b(4),omega_n_half_b(4),omega_full_b(4)
@@ -433,7 +446,7 @@ contains
call lab_to_body(omega_n_half,q_half,omega_n_half_b)
torque_q(1)=zero
- torque_q(2:4)=torque
+ torque_q(2:4)=torque_vec(:)
call lab_to_body(torque_q,q_half,torque_b)
call accel_get(omega_n_half_b,inertia,torque_b,ang_accel_half_b)
@@ -446,6 +459,41 @@ contains
end subroutine ang_full_step
+ subroutine ang_step(q,omega_q,torque_vec,time_step,q1,omega1)
+ use param
+ use ibm_param
+ real(mytype),intent(in) :: q(4),omega_q(4),torque_vec(3),time_step
+ real(mytype),intent(out):: q1(4),omega1(4)
+ real(mytype) :: torque_q(4),omega_b(4),torque_b(4),omega_half_b(4),omega1_b(4)
+ real(mytype) :: ang_accel_b(4), omega_half(4)
+
+ ! write(*,*) 'ang_vel_lab = ', omega_q
+ call lab_to_body(omega_q, q, omega_b) !convert to body frame
+ ! write(*,*) 'ang_vel_b = ', omega_b
+ torque_q(1)=zero
+ torque_q(2:4)=torque_vec(:)
+ ! write(*,*) 'torque = ', torque_q
+ call lab_to_body(torque_q,q, torque_b)
+ ! write(*,*) 'torque_b = ', torque_b
+
+ call accel_get(omega_b,inertia,torque_b,ang_accel_b) !calculate acceleration
+ ! write(*,*) 'acceleration =', ang_accel_b
+
+ call omega_stepper(omega_b, ang_accel_b,time_step*half,omega_half_b)
+ ! write(*,*) 'omega_half_b =', omega_half_b
+ call omega_stepper(omega_b, ang_accel_b,time_step,omega1_b) !calculate omega at half and full timestep
+ ! write(*,*) 'omega_full_b =', omega1_b
+ call body_to_lab(omega_half_b,q,omega_half) !convert back to lab
+ ! write(*,*) 'omega_half_lab', omega_half
+ call orientation_stepper(q,omega_half,time_step,q1) !step forward orientation
+ ! write(*,*) 'time_step =', time_step
+ ! write(*,*) 'orientation1 =', q1
+ call body_to_lab(omega1_b,q1,omega1)
+ ! write(*,*) 'omega_full =', omega1
+
+ end subroutine ang_step
+
+
subroutine lin_step(position,linearVelocity,linearAcceleration,time_step,position_1,linearVelocity_1)
real(mytype),intent(in) :: position(3),linearVelocity(3),linearAcceleration(3),time_step
real(mytype),intent(out) :: position_1(3),linearVelocity_1(3)
@@ -454,4 +502,41 @@ contains
linearVelocity_1 = linearVelocity(:) + time_step*linearAcceleration(:)
end subroutine lin_step
+
+ subroutine ellipMassCalculate(shape,rho_s,mass)
+ use constants, only: pi
+ real(mytype),intent(in) :: shape(3),rho_s
+ real(mytype),intent(out) :: mass
+ real(mytype) :: a,b,c,vol
+
+ a=shape(1)
+ b=shape(2)
+ c=shape(3)
+ vol=(4_mytype/3_mytype)*pi*a*b*c
+ mass=vol*rho_s
+
+ end subroutine ellipMassCalculate
+
+ subroutine ellipInertiaCalculate(shape,rho_s,inertia)
+ real(mytype),intent(in) :: shape(3),rho_s
+ real(mytype),intent(out) :: inertia(3,3)
+ real(mytype) :: a,b,c,i1,i2,i3,mass
+
+ call ellipMassCalculate(shape,rho_s,mass)
+
+ a=shape(1)
+ b=shape(2)
+ c=shape(3)
+
+ i1=mass*(b**2+c**2)*0.2
+ i2=mass*(a**2+c**2)*0.2
+ i3=mass*(a**2+b**2)*0.2
+
+ inertia(:,:)=0_mytype
+ inertia(1,1)=i1
+ inertia(2,2)=i2
+ inertia(3,3)=i3
+
+ end subroutine ellipInertiaCalculate
+
end module ellipsoid_utils
diff --git a/src/module_param.f90 b/src/module_param.f90
index 7a0779692432f45a575c9936168d5a6c3094a9dd..e375dddb011a24d9d2ea816c20419fb0ee1b6fe0 100644
--- a/src/module_param.f90
+++ b/src/module_param.f90
@@ -617,9 +617,9 @@ end module simulation_stats
!############################################################################
module ibm_param
use decomp_2d, only : mytype
- real(mytype) :: cex,cey,cez,shx,shy,shz,oriw,orii,orij,orik,lvx,lvy,lvz,avx,avy,avz,ra,ubcx,ubcy,ubcz,rads, c_air
- real(mytype) :: position(3),orientation(4),linearVelocity(3),angularVelocity(3),linearAcceleration(3),shape(3)
- real(mytype) :: position_1(3),linearVelocity_1(3)
+ real(mytype) :: cex,cey,cez,shx,shy,shz,oriw,orii,orij,orik,lvx,lvy,lvz,avx,avy,avz,ra,ubcx,ubcy,ubcz,rads,rho_s,c_air
+ real(mytype) :: position(3),orientation(4),linearVelocity(3),angularVelocity(4),linearAcceleration(3),torque(3),shape(3),inertia(3,3)
+ real(mytype) :: position_1(3),linearVelocity_1(3),orientation_1(4),angularVelocity_1(4)
real(mytype) :: chord,thickness,omega
integer :: inana ! Analytical BC as Input
integer :: imove
diff --git a/src/parameters.f90 b/src/parameters.f90
index 6ed4b579d052cc1efc513860d36212702b5972b6..3dbc63d09701b4e644c4cebd4ac7dfb8d80bec95 100644
--- a/src/parameters.f90
+++ b/src/parameters.f90
@@ -60,7 +60,7 @@ subroutine parameter(input_i3d)
NAMELIST /LESModel/ jles, smagcst, smagwalldamp, nSmag, walecst, maxdsmagcst, iwall
NAMELIST /WallModel/ smagwalldamp
NAMELIST /Tripping/ itrip,A_tr,xs_tr_tbl,ys_tr_tbl,ts_tr_tbl,x0_tr_tbl
- NAMELIST /ibmstuff/ cex,cey,cez,shx,shy,shz,oriw,orii,orij,orik,lvx,lvy,lvz,avx,avy,avz,ra,nobjmax,nraf,nvol,iforces, npif, izap, ianal, imove, thickness, chord, omega ,ubcx,ubcy,ubcz,rads, c_air
+ NAMELIST /ibmstuff/ cex,cey,cez,shx,shy,shz,oriw,orii,orij,orik,lvx,lvy,lvz,avx,avy,avz,ra,nobjmax,nraf,nvol,iforces, npif, izap, ianal, imove, thickness, chord, omega ,ubcx,ubcy,ubcz,rads,rho_s, c_air
NAMELIST /ForceCVs/ xld, xrd, yld, yud!, zld, zrd
NAMELIST /LMN/ dens1, dens2, prandtl, ilmn_bound, ivarcoeff, ilmn_solve_temp, &
massfrac, mol_weight, imultispecies, primary_species, &
@@ -571,6 +571,7 @@ subroutine parameter_defaults()
use probes, only : nprobes, flag_all_digits, flag_extra_probes
use visu, only : output2D
use forces, only : iforces, nvol
+ use ibm_param, only : oriw,rho_s
implicit none
@@ -610,6 +611,8 @@ subroutine parameter_defaults()
wrotation = zero
irotation = 0
itest=1
+ oriw=one
+ rho_s=one
!! Gravity field
gravx = zero
diff --git a/src/tools.f90 b/src/tools.f90
index 13f352331462dd0594e2e9a82138e4a298f0e131..8088d224bd0cbc92f4ae9fe8d0a032ff991b25e6 100644
--- a/src/tools.f90
+++ b/src/tools.f90
@@ -143,6 +143,7 @@ contains
use simulation_stats
use var
use MPI
+ use ibm_param, only: position,orientation
implicit none
@@ -160,6 +161,8 @@ contains
call cpu_time(time1)
write(*,*) '==========================================================='
write(*,"(' Time step =',i7,'/',i7,', Time unit =',F9.4)") itime,ilast,t
+ write(*,*) 'Centroid position = ',real(position(1),4),real(position(2),4),real(position(3),4)
+ write(*,*) 'Orientation = ',real(orientation(1),4),real(orientation(2),4),real(orientation(3),4),real(orientation(4),4)
endif
else if ((iwhen == 3).and.(itime > ifirst)) then !AT THE END OF A TIME STEP
if (nrank == 0.and.(mod(itime, ilist) == 0 .or. itime == ifirst .or. itime==ilast)) then
diff --git a/src/xcompact3d.f90 b/src/xcompact3d.f90
index 3af729dd3bf2a302e89bf89c767a45ecbc7f93db..6e97086bc9a03bba0e95db6d712fd86dcb6e0ae1 100644
--- a/src/xcompact3d.f90
+++ b/src/xcompact3d.f90
@@ -16,7 +16,7 @@ program xcompact3d
use ibm_param
use ibm, only : body
use genepsi, only : genepsi3d
- use ellipsoid_utils, only: lin_step
+ use ellipsoid_utils, only: lin_step, ang_step
implicit none
@@ -76,11 +76,25 @@ program xcompact3d
call test_flow(rho1,ux1,uy1,uz1,phi1,ep1,drho1,divu3)
+ !Add force calculation here
+
+ linearAcceleration(:)=zero
+ torque(:)=zero
+
call lin_step(position,linearVelocity,linearAcceleration,dt,position_1,linearVelocity_1)
+ call ang_step(orientation,angularVelocity,torque,dt,orientation_1,angularVelocity_1)
position = position_1
linearVelocity = linearVelocity_1
+ orientation = orientation_1
+ angularVelocity = angularVelocity_1
+
+ ! if (nrank==0) then
+ ! write(*,*) 'Centroid position is ', position
+ ! write(*,*) 'Orientation is ', orientation
+ ! end if
+
enddo !! End sub timesteps
call restart(ux1,uy1,uz1,dux1,duy1,duz1,ep1,pp3(:,:,:,1),phi1,dphi1,px1,py1,pz1,rho1,drho1,mu1,1)