Colliding Flows Module

Input Parameters

cloud.data

nCloud = ! Initial density of grid in c.u. (right now this is the entire grid, but could just be for incoming flows if you wanted to set the ambient density to something else, would need an extra parameter or object)
TCloud = ! Initial temperature of grid in real units

! Velocity - in computational units
vx0 = ! Median velocity in the x-direction for the flow coming from the left.  The flow coming from the right has the exact same speed in the opposite direction.
vy0 = ! Same as vx0, but for y-direction
vz0 = ! Same as vx0, but for z-direction

! k values - for geometric (sinusoidal) interface only.  If iseed is not 0, these will not be used.
kx0 = ! x-direction wavenumber 
ky0 = ! y-direction wavenumber 
kz0 = ! z-direction wavenumber 

! For collision area, in c.u.
bmin    =  ! For 3D ONLY - minor axis of elliptical region.  Set this to 0 if you want flow to cover entire domain.  For a circular profile, set bmin = amaj.
amaj    =  ! For 3D ONLY - major axis of elliptical region.  Set this to 0 if you want flow to cover entire domain.
sig0    =  ! For a velocity profile that merges more smoothly with surrounding region.  Set to be positive for a Gaussian profile, and negative for a double tanh as a taper for a viscous profile. Set to 0 for a box profile.
pert    =  ! Overall amplitude of perturbation - currently measured as a fraction of the x-length of the domain.
radius  =  ! For 2D ONLY- radius of 'circular' region.  Set this to 0 if you want flow to cover entire domain.
slope   =  ! Slope of interface - for introducing shear.

! For random perturbative interface
iseed   =  ! Seed for random array generation.  Set this to 0 if you want a geometric interface instead of a random interface
power   =  ! Power of random amplitude generation
kmin    =  ! Minimum wavenumber
kmax    =  ! Maximum wavenumber

Sample cloud.data file

&CloudData     ! NameList declaration
!
!
nCloud  =  1.00004831529d0,        ! Initial Cloud Density - (n,T) at equilibrium point for IICooling curve.
TCloud  =  5006.0d0,               ! Initial Cloud Temperature in (K)
vx0     =  27.4025546254d0,        ! Median Velocity x-direction
vy0     =  0.d0,
vz0     =  0.d0,
kx0     =  0d0,
ky0     =  3d0,
kz0     =  0d0,
!
!
!
iseed   =  -151,
power   =  -1.d-2,
kmin    =  1,
kmax    =  4, 
bmin    =  1.2d0,                ! elliptical input profile
amaj    =  1.6d0,
sig0    =  2.d-1,                ! double tanh profile
pert    =  0.05d0,               ! 0.05 * xlen
radius  =  0.d0,
slope   =  0.176326980708d0      ! Slope of 10 degrees (tan(10))
!
!
/      ! This line MUST be present!

Slope

Having a sloped interface allows us to introduce shear to the flows. In this case, I simply add slope*y to arlo, and subtract slope*y from arhi, where slope = tan(theta), and theta is the angle away from the vertical. Because of my choice of theta as measured from the vertical axis instead of horizontal, and because arlo and arhi are defined w.r.t. x, the 'line' is x = slope * y.

Current Routines

scaleClouds

Routine to read in cloud.data

InitializeClouds

Routine that creates random amplitude perturbations, in the array amplitrand(1:mytot,1:mztot,1). For an interface with velocity in the x-direction, the amplitude variations are dependent on y and z. Calls fldgen routine included in fldgen.f90, which calls genak routine included in genak.f90 and nr_fourn routine included in nr_fourn.f90. genak also calls function nr_randnum, included in nr_randnum.f90. (see bottom of clouds.f90)

xlen (domain%XUpper(1) - domain%XLower(1)) has been important in this routine because the amplitude is scaled by a fraction of this length. However, it might be a good idea to change this to a real number that you input as pert in computational units.

qinitClouds

2 cases: 2D and 3D

  2D:
     IF (iseed.eq.0) THEN
          (Determine "lower" and "upper" area of interface for this cell using geometrical sine formula)
         IF (kz.gt.0) (overlap perturbation at different scale)
     ELSE
          (use amplitrand array to determine arlo and arhi)
     ENDIF

     (Set density, and determine whether px and py is negative or positive by using arlo and arhi to determine whether this cell is to the left or right of boundary)

     IF (radius.ne.0) THEN  ! Limit inflow region to part of a side.
          IF (dabs(y).LE.radius) THEN       !!!! This could be better coded.  This algorithm inherently assumes that y=0 is the center of your "radius" in 2D.  Perhaps a better way would be to introduce a parameter that sets the center.
             wght=1.d0
          ELSE
             wght=0.d0
          END IF
     ELSE 
          wght = 1.d0
     END IF
 
     q(i,j,1,1,2) = q(i,j,1,1,2) * wght

     (Next set tracers.  LO tracer is the tracer to the left (anywhere px > 0).  LO+1 is the tracer to the right (px < 0).  Everywhere else these are set to a Really Low Number(TM), or in this case, 1e-10.)
     (In afterstep routine, this will be the tracer floor everywhere. )

     (Then call CalcComp to determine energy from TCloud) 


  3D:
     (Very similar to 2D setup, with geometrical collision formula in 3D.  Differences are in limited inflow setup:)

     IF ((bmin .NE. 0d0) .AND. (amaj .NE. 0d0)) THEN  ! Limit inflow region to part of a side
        r2 = ((y/amaj)**2+(z/bmin)**2) ! normalized elliptical "radius"

        IF (sig0 .EQ. 0d0) THEN ! box profile
            (straight ellipse, no funny stuff)

        ELSE IF (sig0 .LT. 0d0) THEN ! Gaussian profile
             wght = dexp(-r2/(sig0/xlen)**2)
        ELSE IF (sig0 .GT. 0d0) THEN ! double tanh profile as a taper for viscous profile (BGK!) 
               r2 = dsqrt(r2)
               wght = 0.5d0*(1d0-(dexp((r2-1d0)*xlen/sig0)-dexp(-(r2-1d0)*xlen/sig0))/(dexp((r2-1d0)*xlen/sig0)+dexp(-(r2-1d0)*xlen/sig0))) 
        END IF
         q(i,j,k,1,2) = q(i,j,k,1,2) * wght
     END IF

    (Tracers)

    (CalcComp)

''Note: There is a commented section that provides yet another option to change the interface, however it uses the parameter 'radius', which has been previously used in 2D to mean something totally different.  One of these should probably be changed so this can be an option in the future''
              !IF (radius .GT. 0d0) THEN ! mimicking taper due to expanding blast wave
              !   cosa = dcos(datan(dsqrt(y**2 + z**2)/radius))
              !END IF

b4stepClouds

Nothing important in b4step at this point.

afterstepClouds

For both 2D and 3D, this routine sets a floor on both metallicity tracers, if they exist.

CloudsBC

This is where the boundary conditions are user-defined.

On the left and right, where the boundary is within the original 2D radius or 3D elliptical radius, the original conditions are reset. Outside of this (if it exists), the B.C. are set to be outgoing-only. All incoming velocity is set to 0, with a section to ensure that the energy is appropriately modified (temperature stays the same).

On the top and bottom, and front and back in the case of 3D, the B.C. are again outgoing only.

Note: This setup works with extrapolating boundary conditions too, but Funny Results ensue