From fd1d9699433e7bd18924dc6960da6d4a8f1f9fd7 Mon Sep 17 00:00:00 2001
From: adfboyd <adfboyd@UbuntuVM.myguest.virtualbox.org>
Date: Thu, 21 Mar 2024 17:41:28 +0000
Subject: [PATCH] added some extra debugging parameters
---
src/BC-Ellipsoid.f90 | 8 ++---
src/module_param.f90 | 2 +-
src/parameters.f90 | 2 +-
src/xcompact3d.f90 | 85 +++++++++++++++++++++++++++++++++++++-------
4 files changed, 78 insertions(+), 19 deletions(-)
diff --git a/src/BC-Ellipsoid.f90 b/src/BC-Ellipsoid.f90
index ba4a83f0..ec64576c 100644
--- a/src/BC-Ellipsoid.f90
+++ b/src/BC-Ellipsoid.f90
@@ -313,10 +313,10 @@ subroutine init_ellip (ux1,uy1,uz1,phi1)
if (sine_z.eq.1) then
uz1(i,j,k)=-sin_prec(x*twopi)*-cos_prec(y*twopi)*cos_prec(z*twopi)
endif
- if (iscalar == 1) then
- phi1(i,j,k,1:numscalar)=sin_prec(x)*sin_prec(y)*cos_prec(z)
- endif
- uz1(i,j,k)=zero
+ ! if (iscalar == 1) then
+ ! phi1(i,j,k,1:numscalar)=sin_prec(x)*sin_prec(y)*cos_prec(z)
+ ! endif
+ ! uz1(i,j,k)=zero
enddo
enddo
enddo
diff --git a/src/module_param.f90 b/src/module_param.f90
index 284c8901..788921ef 100644
--- a/src/module_param.f90
+++ b/src/module_param.f90
@@ -623,6 +623,6 @@ module ibm_param
real(mytype) :: chord,thickness,omega
integer :: sine_init, sine_x, sine_y, sine_z
integer :: inana ! Analytical BC as Input
- integer :: imove, nozdrift, force_csv
+ integer :: imove, nozdrift, force_csv, bodies_fixed
end module ibm_param
!############################################################################
diff --git a/src/parameters.f90 b/src/parameters.f90
index c163992f..3d4da04c 100644
--- a/src/parameters.f90
+++ b/src/parameters.f90
@@ -62,7 +62,7 @@ subroutine parameter(input_i3d)
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, 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
+ ubcx,ubcy,ubcz,rads,rho_s, c_air, grav_x,grav_y,grav_z, nozdrift, force_csv, bodies_fixed
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, &
diff --git a/src/xcompact3d.f90 b/src/xcompact3d.f90
index c47b94b6..f59b6bcb 100644
--- a/src/xcompact3d.f90
+++ b/src/xcompact3d.f90
@@ -18,10 +18,13 @@ program xcompact3d
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 dbg_schemes, only : sin_prec, cos_prec
implicit none
+! real(mytype),dimension(xsize(1),xsize(2),xsize(3)):: ux1_exact,uy1_exact,uz1_exact
+
real(mytype) :: dummy,drag,lift,lat,grav_effy,grav_effx,grav_effz
- integer :: iounit,ierr,i
- real, dimension(100) :: x
+ integer :: iounit,ierr,i,j,k
+ real(mytype) :: x,y,z
call init_xcompact3d()
@@ -37,12 +40,37 @@ program xcompact3d
write(*,*) 'Outputting forces'
end if
- if (sine_init.eq.1) then
- write(*,*) 'Generated Sinusoidal flow field.'
- endif
+! if (sine_init.eq.1) then
+! write(*,*) 'Generated Sinusoidal flow field.'
+
+! do k=1,xsize(3)
+! z=real((k+xstart(3)-1-1),mytype)*dz
+! do j=1,xsize(2)
+! y=real((j+xstart(2)-1-1),mytype)*dy
+! do i=1,xsize(1)
+! x=real(i-1,mytype)*dx
+
+! if (sine_x.eq.1) then
+! ux1_exact(i,j,k)=-twopi * cos_prec(x*twopi)*cos_prec(y*twopi)*cos_prec(z*twopi)
+! endif
+! if (sine_y.eq.1) then
+! uy1_exact(i,j,k)=-twopi * cos_prec(x*twopi)*cos_prec(y*twopi)*cos_prec(z*twopi)
+! endif
+! if (sine_z.eq.1) then
+! uz1_exact(i,j,k)= twopi * sin_prec(x*twopi)*-cos_prec(y*twopi)*sin_prec(z*twopi)
+! endif
+! ! if (iscalar == 1) then
+! ! phi1(i,j,k,1:numscalar)=sin_prec(x)*sin_prec(y)*cos_prec(z)
+! ! endif
+! ! uz1(i,j,k)=zero
+! enddo
+! enddo
+! enddo
+! endif
+
! do i = 1,100
! x(i) = i
! enddo
@@ -131,6 +159,10 @@ program xcompact3d
linearForce(3)=zero
endif
+ if (bodies_fixed==1) then
+ 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), '\n'
@@ -169,14 +201,41 @@ program xcompact3d
! 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)
-
- call simu_stats(3)
-
- call postprocessing(rho1,ux1,uy1,uz1,pp3,phi1,ep1)
+ if (sine_init.eq.1) then
+ !x-derivatives
+ call derx (ta1,ux1,di1,sx,ffx,fsx,fwx,xsize(1),xsize(2),xsize(3),0,1) !ubcx is 1. etc
+ call derx (tb1,uy1,di1,sx,ffxp,fsxp,fwxp,xsize(1),xsize(2),xsize(3),1,2)
+ call derx (tc1,uz1,di1,sx,ffxp,fsxp,fwxp,xsize(1),xsize(2),xsize(3),1,3)
+ !y-derivatives
+ call dery (ta2,ux2,di2,sy,ffyp,fsyp,fwyp,ppy,ysize(1),ysize(2),ysize(3),1,1)
+ call dery (tb2,uy2,di2,sy,ffy,fsy,fwy,ppy,ysize(1),ysize(2),ysize(3),0,2)
+ call dery (tc2,uz2,di2,sy,ffyp,fsyp,fwyp,ppy,ysize(1),ysize(2),ysize(3),1,3)
+ !!z-derivatives
+ call derz (ta3,ux3,di3,sz,ffzp,fszp,fwzp,zsize(1),zsize(2),zsize(3),1,1)
+ call derz (tb3,uy3,di3,sz,ffzp,fszp,fwzp,zsize(1),zsize(2),zsize(3),1,2)
+ call derz (tc3,uz3,di3,sz,ffz,fsz,fwz,zsize(1),zsize(2),zsize(3),0,3)
+ !!all back to x-pencils
+ call transpose_z_to_y(ta3,td2)
+ call transpose_z_to_y(tb3,te2)
+ call transpose_z_to_y(tc3,tf2)
+ call transpose_y_to_x(td2,tg1)
+ call transpose_y_to_x(te2,th1)
+ call transpose_y_to_x(tf2,ti1)
+ call transpose_y_to_x(ta2,td1)
+ call transpose_y_to_x(tb2,te1)
+ call transpose_y_to_x(tc2,tf1)
+ !du/dx=ta1 du/dy=td1 and du/dz=tg1
+ !dv/dx=tb1 dv/dy=te1 and dv/dz=th1
+ !dw/dx=tc1 dw/dy=tf1 and dw/dz=ti1
+ endif
+
+ enddo !! End sub timesteps
+
+ call restart(ux1,uy1,uz1,dux1,duy1,duz1,ep1,pp3(:,:,:,1),phi1,dphi1,px1,py1,pz1,rho1,drho1,mu1,1)
+
+ call simu_stats(3)
+
+ call postprocessing(rho1,ux1,uy1,uz1,pp3,phi1,ep1)
enddo !! End time loop
--
GitLab