added shear flow

src/BC-Ellipsoid.f90

@@ -130,7 +130,7 @@ subroutine inflow (phi)
 
     implicit none
 
-    integer  :: j,k,is
+    integer  :: i,j,k,is
     real(mytype),dimension(xsize(1),xsize(2),xsize(3),numscalar) :: phi
 
     !call random_number(bxo)
@@ -144,6 +144,19 @@ subroutine inflow (phi)
         enddo
     enddo
 
+    if (shear_flow_ybc.eq.1) then 
+        do k=1,xsize(3)
+            do i=1,xsize(1)
+                byxn(i,k)=+shear_velocity
+            enddo
+        enddo
+        do k=1,xsize(3)
+            do i=1,xsize(1)
+                byx1(i,k)=-shear_velocity
+            enddo 
+        enddo 
+    endif   
+
     if (iscalar.eq.1) then
         do is=1, numscalar
         do k=1,xsize(3)
@@ -245,7 +258,7 @@ subroutine init_ellip (ux1,uy1,uz1,phi1)
     real(mytype),dimension(xsize(1),xsize(2),xsize(3)) :: ux1,uy1,uz1
     real(mytype),dimension(xsize(1),xsize(2),xsize(3),numscalar) :: phi1
 
-    real(mytype) :: y,um,eqr
+    real(mytype) :: y,um,eqr,temp_vel
     integer :: k,j,i,ii,is,code
 
     ! write(*,*) 'INSIDE INIT ELLIP'
@@ -279,6 +292,8 @@ subroutine init_ellip (ux1,uy1,uz1,phi1)
 
     ux1=zero; uy1=zero; uz1=zero
 
+    
+
     if (iin.ne.0) then
         call system_clock(count=code)
         if (iin.eq.2) code=0
@@ -293,6 +308,7 @@ subroutine init_ellip (ux1,uy1,uz1,phi1)
         do j=1,xsize(2)
             do i=1,xsize(1)
                 ux1(i,j,k)=init_noise*(ux1(i,j,k)-0.5)
+                
                 uy1(i,j,k)=init_noise*(uy1(i,j,k)-0.5)
                 uz1(i,j,k)=init_noise*(uz1(i,j,k)-0.5)
             enddo
@@ -319,6 +335,11 @@ subroutine init_ellip (ux1,uy1,uz1,phi1)
         do j=1,xsize(2)
         do i=1,xsize(1)
             ux1(i,j,k)=ux1(i,j,k)+u1
+            if (shear_flow_ybc.eq.1) then 
+                temp_vel=(real((j+xstart(2)-1-1),mytype)*dy-yly/2.)/(yly/2.0)*shear_velocity
+                ux1(i,j,k)=ux1(i,j,k)+temp_vel
+                ! write(*,*) "shear velocity at y = ", real((j+xstart(2)-1-1),mytype)*dy, " = ",  temp_vel
+            endif
             uy1(i,j,k)=uy1(i,j,k)
             uz1(i,j,k)=uz1(i,j,k)
         enddo

src/module_param.f90

@@ -620,8 +620,8 @@ module ibm_param
   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) :: chord,thickness,omega, tconv2_sign
+  real(mytype) :: chord,thickness,omega, tconv2_sign, shear_velocity
   integer :: inana ! Analytical BC as Input
-  integer :: imove, nozdrift, force_csv, bodies_fixed, cube_flag
+  integer :: imove, nozdrift, force_csv, bodies_fixed, cube_flag, shear_flow_ybc
 end module ibm_param
 !############################################################################

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, bodies_fixed, cube_flag, tconv2_sign
+      ubcx,ubcy,ubcz,rads,rho_s, c_air, grav_x,grav_y,grav_z, nozdrift, force_csv, bodies_fixed, cube_flag, tconv2_sign, shear_flow_ybc, shear_velocity
   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, &