Fast Marching Method

.gitignore

@@ -18,6 +18,7 @@
 # Mac stuff
 .DS_Store
 *.dSYM/
+.history/
 
 # Temporary files/folders
 _*

examples/amr_ljcf/input

@@ -1,27 +1,34 @@
 # Mesh parameters
 Base nx:   16
 Base ny:   16
-Base nz:   16
-Max level: 6
+Base nz:    8
+Max level:  3
+
+# Domain
+Lx: 16 ! Jet injection direction
+Ly: 16 ! Crossflow
+Lz:  8 ! Spanwise
+Ly offset: 2 ! Distance from edge of domain to jet
 
 # Regridding parameters
 Regrid nsteps:    10
 Tagging Reynolds: 20
 
 # Flow parameters
-Density ratio:   800
-Mom flux ratio:  50
-Weber number:    100
-Reynolds number: 4e4
-Viscosity ratio: 50
-Froude number:   15
+Density ratio:   833.3
+#Mom flux ratio:  50
+Weber number:    108.9
+Reynolds number: 39600 !4e4
+Viscosity ratio: 55.5
+Froude number:   15.12
 
 # Time integration
 Max time: 100
 Max dt:   2.5e-3
-Max CFL:  0.7
+Max CFL:  0.5
 
 # Output
 Output period:     0.5
 Checkpoint period: 5.0
-Restart from:      !restart/jet_3.00024E+01
+CCLabel period:    0.1
+Restart from:      !restart/jet_4.50006E+01

examples/amr_ljcf/src/simulation.f90

@@ -8,6 +8,7 @@ module simulation
    use amrdata_class,     only: amrdata
    use timetracker_class, only: timetracker
    use event_class,       only: event
+   use amrcclabel_class,  only: amrcclabel,stats_type
    use monitor_class,     only: monitor
    use messager,          only: log
    use amrio_class,       only: amrio
@@ -35,6 +36,12 @@ module simulation
    type(event) :: regrid_evt
    real(WP) :: Re_tag=huge(1.0_WP)
 
+   ! CCLabel
+   type(event) :: cclabel_evt
+   type(amrcclabel) :: cclabel
+   type(stats_type), dimension(:), allocatable :: stats 
+   type(monitor) :: cclabel_file 
+
    ! Monitoring
    type(monitor) :: mfile,cflfile,gridfile
 
@@ -50,6 +57,10 @@ module simulation
    real(WP) :: viscL_mol,viscG_mol
    real(WP), dimension(3) :: gravity
 
+   ! Sponge layer parameters for outflow damping
+   real(WP) :: y_spg_start
+   real(WP) :: L_spg
+   real(WP) :: max_cfl_spg
 contains
 
    !> Levelset function for sphere
@@ -78,6 +89,10 @@ subroutine get_viscosity()
       type(amrex_box) :: bx
       real(WP), dimension(:,:,:,:), contiguous, pointer :: pVF,pVisc
       real(WP), parameter :: myeps=1.0e-15_WP
+      real(WP) :: blend,mu_spg,y_loc,nu_spg
+      ! Compute maximum allowable kinematic viscosity in the sponge at finest level
+      ! This ensures viscous terms satisfy CFL: dt <= CFL*dx^2/(4*nu) => nu <= CFL*dx^2/(4*dt)
+      nu_spg=max_cfl_spg*amr%min_meshsize(amr%clvl())**2/(4.0_WP*time%dt)
       ! Loop over levels
       do lvl=0,amr%clvl()
          ! Loop over domain
@@ -91,12 +106,52 @@ subroutine get_viscosity()
             do k=bx%lo(3),bx%hi(3); do j=bx%lo(2),bx%hi(2); do i=bx%lo(1),bx%hi(1)
                ! Use harmonic averaging
                pVisc(i,j,k,1)=1.0_WP/(pVF(i,j,k,1)/max(viscL_mol,myeps)+(1.0_WP-pVF(i,j,k,1))/max(viscG_mol,myeps))
+               ! Apply sponge layer viscosity damping in outflow region (y+ boundary)
+               y_loc=amr%ylo+(real(j,WP)+0.5_WP)*amr%dy(lvl)
+               if (y_loc.gt.y_spg_start) then
+                  ! Smooth quadratic blend: 0 at y_spg_start, 1 at y_spg_start+L_spg
+                  blend=min((y_loc-y_spg_start)/L_spg,1.0_WP)**2
+                  ! Compute sponge viscosity (convert kinematic to dynamic)
+                  mu_spg=nu_spg/(pVF(i,j,k,1)/max(1.0_WP,myeps)+(1.0_WP-pVF(i,j,k,1))/max(1.0_WP,myeps))
+                  ! Only increase viscosity if beneficial
+                  pVisc(i,j,k,1)=max(pVisc(i,j,k,1),blend*mu_spg)
+               end if
             end do; end do; end do
          end do
          call amr%mfiter_destroy(mfi)
       end do
    end subroutine get_viscosity
 
+   !> Clip VOF near outflow to prevent unphysical values and instabilities
+   subroutine clip_vof_outflow()
+      use amrex_amr_module, only: amrex_mfiter,amrex_box
+      integer :: lvl,i,j,k
+      type(amrex_mfiter) :: mfi
+      type(amrex_box) :: bx
+      real(WP), dimension(:,:,:,:), contiguous, pointer :: pVF
+      real(WP) :: dy_loc
+      ! Loop over levels
+      do lvl=0,amr%clvl()
+         ! Loop over domain
+         call amr%mfiter_build(lvl,mfi)
+         do while (mfi%next())
+            ! Get pointers to data
+            pVF=>fs%VF%mf(lvl)%dataptr(mfi)
+            ! Get tilebox (grow for safety, clip only at boundaries)
+            bx=mfi%tilebox()
+            do k=bx%lo(3),bx%hi(3); do j=bx%lo(2),bx%hi(2); do i=bx%lo(1),bx%hi(1)
+               ! Check if in sponge layer outflow region (y+ boundary)
+               dy_loc=amr%ylo+(real(j,WP)+0.5_WP)*amr%dy(lvl)
+               ! Remove liquid in sponge region to prevent outflow instabilities
+               if (dy_loc.gt.y_spg_start) then
+                  pVF(i,j,k,1)=0.0_WP
+               end if
+            end do; end do; end do
+         end do
+         call amr%mfiter_destroy(mfi)
+      end do
+   end subroutine clip_vof_outflow
+
    !> Tagger for this case based on velocity gradient magnitude
    subroutine my_tagger(solver,lvl,time,tags_ptr)
       use iso_c_binding,    only: c_ptr,c_char
@@ -130,6 +185,8 @@ subroutine my_tagger(solver,lvl,time,tags_ptr)
             ! Prevent maximum Re-driven refinement near wall
             near_wall=(solver%amr%xlo+(real(i,WP)+0.5_WP)*dx.lt.solver%amr%xlo+2.0_WP*dx)
             if (near_wall.and.lvl.ge.solver%amr%maxlvl-1) cycle
+            ! Prevent refinement in outflow sponge region
+            if (solver%amr%ylo+(real(j,WP)+0.5_WP)*dy.gt.y_spg_start) cycle
             ! Laplacian of velocity Q=UVW
             lapU=(pQ(i+1,j,k,1)-2.0_WP*pQ(i,j,k,1)+pQ(i-1,j,k,1))*dxi2+(pQ(i,j+1,k,1)-2.0_WP*pQ(i,j,k,1)+pQ(i,j-1,k,1))*dyi2+(pQ(i,j,k+1,1)-2.0_WP*pQ(i,j,k,1)+pQ(i,j,k-1,1))*dzi2
             lapV=(pQ(i+1,j,k,2)-2.0_WP*pQ(i,j,k,2)+pQ(i-1,j,k,2))*dxi2+(pQ(i,j+1,k,2)-2.0_WP*pQ(i,j,k,2)+pQ(i,j-1,k,2))*dyi2+(pQ(i,j,k+1,2)-2.0_WP*pQ(i,j,k,2)+pQ(i,j,k-1,2))*dzi2
@@ -160,6 +217,7 @@ subroutine dirichlet_velocity(solver,lvl,time,face,bx,comp,p)
        case (1)  ! Inflow in X-
          select case (comp)
           case ('U')  ! Staggered U=Ujet
+            Ujet = gravity(1)*time
             do k=bx%lo(3),bx%hi(3); do j=bx%lo(2),bx%hi(2); do i=bx%lo(1),bx%hi(1)
                rad=sqrt((amr%ylo+(real(j,WP)+0.5_WP)*amr%dy(lvl))**2+(amr%zlo+(real(k,WP)+0.5_WP)*amr%dz(lvl))**2)
                if (amr%nz.eq.1) rad=sqrt((amr%ylo+(real(j,WP)+0.5_WP)*amr%dy(lvl))**2)
@@ -318,20 +376,143 @@ subroutine jet_init(solver,lvl,time,ba,dm)
       call amrex_mfiter_destroy(mfi)
    end subroutine jet_init
 
+   !> Function that identifies cells within a structure
+   logical function make_label(pVF,lo,i,j,k)
+      implicit none
+      real(WP), dimension(:,:,:,:), intent(in) :: pVF
+      integer, dimension(3), intent(in) :: lo
+      integer, intent(in) :: i,j,k
+      integer :: il,jl,kl
+      il = i - lo(1) + 1
+      jl = j - lo(2) + 1
+      kl = k - lo(3) + 1
+      if (pVF(il,jl,kl,1).gt.0.0_WP) then
+         make_label=.true.
+      else
+         make_label=.false.
+      end if
+   end function make_label
+
+   !> Function that identifies if neighbors are within the same structure
+   logical function same_label(pVF,lo,i,j,k,ii,jj,kk)
+      implicit none
+      real(WP), dimension(:,:,:,:), intent(in) :: pVF
+      integer, dimension(3), intent(in) :: lo
+      integer, intent(in) :: i,j,k,ii,jj,kk
+      integer :: il,jl,kl,iil,jjl,kkl
+      il  = i  - lo(1) + 1
+      jl  = j  - lo(2) + 1
+      kl  = k  - lo(3) + 1
+      iil = ii - lo(1) + 1
+      jjl = jj - lo(2) + 1
+      kkl = kk - lo(3) + 1
+      if (pVF(il,jl,kl,1).gt.0.0_WP .and. pVF(iil,jjl,kkl,1).gt.0.0_WP) then
+         same_label=.true.
+      else
+         same_label=.false.
+      end if
+   end function same_label
+
+   !> Function that identifies cells within a structure on coarse level
+   logical function coarse_make_label(pVF,lo,i,j,k)
+      use amrmpinc_class,   only: VFhi
+      implicit none
+      real(WP), dimension(:,:,:,:), intent(in) :: pVF
+      integer, dimension(3), intent(in) :: lo
+      integer, intent(in) :: i,j,k
+      integer :: il,jl,kl
+      il = i - lo(1) + 1
+      jl = j - lo(2) + 1
+      kl = k - lo(3) + 1
+      if (pVF(il,jl,kl,1).gt.VFhi) then
+         coarse_make_label=.true.
+      else
+         coarse_make_label=.false.
+      end if
+   end function coarse_make_label
+
+   !> Function that identifies if neighbors are within the same structure on coarse level
+   logical function coarse_same_label(pVF,lo,i,j,k,ii,jj,kk)
+      use amrmpinc_class,   only: VFhi
+      implicit none
+      real(WP), dimension(:,:,:,:), intent(in) :: pVF
+      integer, dimension(3), intent(in) :: lo
+      integer, intent(in) :: i,j,k,ii,jj,kk
+      integer :: il,jl,kl,iil,jjl,kkl
+      il  = i  - lo(1) + 1
+      jl  = j  - lo(2) + 1
+      kl  = k  - lo(3) + 1
+      iil = ii - lo(1) + 1
+      jjl = jj - lo(2) + 1
+      kkl = kk - lo(3) + 1
+      if (pVF(il,jl,kl,1).gt.VFhi .and. pVF(iil,jjl,kkl,1).gt.VFhi) then
+         coarse_same_label=.true.
+      else
+         coarse_same_label=.false.
+      end if
+   end function coarse_same_label
+
+   !> Write droplet statistics to monitor files
+   subroutine write_stats()
+      use monitor_class, only: iformat,rformat
+      use string,    only: str_medium
+      implicit none
+      character(len=str_medium) :: filename,struct_name
+      type(stats_type) :: buf   ! single-structure buffer monitor points into
+      integer :: n
+      ! Create a file to write Weber numbers
+      write(filename, rformat) time%t
+      filename = 'structStats_'//trim(adjustl(filename))
+      cclabel_file=monitor(fs%amr%amRoot,filename)
+
+      ! Register columns with buffer
+      call cclabel_file%add_column(buf%id,      'Structure ID')
+      call cclabel_file%add_column(buf%vol,     'Drop Volume')
+      call cclabel_file%add_column(buf%Deq,     'Equiv Diameter')
+      call cclabel_file%add_column(buf%com(1),  'X Drop Pos')
+      call cclabel_file%add_column(buf%com(2),  'Y Drop Pos')
+      call cclabel_file%add_column(buf%com(3),  'Z Drop Pos')
+      call cclabel_file%add_column(buf%vel(1),  'X Drop Vel')
+      call cclabel_file%add_column(buf%vel(2),  'Y Drop Vel')
+      call cclabel_file%add_column(buf%vel(3),  'Z Drop Vel')
+      call cclabel_file%add_column(buf%gvel(1), 'X Gas Vel')
+      call cclabel_file%add_column(buf%gvel(2), 'Y Gas Vel')
+      call cclabel_file%add_column(buf%gvel(3), 'Z Gas Vel')
+      call cclabel_file%add_column(buf%moi(1,1),'Ixx')
+      call cclabel_file%add_column(buf%moi(2,2),'Iyy')
+      call cclabel_file%add_column(buf%moi(3,3),'Izz')
+      call cclabel_file%add_column(buf%moi(1,2),'Ixy')
+      call cclabel_file%add_column(buf%moi(1,3),'Ixz')
+      call cclabel_file%add_column(buf%moi(2,3),'Iyz')
+      call cclabel_file%add_column(buf%weber,   'Weber')
+      do n=1,cclabel%nstruct
+         ! Set buffer and write the data for this structure
+         buf = stats(n)
+         call cclabel_file%write()
+      end do
+      ! Close file
+      call cclabel_file%close()
+   end subroutine write_stats
+
    !> Initialization hook
    subroutine simulation_init()
       use param, only: param_read
       implicit none
 
       ! Create amrgrid
       create_amrgrid: block
+         real(WP) :: Lx,Ly,Lz,Ly_offset
          amr%name='LJCF'
          call param_read('Base nx',amr%nx)
          call param_read('Base ny',amr%ny)
          call param_read('Base nz',amr%nz)
-         amr%xlo= 00.0_WP; amr%xhi=+20.0_WP
-         amr%ylo=-05.0_WP; amr%yhi=+15.0_WP
-         amr%zlo=-10.0_WP; amr%zhi=+10.0_WP
+         call param_read("Lx",Lx)
+         call param_read("Ly",Ly)
+         call param_read("Lz",Lz)
+         call param_read("Ly offset",Ly_offset)
+         amr%xlo= 00.0_WP; amr%xhi=+Lx
+         amr%ylo=-Ly_offset; amr%yhi=Ly-Ly_offset
+         amr%zlo=-Lz/2.0_WP; amr%zhi=+Lz/2.0_WP
          amr%xper=.false.; amr%yper=.false.; amr%zper=.true.
          call param_read('Max level',amr%maxlvl)
          ! Handle 2D case
@@ -384,14 +565,18 @@ subroutine simulation_init()
          ! Set densities
          fs%rhoG=1.0_WP; call param_read('Density ratio',fs%rhoL)
          ! Read in momentum flux ratio and set liquid velocity
-         call param_read('Mom flux ratio',Ujet); Ujet=sqrt(Ujet/fs%rhoL)
+         ! call param_read('Mom flux ratio',Ujet); Ujet=sqrt(Ujet/fs%rhoL)
          ! Set surface tension coefficient
          call param_read('Weber number',fs%sigma); fs%sigma=1.0_WP/fs%sigma
          ! Set molecular viscosities
          call param_read('Reynolds number',viscG_mol); viscG_mol=1.0_WP/viscG_mol
          call param_read('Viscosity ratio',viscL_mol); viscL_mol=viscG_mol*viscL_mol
          ! Set gravity
          gravity=0.0_WP; call param_read('Froude number',gravity(1),default=1.0e30_WP); gravity(1)=1.0_WP/gravity(1)**2
+         ! Set sponge layer parameters (optional for outflow damping at y+ boundary)
+         call param_read('Sponge y-start',y_spg_start,default=12.0_WP)
+         call param_read('Sponge thickness',L_spg,default=4.0_WP)
+         call param_read('Sponge max CFL',max_cfl_spg,default=0.5_WP)
          ! Set pressure convergence
          fs%psolver%outer_solver=amrmg_outer_pcg_mlmg
          fs%psolver%tol_rel=1.0e-5_WP
@@ -460,6 +645,14 @@ subroutine simulation_init()
          call io%add_scalar(name='dt',value=time%dt)
       end block init_checkpoint
 
+      ! Initialize CClabel
+      cclabel_evt=event(time=time,name="CCLabel output")
+      call param_read('CCLabel period',cclabel_evt%tper)
+      call cclabel%initialize(amr,name='amr_ljcf')
+      call cclabel%build(make_label,same_label,coarse_make_label,coarse_same_label,fs%VF) 
+      call cclabel%compute_stats(fs%VF, fs%Q, fs%rhoG, fs%sigma, stats)
+      call write_stats()
+
       ! Initialize visualization
       create_visualization: block
          ! Create visualization object
@@ -471,6 +664,7 @@ subroutine simulation_init()
          call viz%add_scalar(fs%visc,1,'visc')
          call viz%add_scalar(fs%P,1,'pressure')
          call viz%add_scalar(fs%VF,1,'VF')
+         call viz%add_scalar(cclabel%id,1,'ID')
          call viz%add_surfmesh(fs%smesh,'plic')
          ! Create visualization output event
          viz_evt=event(time=time,name='Visualization output')
@@ -563,6 +757,9 @@ subroutine simulation_run()
             call fs%build_plic(time%t)
             call fs%build_subVF()
 
+            ! Clip VOF in outflow region to prevent instabilities
+            ! call clip_vof_outflow()
+
             ! Interpolate velocity to the faces
             call fs%get_face_velocity()
 
@@ -612,6 +809,13 @@ subroutine simulation_run()
          ! Compute Umag
          call Umag%get_magnitude(srcX=fs%Q,srcY=fs%Q,srcZ=fs%Q,compX=1,compY=2,compZ=3)
 
+         ! Construct CCLabel then compute & write stats
+         if (cclabel_evt%occurs()) then 
+            call cclabel%build(make_label,same_label,coarse_make_label,coarse_same_label,fs%VF)
+            call cclabel%compute_stats(fs%VF, fs%Q, fs%rhoG, fs%sigma, stats)
+            call write_stats()
+         end if
+
          ! Monitor output
          call fs%get_info()
          call mfile%write()