From bafdf492c8fc82107a7a29d47afaf44824d50b31 Mon Sep 17 00:00:00 2001
From: s2006749 <s2006749@ed.ac.uk>
Date: Wed, 9 Oct 2024 16:25:35 +0100
Subject: [PATCH] vectorising all global body parameters to handle up to 10
bodies
---
src/BC-Ellipsoid.f90 | 75 ++++++++++++------------
src/ellip_utils.f90 | 17 ++++++
src/forces.f90 | 132 +++++++++++++++++++++----------------------
src/genepsi3d.f90 | 70 +++++++++++++++++++----
src/module_param.f90 | 9 +--
src/parameters.f90 | 8 ++-
src/xcompact3d.f90 | 63 ++++++++++++---------
7 files changed, 228 insertions(+), 146 deletions(-)
diff --git a/src/BC-Ellipsoid.f90 b/src/BC-Ellipsoid.f90
index c54e03a4..94fc9562 100644
--- a/src/BC-Ellipsoid.f90
+++ b/src/BC-Ellipsoid.f90
@@ -31,7 +31,7 @@ subroutine geomcomplex_ellip(epsi,nxi,nxf,ny,nyi,nyf,nzi,nzf,dx,yp,dz,remp)
real(mytype),dimension(ny) :: yp
real(mytype) :: dx,dz
real(mytype) :: remp
- integer :: i,j,k
+ integer :: i,j,k, i_body
real(mytype) :: xm,ym,zm,r,rads2,kcon
real(mytype) :: zeromach
real(mytype) :: cexx,ceyy,cezz,dist_axi
@@ -81,27 +81,30 @@ subroutine geomcomplex_ellip(epsi,nxi,nxf,ny,nyi,nyf,nzi,nzf,dx,yp,dz,remp)
xm=real(i-1,mytype)*dx
point=[xm, ym, zm]
! call EllipsoidalRadius(point, position, orientation, shape, r)
- if (cube_flag.eq.0) then
- call EllipsoidalRadius(point,position,orientation,shape,r)
- is_inside = (r-ra).lt.zeromach
-
- if (ra /= ra) then
- write(*,*) "Nrank = ", nrank
- write(*,*) "Point = ", point
+ do i_body = 1,nbody
+ if (cube_flag.eq.0) then
+ call EllipsoidalRadius(point,position(i_body,:),orientation(i_body,:),shape(i_body,:),r)
+ is_inside = (r-ra(i_body)).lt.zeromach
+
+ if (ra(i_body) /= ra(i_body)) then
+ write(*,*) "Nrank = ", nrank
+ write(*,*) "Point = ", point
+ endif
+ else if (cube_flag.eq.1) then
+ is_inside = (abs(xm-position(i_body,1)).lt.ra(i_body)).and.(abs(ym-position(i_body,2)).lt.ra(i_body)).and.(abs(zm-position(i_body,3)).lt.ra(i_body))
endif
- else if (cube_flag.eq.1) then
- is_inside = (abs(xm-position(1)).lt.ra).and.(abs(ym-position(2)).lt.ra).and.(abs(zm-position(3)).lt.ra)
- endif
- ! r=sqrt_prec((xm-cexx)**two+(ym-ceyy)**two+(zm-cezz)**two)
- ! r=sqrt_prec((xm-cexx)**two+(ym-ceyy)**two)
- if (.not.is_inside) then
- ! write(*,*) i, j, k
- cycle
- endif
+ ! r=sqrt_prec((xm-cexx)**two+(ym-ceyy)**two+(zm-cezz)**two)
+ ! r=sqrt_prec((xm-cexx)**two+(ym-ceyy)**two)
+ if (is_inside) then
+ ! write(*,*) i, j, k
+ epsi(i,j,k)=remp
+ cycle
+ endif
+ enddo
! write(*,*) is_inside
! write(*,*) i, j, k, zm
- epsi(i,j,k)=remp
+ ! epsi(i,j,k)=remp
! write(*,*) remp
enddo
enddo
@@ -290,26 +293,26 @@ subroutine init_ellip (ux1,uy1,uz1,phi1)
! write(*,*) 'INSIDE INIT ELLIP'
- eqr=(shx*shy*shz)**(1.0/3.0)
- shape=[shx/eqr,shy/eqr,shz/eqr]
+ ! eqr=(sh(1)*sh(2)*sh(3))**(1.0/3.0)
+ ! shape=sh(:)/eqr
- orientation=[oriw,orii,orij,orik]
- call NormalizeQuaternion(orientation)
- position=[cex,cey,cez]
- linearVelocity=[lvx,lvy,lvz]
- angularVelocity=[zero,avx,avy,avz]
- call ellipInertiaCalculate(shape,rho_s,inertia)
- call ellipMassCalculate(shape,rho_s,ellip_m)
+ ! orientation=ori
+ ! call NormalizeQuaternion(orientation)
+ ! position=ce
+ ! linearVelocity=lv
+ ! angularVelocity=[zero, av(1), av(2), av(3)]
+ ! call ellipInertiaCalculate(shape,rho_s,inertia)
+ ! call ellipMassCalculate(shape,rho_s,ellip_m)
- 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 (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 db683465..ee2cf0b5 100644
--- a/src/ellip_utils.f90
+++ b/src/ellip_utils.f90
@@ -7,6 +7,7 @@ module ellipsoid_utils
use decomp_2d, only: mytype
use param, only: zero, one, two
use dbg_schemes, only: sqrt_prec, cos_prec, exp_prec, sin_prec
+ use ibm_param
implicit none
! public QuatRot, cross, IsoKernel, AnIsoKernel, int2str
@@ -257,6 +258,22 @@ contains
pointVelocity = crossed + linearVelocity
end subroutine CalculatePointVelocity
+ subroutine CalculatePointVelocity_Multi(point, pointVelocity)
+ real(mytype), intent(in) :: point(3)
+ real(mytype), intent(out):: pointVelocity(3)
+ real(mytype) :: radii(10),r
+ integer :: i,i_closest
+
+ radii(:) = 10000000.
+ do i = 1,nbody
+ call EllipsoidalRadius(point, position(i,:), orientation(i,:), shape(i,:), r)
+ radii(i) = r
+ enddo
+ i_closest = minloc(radii)
+
+ call CalculatePointVelocity(point, position(i_closest,:), angularVelocity(i_closest,:), linearVelocity(i_closest, :), pointVelocity)
+ end subroutine
+
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
diff --git a/src/forces.f90 b/src/forces.f90
index 74f60675..044673aa 100644
--- a/src/forces.f90
+++ b/src/forces.f90
@@ -388,7 +388,7 @@ contains
real(mytype), dimension(xsize(1),xsize(2),xsize(3)),intent(in) :: ux1, uy1, uz1
real(mytype), dimension(xsize(1),xsize(2),xsize(3)),intent(in) :: ep1
integer, intent(in) ::record_var
- real(mytype), intent(out) :: dra1,dra2,dra3
+ real(mytype), intent(out) :: dra1(10),dra2(10),dra3(10)
real(mytype), dimension(ysize(1),ysize(2),ysize(3)) :: ppi2
real(mytype), dimension(zsize(1),zsize(2),zsize(3)) :: ppi3
@@ -552,8 +552,8 @@ contains
fac2 = (onepfive*uy1(i,j,k)-two*uy01(i,j,k)+half*uy11(i,j,k))*(one-ep1(i,j,k))
fac3 = (onepfive*uz1(i,j,k)-two*uz01(i,j,k)+half*uz11(i,j,k))*(one-ep1(i,j,k))
- call coriolis_force(angularVelocity(2:4),[fac1,fac2,fac3],coriolis)
- call centrifugal_force(angularVelocity(2:4), [xm,ym,zm]-position,centrifugal)
+ call coriolis_force(angularVelocity(iv,2:4),[fac1,fac2,fac3],coriolis)
+ call centrifugal_force(angularVelocity(iv,2:4), [xm,ym,zm]-position(iv,:),centrifugal)
! The velocity time rate has to be relative to the cell center,
! and not to the nodes, because, here, we have an integral
! relative to the volume, and, therefore, this has a sense
@@ -642,10 +642,10 @@ contains
! write(*,*) 'Calculating force at upper y boundary', [xm,ym,zm]
!momentum flux
- call crossProduct(angularVelocity(2:4),[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)
+ call crossProduct(angularVelocity(iv,2:4),[xm,ym,zm]-position(iv,:),rotationalComponent)
+ uxmid = half*(ux1(i,j,k)+ux1(i+1,j,k)) - linearVelocity(iv,1) + rotationalComponent(1)
+ uymid = half*(uy1(i,j,k)+uy1(i+1,j,k)) - linearVelocity(iv,2) + rotationalComponent(2)
+ uzmid = half*(uz1(i,j,k)+uz1(i+1,j,k)) - linearVelocity(iv,3) + rotationalComponent(3)
fcvx = fcvx -uxmid*uymid*dx*dz
fcvy = fcvy -uymid*uymid*dx*dz
@@ -701,10 +701,10 @@ contains
! write(*,*) 'Calculating force at lower y boundary', [xm,ym,zm]
!momentum flux
- call crossProduct(angularVelocity(2:4),[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)
+ call crossProduct(angularVelocity(iv,2:4),[xm,ym,zm]-position(iv,:),rotationalComponent)
+ uxmid = half*(ux1(i,j,k)+ux1(i+1,j,k)) - linearVelocity(iv,1) + rotationalComponent(1)
+ uymid = half*(uy1(i,j,k)+uy1(i+1,j,k)) - linearVelocity(iv,2) + rotationalComponent(2)
+ uzmid = half*(uz1(i,j,k)+uz1(i+1,j,k)) - linearVelocity(iv,3) + rotationalComponent(3)
fcvx = fcvx +uxmid*uymid*dx*dz
fcvy = fcvy +uymid*uymid*dx*dz
@@ -760,10 +760,10 @@ contains
ym=real(jj,mytype)*dz
! write(*,*) 'Calculating force at left x boundary', [xm,ym,zm]
!momentum flux
- call crossProduct(angularVelocity(2:4),[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)
+ call crossProduct(angularVelocity(iv,2:4),[xm,ym,zm]-position(iv,:),rotationalComponent)
+ uxmid = half*(ux2(i,j,k)+ux2(i,j+1,k)) - linearVelocity(iv,1) + rotationalComponent(1)
+ uymid = half*(uy2(i,j,k)+uy2(i,j+1,k)) - linearVelocity(iv,2) + rotationalComponent(2)
+ uzmid = half*(uz2(i,j,k)+uz2(i,j+1,k)) - linearVelocity(iv,3) + rotationalComponent(3)
fcvx = fcvx -uxmid*uxmid*del_y(j)*dz
@@ -818,10 +818,10 @@ contains
! write(*,*) 'Calculating force at right x boundary', [xm,ym,zm]
!momentum flux
- call crossProduct(angularVelocity(2:4),[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)
+ call crossProduct(angularVelocity(iv,2:4),[xm,ym,zm]-position(iv,:),rotationalComponent)
+ uxmid = half*(ux2(i,j,k)+ux2(i,j+1,k)) - linearVelocity(iv,1) + rotationalComponent(1)
+ uymid = half*(uy2(i,j,k)+uy2(i,j+1,k)) - linearVelocity(iv,2) + rotationalComponent(2)
+ uzmid = half*(uz2(i,j,k)+uz2(i,j+1,k)) - linearVelocity(iv,3) + rotationalComponent(3)
fcvx = fcvx + uxmid*uxmid*del_y(j)*dz
@@ -882,10 +882,10 @@ contains
! write(*,*) 'Calculating force at left z boundary', [xm,ym,zm]
!momentum flux
- call crossProduct(angularVelocity(2:4),[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)
+ call crossProduct(angularVelocity(iv,2:4),[xm,ym,zm]-position(iv,:),rotationalComponent)
+ uxmid = half*(ux1(i,j,k)+ux1(i+1,j,k)) - linearVelocity(iv,1) + rotationalComponent(1)
+ uymid = half*(uy1(i,j,k)+uy1(i+1,j,k)) - linearVelocity(iv,2) + rotationalComponent(2)
+ uzmid = half*(uz1(i,j,k)+uz1(i+1,j,k)) - linearVelocity(iv,3) + rotationalComponent(3)
fcvx= fcvx +uxmid*uzmid*dx*dy
fcvy= fcvy +uymid*uzmid*dx*dy
@@ -941,12 +941,12 @@ contains
ii=xstart(1)+i-1
xm=real(ii,mytype)*dx
!momentum flux
- call crossProduct(angularVelocity(2:4),[xm,ym,zm]-position,rotationalComponent)
+ call crossProduct(angularVelocity(iv,2:4),[xm,ym,zm]-position(iv,:),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(iv,1) + rotationalComponent(1)
+ uymid = half*(uy1(i,j,k)+uy1(i+1,j,k)) - linearVelocity(iv,2) + rotationalComponent(2)
+ uzmid = half*(uz1(i,j,k)+uz1(i+1,j,k)) - linearVelocity(iv,3) + rotationalComponent(3)
fcvx= fcvx -uxmid*uzmid*dx*dy
fcvy= fcvy -uymid*uzmid*dx*dy
@@ -1006,9 +1006,9 @@ contains
mom2 = sum(tunsty(:)) + sum(tconvy(:)) + (-2.0*tconv2_sign+1.0)*sum(tconvy2(:))
mom3 = sum(tunstz(:)) + sum(tconvz(:)) + (-2.0*tconv2_sign+1.0)*sum(tconvz2(:))
- dra1 = (sum(tdiffx) + sum(tdiffx2) + tp1 - mom1)
- dra2 = (sum(tdiffy) + sum(tdiffy2) + tp2 - mom2)
- dra3 = (sum(tdiffz) + sum(tdiffz2) + tp3 - mom3)
+ dra1(iv) = (sum(tdiffx) + sum(tdiffx2) + tp1 - mom1)
+ dra2(iv) = (sum(tdiffy) + sum(tdiffy2) + tp2 - mom2)
+ dra3(iv) = (sum(tdiffz) + sum(tdiffz2) + tp3 - mom3)
do k=1,zsize(3)
@@ -1096,7 +1096,7 @@ contains
real(mytype), dimension(xsize(1),xsize(2),xsize(3)),intent(in) :: ux1, uy1, uz1
real(mytype), dimension(xsize(1),xsize(2),xsize(3)),intent(in) :: ep1
integer, intent(in) ::record_var
- real(mytype), intent(out) :: dra1,dra2,dra3
+ real(mytype), intent(out) :: dra1(10),dra2(10),dra3(10)
real(mytype), dimension(ysize(1),ysize(2),ysize(3)) :: ppi2
real(mytype), dimension(zsize(1),zsize(2),zsize(3)) :: ppi3
@@ -1165,7 +1165,7 @@ contains
open(45+(iv-1),file=filename,status='unknown',form='formatted')
! write(*,*) 'Opened file: ', filename, 'number = ', 38+(iv-1)
endif
- if ((nrank.eq.0).and.(torq_debug.eq.1)) then
+ if ((nrank.eq.0).and.(torq_debug.eq.1).and.(iv.eq.1)) then
open(100,file="backside.dat",status='unknown',form='formatted')
open(101,file="frontside.dat",status='unknown',form='formatted')
open(102,file="lowerside.dat",status='unknown',form='formatted')
@@ -1276,7 +1276,7 @@ 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)
- call crossProduct(([fac1,fac2,fac3]),[xm,ym,zm]-position, angular_velocity_result)
+ call crossProduct(([fac1,fac2,fac3]),[xm,ym,zm]-position(iv,:), 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
@@ -1360,13 +1360,13 @@ contains
ii=xstart(1)+i-1
xm=real(ii,mytype)*dx
! write(*,*) 'Calculating force at upper y boundary', [xm,ym,zm]
- radial = [xm,ym,zm]-position
+ radial = [xm,ym,zm]-position(iv,:)
!momentum flux
- call crossProduct(angularVelocity(2:4),[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)
+ call crossProduct(angularVelocity(iv,2:4),[xm,ym,zm]-position(iv,:),rotationalComponent)
+ uxmid = half*(ux1(i,j,k)+ux1(i+1,j,k)) - linearVelocity(iv,1) + rotationalComponent(1)
+ uymid = half*(uy1(i,j,k)+uy1(i+1,j,k)) - linearVelocity(iv,2) + rotationalComponent(2)
+ uzmid = half*(uz1(i,j,k)+uz1(i+1,j,k)) - linearVelocity(iv,3) + rotationalComponent(3)
fcvx = fcvx -(uymid*radial(3)-uzmid*radial(2))*uymid*dx*dz
fcvy = fcvy -(uzmid*radial(1)-uxmid*radial(3))*uymid*dx*dz
@@ -1435,11 +1435,11 @@ contains
! write(*,*) 'Calculating force at lower y boundary', [xm,ym,zm]
!momentum flux
- call crossProduct(angularVelocity(2:4),[xm,ym,zm]-position,rotationalComponent)
- radial = [xm,ym,zm]-position
- 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)
+ call crossProduct(angularVelocity(iv,2:4),[xm,ym,zm]-position(iv,:),rotationalComponent)
+ radial = [xm,ym,zm]-position(iv,:)
+ uxmid = half*(ux1(i,j,k)+ux1(i+1,j,k)) - linearVelocity(iv,1) + rotationalComponent(1)
+ uymid = half*(uy1(i,j,k)+uy1(i+1,j,k)) - linearVelocity(iv,2) + rotationalComponent(2)
+ uzmid = half*(uz1(i,j,k)+uz1(i+1,j,k)) - linearVelocity(iv,3) + rotationalComponent(3)
fcvx = fcvx +(uymid*radial(3)-uzmid*radial(2))*uymid*dx*dz
fcvy = fcvy +(uzmid*radial(1)-uxmid*radial(3))*uymid*dx*dz
@@ -1506,12 +1506,12 @@ contains
ym=real(jj,mytype)*dz
! write(*,*) 'Calculating force at left x boundary', [xm,ym,zm]
!momentum flux
- radial = [xm,ym,zm]-position
+ radial = [xm,ym,zm]-position(iv,:)
- call crossProduct(angularVelocity(2:4),[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)
+ call crossProduct(angularVelocity(iv,2:4),[xm,ym,zm]-position(iv,:),rotationalComponent)
+ uxmid = half*(ux2(i,j,k)+ux2(i,j+1,k)) - linearVelocity(iv,1) + rotationalComponent(1)
+ uymid = half*(uy2(i,j,k)+uy2(i,j+1,k)) - linearVelocity(iv,2) + rotationalComponent(2)
+ uzmid = half*(uz2(i,j,k)+uz2(i,j+1,k)) - linearVelocity(iv,3) + rotationalComponent(3)
fcvx = fcvx -(uymid*radial(3)-uzmid*radial(2))*uxmid*del_y(j)*dz
fcvy = fcvy -(uzmid*radial(1)-uxmid*radial(3))*uxmid*del_y(j)*dz
@@ -1577,13 +1577,13 @@ contains
ym=real(jj,mytype)*dy
! write(*,*) 'Calculating force at right x boundary', [xm,ym,zm]
- radial = [xm,ym,zm]-position
+ radial = [xm,ym,zm]-position(iv,:)
!momentum flux
- call crossProduct(angularVelocity(2:4),[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)
+ call crossProduct(angularVelocity(iv,2:4),[xm,ym,zm]-position(iv,:),rotationalComponent)
+ uxmid = half*(ux2(i,j,k)+ux2(i,j+1,k)) - linearVelocity(iv,1) + rotationalComponent(1)
+ uymid = half*(uy2(i,j,k)+uy2(i,j+1,k)) - linearVelocity(iv,2) + rotationalComponent(2)
+ uzmid = half*(uz2(i,j,k)+uz2(i,j+1,k)) - linearVelocity(iv,3) + rotationalComponent(3)
fcvx = fcvx +(uymid*radial(3)-uzmid*radial(2))*uxmid*del_y(j)*dz
@@ -1657,12 +1657,12 @@ contains
! write(*,*) 'Calculating force at left z boundary', [xm,ym,zm]
!momentum flux
- radial = [xm,ym,zm]-position
+ radial = [xm,ym,zm]-position(iv,:)
- call crossProduct(angularVelocity(2:4),[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)
+ call crossProduct(angularVelocity(iv,2:4),[xm,ym,zm]-position(iv,:),rotationalComponent)
+ uxmid = half*(ux1(i,j,k)+ux1(i+1,j,k)) - linearVelocity(iv,1) + rotationalComponent(1)
+ uymid = half*(uy1(i,j,k)+uy1(i+1,j,k)) - linearVelocity(iv,2) + rotationalComponent(2)
+ uzmid = half*(uz1(i,j,k)+uz1(i+1,j,k)) - linearVelocity(iv,3) + rotationalComponent(3)
fcvx = fcvx -(uymid*radial(3)-uzmid*radial(2))*uzmid*dx*del_y(j) !!!CHANGE
fcvy = fcvy -(uzmid*radial(1)-uxmid*radial(3))*uzmid*dx*del_y(j)
@@ -1729,14 +1729,14 @@ contains
ii=xstart(1)+i-1
xm=real(ii,mytype)*dx
!momentum flux
- radial = [xm,ym,zm]-position
+ radial = [xm,ym,zm]-position(iv,:)
- call crossProduct(angularVelocity(2:4),[xm,ym,zm]-position,rotationalComponent)
+ call crossProduct(angularVelocity(iv,2:4),[xm,ym,zm]-position(iv,:),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(iv,1) + rotationalComponent(1)
+ uymid = half*(uy1(i,j,k)+uy1(i+1,j,k)) - linearVelocity(iv,2) + rotationalComponent(2)
+ uzmid = half*(uz1(i,j,k)+uz1(i+1,j,k)) - linearVelocity(iv,3) + rotationalComponent(3)
fcvx = fcvx +(uymid*radial(3)-uzmid*radial(2))*uzmid*dx*del_y(j)
fcvy = fcvy +(uzmid*radial(1)-uxmid*radial(3))*uzmid*dx*del_y(j)
@@ -1808,9 +1808,9 @@ contains
mom2 = sum(tunsty(:)) + sum(tconvy(:)) + (-2.0*tconv2_sign+1.0)*sum(tconvy2(:))
mom3 = sum(tunstz(:)) + sum(tconvz(:)) + (-2.0*tconv2_sign+1.0)*sum(tconvz2(:))
- dra1 = -(sum(tdiffx) + sum(tdiffx2) + tp1 - mom1)*(1.0-2.0*torq_flip)
- dra2 = -(sum(tdiffy) + sum(tdiffy2) + tp2 - mom2)*(1.0-2.0*torq_flip)
- dra3 = -(sum(tdiffz) + sum(tdiffz2) + tp3 - mom3)*(1.0-2.0*torq_flip)
+ dra1(iv) = -(sum(tdiffx) + sum(tdiffx2) + tp1 - mom1)*(1.0-2.0*torq_flip)
+ dra2(iv) = -(sum(tdiffy) + sum(tdiffy2) + tp2 - mom2)*(1.0-2.0*torq_flip)
+ dra3(iv) = -(sum(tdiffz) + sum(tdiffz2) + tp3 - mom3)*(1.0-2.0*torq_flip)
! do k=1,zsize(3)
diff --git a/src/genepsi3d.f90 b/src/genepsi3d.f90
index 23475ff5..dcc3e91d 100644
--- a/src/genepsi3d.f90
+++ b/src/genepsi3d.f90
@@ -91,16 +91,66 @@ contains
use ibm_param
use ellipsoid_utils, only: NormalizeQuaternion,ellipInertiaCalculate,ellipMassCalculate
use param
- real(mytype) :: eqr
-
- eqr=(shx*shy*shz)**(1.0/3.0)
- shape=[shx/eqr,shy/eqr,shz/eqr]
-
- orientation=[oriw,orii,orij,orik]
- call NormalizeQuaternion(orientation)
- position=[cex,cey,cez]
- linearVelocity=[lvx,lvy,lvz]
- angularVelocity=[zero,avx,avy,avz]
+ real(mytype) :: eqr, ori_dummy(4)
+ integer :: i,ii,j
+
+ do i =1,nbody
+ ii = (i-1)*3
+ write(*,*) sh(ii+1), sh(ii+2), sh(ii+3)
+ eqr=(sh(ii+1)*sh(ii+2)*sh(ii+3))**(1.0/3.0)
+ if (eqr.lt.0.001) then
+ eqr=1.0
+ endif
+ write(*,*) "Body ", i, ", eqr = ", eqr
+ do j = 1,3
+ shape(i,j) = sh(ii+j)/eqr
+ enddo
+ enddo
+ do i = 1,nbody
+ write(*,*) i, "'s shape = ", shape(i,:)
+ enddo
+
+
+ do i = 1,nbody
+ ii = (i-1)*4
+ do j = 1,4
+ ori_dummy(j) = ori(ii+j)
+ enddo
+
+ ! write(*,*) "Body, ", i, "orientation = ", ori_dummy
+ call NormalizeQuaternion(ori_dummy)
+ write(*,*) "Body, ", i, "orientation = ", ori_dummy
+
+ orientation(i,:) = ori_dummy
+ enddo
+ ! call NormalizeQuaternion(orientation)
+ do i = 1,nbody
+ ii = (i-1)*3
+ do j = 1,3
+ position(i,j) = ce(ii+j)
+ ! write(*,*) ce(i,j), position(j,i)
+ enddo
+ write(*,*) "Nbody", i, "position = ", position(i, :)
+ enddo
+ ! write(*,*) "CE = ", ce
+ ! write(*,*) "Position = ", position
+ ! position=ce
+ do i = 1,nbody
+ ii = (i-1)*3
+ do j = 1,3
+ linearVelocity(i,j) = lv(ii+j)
+ enddo
+ enddo
+
+ do i = 1,nbody
+ ii = (i-1)*3
+ angularVelocity(i,1)=zero
+ do j = 1,3
+ angularVelocity(i,j+1)=av(ii+j)
+ enddo
+ write(*,*) "Nbody", i, "angvel = ", angularVelocity(i, :)
+ enddo
+ write(*,*) "Ra = ", ra
call ellipInertiaCalculate(shape,rho_s,inertia)
call ellipMassCalculate(shape,rho_s,ellip_m)
diff --git a/src/module_param.f90 b/src/module_param.f90
index 90386e2a..49be2679 100644
--- a/src/module_param.f90
+++ b/src/module_param.f90
@@ -617,11 +617,12 @@ 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,rho_s,ellip_m,c_air,cvl_scalar,grav_y,grav_x,grav_z
- real(mytype) :: position(3),orientation(4),linearVelocity(3),angularVelocity(4),linearAcceleration(3),linearForce(3),torque(3),shape(3),inertia(3,3)
- real(mytype) :: position_1(3),linearVelocity_1(3),orientation_1(4),angularVelocity_1(4)
+ real(mytype) :: cex,cey,cez,shx,shy,shz,oriw,orii,orij,orik,lvx,lvy,lvz,avx,avy,avz,ubcx,ubcy,ubcz,rads,ellip_m,c_air,cvl_scalar,grav_y,grav_x,grav_z
+ real(mytype) :: position(10,3),orientation(10,4),linearVelocity(10,3),angularVelocity(10,4),linearAcceleration(3),linearForce(10,3),torque(10,3),shape(10,3),inertia(3,3)
+ real(mytype) :: position_1(10,3),linearVelocity_1(10,3),orientation_1(10,4),angularVelocity_1(10,4),ra(10),rho_s(10)
real(mytype) :: chord,thickness,omega, tconv2_sign, shear_velocity
+ real(mytype) :: ce(30),sh(30),ori(40), lv(30), av(30)
integer :: inana ! Analytical BC as Input
- integer :: imove, nozdrift, force_csv, bodies_fixed, cube_flag, torques_flag,orientations_free, shear_flow_ybc, shear_flow_zbc,torq_debug, torq_flip, ztorq_only
+ integer :: imove, nozdrift, force_csv, bodies_fixed, cube_flag, torques_flag,orientations_free, shear_flow_ybc, shear_flow_zbc,torq_debug, torq_flip, ztorq_only, nbody
end module ibm_param
!############################################################################
diff --git a/src/parameters.f90 b/src/parameters.f90
index bdee3246..6dab4241 100644
--- a/src/parameters.f90
+++ b/src/parameters.f90
@@ -60,10 +60,10 @@ 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, &
+ NAMELIST /ibmstuff/ ce,sh,ori,lv,av,ra, &
nobjmax,nraf,nvol,iforces, cvl_scalar, npif, izap, ianal, imove, thickness, chord, omega , &
ubcx,ubcy,ubcz,rads,rho_s, c_air, grav_x,grav_y,grav_z, nozdrift, force_csv, bodies_fixed, cube_flag, tconv2_sign, &
- torques_flag, orientations_free, shear_flow_ybc, shear_flow_zbc, shear_velocity, torq_debug, torq_flip, ztorq_only
+ torques_flag, orientations_free, shear_flow_ybc, shear_flow_zbc, shear_velocity, torq_debug, torq_flip, ztorq_only, nbody
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, &
@@ -574,7 +574,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,cvl_scalar,grav_y,grav_x,grav_z,nozdrift,force_csv
+ use ibm_param, only : oriw,rho_s,cvl_scalar,grav_y,grav_x,grav_z,nozdrift,force_csv,nbody,ra
implicit none
@@ -622,6 +622,8 @@ subroutine parameter_defaults()
grav_z=zero
nozdrift=0
force_csv=0
+ nbody=1
+ ra(:) = 1.0
!! Gravity field
gravx = zero
diff --git a/src/xcompact3d.f90 b/src/xcompact3d.f90
index ffec99fb..4127f9bd 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,torque_calc
+ use forces, only : force, init_forces, iforces,update_forces, xld,xrd,yld,yud,zld,zrd,torque_calc,nvol
implicit none
- real(mytype) :: dummy,drag,lift,lat,grav_effy,grav_effx,grav_effz,xtorq,ytorq,ztorq,maxrad
+ real(mytype) :: dummy,drag(10),lift(10),lat(10),grav_effy(10),grav_effx(10),grav_effz(10),xtorq(10),ytorq(10),ztorq(10),maxrad
integer :: iounit,ierr,i
real, dimension(100) :: x
@@ -76,20 +76,24 @@ program xcompact3d
if (imove.eq.1) then ! update epsi for moving objects
if ((iibm.eq.2).or.(iibm.eq.3)) then
call genepsi3d(ep1)
- maxrad = max(shape(1),shape(2),shape(3))
- if (iforces.eq.1) then
- xld(1) = position(1) - maxrad * ra * cvl_scalar
- xrd(1) = position(1) + maxrad * ra * cvl_scalar
- yld(1) = position(2) - maxrad * ra * cvl_scalar
- yud(1) = position(2) + maxrad * ra * cvl_scalar
- zld(1) = position(3) - maxrad * ra * cvl_scalar
- zrd(1) = position(3) + maxrad * ra * cvl_scalar
- if (itime.eq.ifirst) then
- call init_forces()
- else
- call update_forces()
+ do i = 1,nobjmax
+ maxrad = max(shape(i,1),shape(i,2),shape(i,3))
+ if (iforces.eq.1) then
+ xld(i) = position(i,1) - maxrad * ra(i) * cvl_scalar
+ xrd(i) = position(i,1) + maxrad * ra(i) * cvl_scalar
+ yld(i) = position(i,2) - maxrad * ra(i) * cvl_scalar
+ yud(i) = position(i,2) + maxrad * ra(i) * cvl_scalar
+ zld(i) = position(i,3) - maxrad * ra(i) * cvl_scalar
+ zrd(i) = position(i,3) + maxrad * ra(i) * cvl_scalar
+ ! write(*,*) "CV bounds = ", xld(i), xrd(i), yld(i), yud(i), zld(i), zrd(i)
+
endif
- endif
+ enddo
+ if (itime.eq.ifirst) then
+ call init_forces()
+ else
+ call update_forces()
+ endif
else if (iibm.eq.1) then
call body(ux1,uy1,uz1,ep1)
endif
@@ -127,19 +131,21 @@ program xcompact3d
grav_effx = grav_x*(rho_s-1.0)
grav_effy = grav_y*(rho_s-1.0)
grav_effz = grav_z*(rho_s-1.0)
- linearForce=[drag-grav_effx,lift-grav_effy,lat-grav_effz]
+ do i = 1,nbody
+ linearForce(i,:) = [drag(i)-grav_effx(i), lift(i)-grav_effy(i), lat(i)-grav_effz(i)]
+ enddo
if (nozdrift==1) then
- linearForce(3)=zero
+ linearForce(:,3)=zero
endif
if (bodies_fixed==1) then
- linearForce(:)=zero
+ linearForce(:,:)=zero
endif
if ((nrank==0).and.(force_csv.eq.1)) then
! open(unit=20, file='force_out.dat', action='write')
- write(20, *) linearForce(1), linearForce(2), linearForce(3)
- write(*,*) 'Writing forces', linearForce(1), linearForce(2), linearForce(3)
+ write(20, *) linearForce(1,1), linearForce(1,2), linearForce(1,3)
+ write(*,*) 'Writing forces', linearForce(1,1), linearForce(1,2), linearForce(1,3)
flush(20)
endif
@@ -148,12 +154,15 @@ program xcompact3d
call torque_calc(ux1,uy1,uz1,ep1,xtorq,ytorq,ztorq,1)
endif
if (orientations_free.eq.1) then
- torque = [xtorq,ytorq,ztorq]
+ do i = 1,nvol
+ torque(i,:) = [xtorq(i), ytorq(i), ztorq(i)]
+ enddo
if (ztorq_only.eq.1) then
- torque = [0.0_mytype, 0.0_mytype, ztorq]
+ torque(:,1) = zero
+ torque(:,2) = zero
endif
else
- torque(:) = zero
+ torque(:,:) = zero
endif
! if (nrank==0) then
@@ -169,10 +178,10 @@ program xcompact3d
- if (nrank==0) then
- write(12 ,*) t, position(1), position(2), position(3), orientation(1), orientation(2), orientation(3), orientation(4), linearVelocity(1), linearVelocity(2), linearVelocity(3), angularVelocity(2), angularVelocity(3), angularVelocity(4), drag, lift, lat, xtorq, ytorq, ztorq
- flush(12)
- endif
+ ! if (nrank==0) then
+ ! write(12 ,*) t, position(1), position(2), position(3), orientation(1), orientation(2), orientation(3), orientation(4), linearVelocity(1), linearVelocity(2), linearVelocity(3), angularVelocity(2), angularVelocity(3), angularVelocity(4), drag, lift, lat, xtorq, ytorq, ztorq
+ ! flush(12)
+ ! endif
! endif
if ((nrank==0).and.(mod(itime,ilist)==0)) then
--
GitLab