Turbulent Round Jet

[Liuhaosheng1]

Hello author,
I couldn't find the runtime code for Turbulent Round Jet in the FLOWVPM case. It only contains its functions and simulation. Could you please provide me with its runtime code.

[AndreaRapi]

Hi @Liuhaosheng1

Out of curiosity, which case are you talking about in particular?

[EdoAlvarezR]

@AndreaRapi, he is probably referring to the round jet validation case from my dissertation:

I did some digging through my dissertation files, and it seems like this is the driver file that I used to run the simulation (the first line only works if you copy and paste this code into a file). It should still run if you roll back FLOWVPM and possibly FLOWExaFMM to the version that I released when I defended my dissertation. Have fun!

this_file_name = splitext(splitdir(@__FILE__)[2])[1]

#=

    CHANGE LOG
    * max_zd=0.85
    * Reformulated VPM
    * Two-Level SFS
    * Classic VPM, No SFS
    * Remove weak particles
    * SFS two-level
    * max_zd=0.925, U1=40.0
    * max_zd=1.0, U1=42.5
    * max_zd=0.5
    * Quinn's experimental conditions
    * Full cylinder RBF, max_zd=1.0
    * Pedrizzetti's relaxation

    TODO
    * []
=#

# ------------ IMPORT PACKAGES -------------------------------------------------
println("Importing packages...")
import Pkg

using Revise
import FLOWVPM
vpm = FLOWVPM

vpm_path = joinpath(dirname(pathof(FLOWVPM)), "..")*"/"
include(vpm_path*"examples/roundjet/roundjet.jl")

# extdrive_path = "/media/flowlab/Storage/ealvarez/simulations/";
# extdrive_path = "/media/edoalvar/Samsung_T51/simulationdata202108/"
extdrive_path = "/fslhome/edoalvar/simulationdata/"

restart_file = nothing
# restart_file = extdrive_path*"roundjet_val030/"
restart_sigma = nothing

if restart_file!=nothing
    maxnum = maximum([parse(Int, f[findfirst('.', f)+1:findlast('.', f)-1])
                                for f in readdir(restart_file) if contains(f, ".h5")])
    restart_file *= "roundjet.$(maxnum)"
end

# ------------ INPUTS AND OUTPUTS ----------------------------------------------
save_path = extdrive_path*this_file_name*"/"  # Where to save the simulation

run_name  = "roundjet"
prompt    = true
debug     = true
verbose   = true
v_lvl     = 0


# -------------- SIMULATION PARAMETERS -----------------------------------------
# Jet geometry
d         = 45.4e-3                     # (m) jet diameter
O         = zeros(3)                    # Origin of jet
Oaxis     = Float64[i==j for i in 1:3, j in 1:3] # Orientation of jet

# Jet characteristics
# U1        = 60.0                        # (m/s) Inflow centerline velocity
U1        = 40.0
U2        = 0.0                         # (m/s) Coflow velocity
U2angle   = [0, 0, 0]                   # (deg) Coflow angle from centerline
thetaod   = 0.025                       # Ratio of inflow momentum thickness of shear layer to diameter, θ/d
Re        = 2.08e5                      # Reynolds number Re = (U1-U2z)*d / ν

# Temporal discretization
steps_per_d  = 50                       # Number of time steps for the centerline at U1 to travel one diameter
# steps_per_d  = 100
d_travel_tot = 60                       # Run simulation for an equivalent of this many diameters

# Spatial discretization
dxotheta   = 1/4                        # Distance Δx between particles over momentum thickness θ
overlap    = 2.4                        # Overlap between particles
minWfraction = 1e-2                     # Threshold at which not to add particles


# Initialization parameters
# max_zd     = 1.0                        # Maximum diameters in z-direction to create annulis for defining BC
max_zd     = 0.85
max_zsigma = max_zd/(overlap*dxotheta*thetaod)  # Maximum sigmas in z-direction to create annulis for defining BC
rbf        = true                       # RBF on the boundary condition for increased accuracy
rbf_optargs= [(:itmax,200), (:tol,2.5e-2), (:iterror,true), (:verbose,true), (:debug,true)]
rbf_Wf     = 1.0                        # Scales target RBF vorticity by this factor
# rbf_Wf     = 0.5
# rbf_fullcylinder = false                # Use a cylinder of particles for the RBF with half of it in inside the jet
rbf_fullcylinder = true
keep_rbfparticles = true              # Keep particles used for RBF as static particles
# keep_rbfparticles = false
keep_rbffullcylinder = false            # Keep the half of the cylinder that is inside the jet as free particles
# keep_rbffullcylinder = true

# remove_particles_weak = false           # Whether to remove particles with a strength lower than certain threshold
remove_particles_weak = true

# -------------- SFS SCHEME ------------------------------------------------

# SFS = vpm.SFS_none
# SFS = vpm.SFS_Cs_nobackscatter
SFS = vpm.SFS_Cd_twolevel_nobackscatter
# SFS = vpm.SFS_Cd_threelevel_nobackscatter

# SFS = vpm.DynamicSFS(vpm.Estr_fmm, vpm.pseudo3level_positive; clippings=[],
# SFS = vpm.DynamicSFS(vpm.Estr_fmm, vpm.pseudo3level_positive; clippings=[vpm.clipping_backscatter],
#                                                         alpha=0.833, rlxf=0.005, minC=0, maxC=1)

# SFS = vpm.DynamicSFS(vpm.Estr_fmm, vpm.sensorfunction; clippings=[vpm.clipping_backscatter],
#                                                         alpha=1.1, rlxf=0.5, minC=0, maxC=1)

# sensorfunction(args...; optargs...) = vpm.sensorfunction(args...; sensor=Lmbd->exp(-Lmbd^6),
#                                            Lambda=(lmbd, lmbdcrit) -> lmbd/lmbdcrit, optargs...)
# SFS = vpm.DynamicSFS(vpm.Estr_fmm, sensorfunction; clippings=[vpm.clipping_backscatter],
                                                        # alpha=1.1, rlxf=0.65, minC=0, maxC=1)

# -------------- SOLVER SETTINGS -------------------------------------------

solver = (
    # formulation   = vpm.cVPM,
    formulation   = vpm.rVPM,
    SFS           = SFS,
    kernel        = vpm.gaussianerf,
    # viscous       = vpm.Inviscid(),
    viscous       = vpm.CoreSpreading(-1, -1, vpm.zeta_fmm; beta=100.0, itmax=20, tol=1e-1,
                                                         iterror=true, verbose=false, debug=debug),
    transposed    = true,
    integration   = vpm.rungekutta3,
    # relaxation    = vpm.correctedpedrizzetti,
    relaxation    = vpm.pedrizzetti,
    UJ            = vpm.UJ_fmm,
    fmm           = vpm.FMM(; p=4, ncrit=50, theta=0.4, phi=0.3),
    # fmm           = vpm.FMM(; p=4, ncrit=50, theta=0.4, phi=0.5),
)

maxparticles = Int(15e6)               # Maximum number of particles

# --------------- DEFINE MONITORS ----------------------------------------------

axisline = [(O - 2*d*Oaxis[:, 3], O + 30*d*Oaxis[:, 3])]

transverselines = [
            (O - 10.00*d/2*Oaxis[:, 1] + zod*d*Oaxis[:, 3],
             O + 10.00*d/2*Oaxis[:, 1] + zod*d*Oaxis[:, 3])
            for zod in [0, 0.025, 0.05, 0.1, 0.28, 0.5, 1.12, 2.66, 4.48, 10, 15, 20, 25, 30, 1, 2, 3, 4, 5]
        ]

lines = vcat(axisline, transverselines)

probelines = generate_probelines(lines; nprobes=1000,
                                    save_path=save_path, fname_pref=run_name)

function monitor_enstrophy(args...; optargs...)
    bool1 = vpm.monitor_enstrophy_Gamma2(args...; save_path=save_path, suff="enstrophy-Gamma2.log", optargs...)
    bool2 = vpm.monitor_enstrophy_Gammaomega(args...; save_path=save_path, suff="enstrophy-Gammaomega.log", optargs...)
    return bool1 || bool2
end

remove_Gammapeak = nothing
function remove_particles(pfield, args...; optargs...)

    global remove_particles_weak
    global remove_Gammapeak

    if remove_particles_weak

        # Find maximum strength at the boundary condition
        if remove_Gammapeak == nothing
            remove_Gammapeak = 0
            for P in vpm.iterate(pfield; include_static=true)
                # NOTE: Out of laziness, this is hardcode to work only on jet at z=0 and aligned with z
                if P.static[1] && P.X[3]<0 && abs(P.X[3])<=0.1*d
                    if norm(P.Gamma) > remove_Gammapeak
                        remove_Gammapeak = norm(P.Gamma)
                    end
                end
            end
        end

        # Remove weak particles
        for pi in vpm.get_np(pfield):-1:1
            P = vpm.get_particle(pfield, pi)
            if P.static[1]==false && norm(P.Gamma) <= 0.05*remove_Gammapeak
                vpm.remove_particle(pfield, pi)
            end
        end

    end

    return false
end

runtime_functions = [
                        monitor_enstrophy,
                        probelines,
                        (args...; optargs...) -> vpm.monitor_Cd(args...; save_path=save_path, optargs...),
                        remove_particles
                    ]


# -------------- SIMULATION ------------------------------------------------
println("\nRunning simulation...")

run_roundjet_simulation(    # ------- SIMULATION PARAMETERS --------
                            d, U1, U2;
                            U2angle=U2angle,
                            O=O,
                            Oaxis=Oaxis,
                            thetaod=thetaod,
                            Re=Re,
                            # ------- SOLVER OPTIONS ---------------
                            maxparticles=maxparticles, pfieldargs=solver,
                            steps_per_d=steps_per_d,
                            d_travel_tot=d_travel_tot,
                            dxotheta=dxotheta,
                            overlap=overlap,
                            minWfraction=minWfraction,
                            max_zsigma=max_zsigma,
                            rbf=rbf,
                            rbf_optargs=rbf_optargs,
                            rbf_Wf=rbf_Wf,
                            rbf_fullcylinder=rbf_fullcylinder,
                            keep_rbfparticles=keep_rbfparticles,
                            keep_rbffullcylinder=keep_rbffullcylinder,
                            restart_file=restart_file,
                            restart_sigma=restart_sigma,
                            # ------- OUTPUT OPTIONS ---------------
                            verbose=verbose,
                            v_lvl=v_lvl,
                            runtime_functions=runtime_functions,
                            run_name=run_name,
                            verbose_nsteps=25,
                            save_code=@__FILE__
                            )
nothing

[Liuhaosheng1]

Hello author, thank you very much for your answer. I have another question. I pasted the above code into a file and ran it, but an error occurred at the beginning of the run. Secondly, is the circular jet validation case an independent model case, or does it require first calculating the hovering rotor and then importing the calculation results. Put the H5 file into this case to run it

[Liuhaosheng1]

I want to replicate the example you ran in your video using the validation case for round jet in FLOWVPM. I've tried running the roundjet.jl file, but I'm not getting the same results as in the video. Could there be an issue with the file I'm running?