|
|
@@ -1,3292 +0,0 @@
|
|
|
-"""
|
|
|
-
|
|
|
-Main function of Flower.jl code to run a simulation
|
|
|
-
|
|
|
-"""
|
|
|
-function run_forward!(
|
|
|
- num::Numerical{Float64, Int64},
|
|
|
- grid_p::Mesh{Flower.GridCC, Float64, Int64},
|
|
|
- grid_u::Mesh{Flower.GridFCx, Float64, Int64},
|
|
|
- grid_v::Mesh{Flower.GridFCy, Float64, Int64},
|
|
|
- op::DiscreteOperators{Float64, Int64},
|
|
|
- phS::Union{Phase{Float64},Nothing},
|
|
|
- phL::Phase{Float64};
|
|
|
- periodic_x::Bool = false,
|
|
|
- periodic_y::Bool = false,
|
|
|
- BC_TS = Boundaries(),
|
|
|
- BC_TL = Boundaries(),
|
|
|
- BC_pS = Boundaries(),
|
|
|
- BC_pL = Boundaries(),
|
|
|
- BC_uS = BoundariesInt(),
|
|
|
- BC_uL = BoundariesInt(),
|
|
|
- BC_vS = BoundariesInt(),
|
|
|
- BC_vL = BoundariesInt(),
|
|
|
- BC_u = Boundaries(),
|
|
|
- BC_trans_scal = Vector{BoundariesInt}(),
|
|
|
- BC_phi_ele = BoundariesInt(),
|
|
|
- BC_int::Vector{<:BoundaryCondition} = [WallNoSlip()],
|
|
|
- time_scheme::TemporalIntegration = CN,
|
|
|
- ls_scheme::LevelsetDiscretization = weno5,
|
|
|
- auto_reinit::Int64 = 0,
|
|
|
- heat::Bool = false,
|
|
|
- heat_convection::Bool = false,
|
|
|
- heat_liquid_phase::Bool = false,
|
|
|
- heat_solid_phase::Bool = false,
|
|
|
- navier_stokes ::Bool= false,
|
|
|
- ns_advection::Bool = false,
|
|
|
- ns_liquid_phase::Bool = false,
|
|
|
- ns_solid_phase::Bool = false,
|
|
|
- hill::Bool = false,
|
|
|
- Vmean::Bool = false,
|
|
|
- levelset::Bool = true,
|
|
|
- analytical::Bool = false,
|
|
|
- verbose::Bool = false,
|
|
|
- show_every::Int64 = 100,
|
|
|
- save_radius::Bool = false,
|
|
|
- adaptative_t::Bool = false,
|
|
|
- breakup::Int64 = 0,
|
|
|
- Ra::Float64 = 0.0,
|
|
|
- electrolysis::Bool = false,
|
|
|
- electrolysis_convection::Bool = false,
|
|
|
- electrolysis_liquid_phase::Bool = false,
|
|
|
- electrolysis_solid_phase::Bool = false,
|
|
|
- electrolysis_phase_change_case::String = "Khalighi",
|
|
|
- imposed_velocity::String = "none",
|
|
|
- adapt_timestep_mode::Int64 = 0,
|
|
|
- non_dimensionalize::Int64=1,
|
|
|
- mode_2d::Int64=0,
|
|
|
- test_laplacian::Bool = false,
|
|
|
- )
|
|
|
-
|
|
|
- #region Initialize simulation parameters
|
|
|
-
|
|
|
- λ = 1 #for Stefan velocity
|
|
|
- speed = 0.0
|
|
|
-
|
|
|
- radius_pdi = [0.0] #radius for pdi, list so that value is mutable
|
|
|
-
|
|
|
- if num.time > num.nucleation_time && num.phase_change_method >0 #TODO more precisely no mass transfer but velocity
|
|
|
- num.phase_change_currently_activated = 1
|
|
|
- end
|
|
|
-
|
|
|
- if num.mu1 == 0.0 && num.mu2 == 0.0 && num.mu_one_fluid_average !=0
|
|
|
- @error ("Resetting num.mu_one_fluid_average to 0")
|
|
|
- num.mu_one_fluid_average = 0
|
|
|
- end
|
|
|
-
|
|
|
- # if (num.one_fluid_model == 1 && num.solve_Navier_Stokes_liquid_phase == 1)
|
|
|
- # if CL BC
|
|
|
- # # grid_u.LS[1].cl
|
|
|
- # @error("BC CL with one fluid")
|
|
|
- # end
|
|
|
- # end
|
|
|
-
|
|
|
- #TODO Re
|
|
|
- #TODO homogenize set_poisson... functions which use -b - in the original implementation and not +b + for the Robin BC
|
|
|
-
|
|
|
- # Defining epsilon length and volume to handle special cells
|
|
|
- if num.epsilon_mode == 1 || num.epsilon_mode ==2
|
|
|
- num.epsilon_dist = eps(0.01) * num.Δ
|
|
|
- num.epsilon_vol = (eps(0.01)*num.Δ)^2
|
|
|
- num.epsilon_dist_mass_transfer_rate = num.Δ / 10 #TODO improve ex crit vol cf num.epsilon_volume_fraction_phase_change
|
|
|
- #TODO kill dead cells
|
|
|
- #TODO 1e-...
|
|
|
- end
|
|
|
-
|
|
|
- electrode_definition_function = (num.electrolysis_reaction_symb === :Butler_no_concentration) ? vecb_L : vecb_B
|
|
|
- # Temp = heat ? T : num.temperature0
|
|
|
-
|
|
|
-
|
|
|
- # if monophasic
|
|
|
- # if length(BC_int) != num.nLS
|
|
|
- # @error ("You have to specify $(num.nLS) boundary conditions.")
|
|
|
- # return nothing
|
|
|
- # end
|
|
|
-
|
|
|
- num.status = 0 # used to stop the simulation
|
|
|
- num.current_iter = 0 #1
|
|
|
- num.time = 0.
|
|
|
- free_surface = false
|
|
|
- stefan = false
|
|
|
- navier = false
|
|
|
- ls_advection = true
|
|
|
- if any(is_fs, BC_int)
|
|
|
- free_surface = true
|
|
|
- end
|
|
|
- if any(is_stefan, BC_int)
|
|
|
- stefan = true
|
|
|
- end
|
|
|
- if any(is_navier_cl, BC_int) || any(is_navier, BC_int)
|
|
|
- navier = true
|
|
|
- end
|
|
|
-
|
|
|
- if num.nNavier > 1
|
|
|
- @warn ("When using more than 1 Navier BC, the interfaces shouldn't cross")
|
|
|
- end
|
|
|
-
|
|
|
-
|
|
|
-
|
|
|
- if free_surface && stefan
|
|
|
- @error ("Cannot advect the levelset using both free-surface and stefan condition.")
|
|
|
- return nothing
|
|
|
- elseif free_surface || stefan || (num.one_fluid_model == 1 && num.solve_Navier_Stokes_liquid_phase == 1) && (num.activate_interface != 0) #|| electrolysis_phase_change_case !="none"
|
|
|
- # num.activate_interface != 0 : no interface, do not advect (monophasic)
|
|
|
- advection = true
|
|
|
- else
|
|
|
- advection = false
|
|
|
- end
|
|
|
-
|
|
|
- if num.solve_Navier_Stokes_liquid_phase == 0
|
|
|
- advection = false
|
|
|
- electrolysis_advection = false
|
|
|
-
|
|
|
- else
|
|
|
- if electrolysis
|
|
|
- electrolysis_advection = true
|
|
|
- end
|
|
|
- end
|
|
|
-
|
|
|
- if num.advection_LS_mode == 16 || num.advection_LS_mode == 16 #test Cipriano 2024 's method
|
|
|
- extend_liquid_velocity = true
|
|
|
- else
|
|
|
- extend_liquid_velocity = false
|
|
|
- end
|
|
|
-
|
|
|
- # The count threshold shouldn't be smaller than 2
|
|
|
- count_limit_breakup = 6
|
|
|
-
|
|
|
- if num.verbosity > 0
|
|
|
- printstyled(color=:green, @sprintf "\n num.CFL : %.2e dt : %.2e\n" num.CFL num.timestep_n)
|
|
|
- end
|
|
|
-
|
|
|
- #compute initial curvature (TODO other interfaces than circle)
|
|
|
- num.mean_curvature = 1.0/num.R
|
|
|
- num.current_radius = num.R
|
|
|
-
|
|
|
- # if num.surface_tension == 0
|
|
|
- # compute_surface_tension_VOF!(num,grid_p, grid_u, grid_v, opC_p, opC_u, opC_v,
|
|
|
- # volume_fraction,levelset_one_fluid,volumic_surface_tension_u,volumic_surface_tension_v,tmp_vec_p,tmp_vec_p0)
|
|
|
- # elseif num.surface_tension == 1
|
|
|
- # compute_surface_tension_LS!(num,grid_p, grid_u, grid_v, opC_p, opC_u, opC_v,
|
|
|
- # volume_fraction,levelset_one_fluid,volumic_surface_tension_u,volumic_surface_tension_v,tmp_vec_p,tmp_vec_p0,
|
|
|
- # levelset_1D, levelset_heavyside_2D, normal_and_dirac_u, normal_and_dirac_v,
|
|
|
- # normal_u, normal_v, curvature_u, curvature_v)
|
|
|
- # end
|
|
|
-
|
|
|
-
|
|
|
- num.timestep_n = adapt_timestep!(num, phL, phS, grid_u, grid_v,adapt_timestep_mode)
|
|
|
-
|
|
|
- pres_free_surfaceS = 0.0
|
|
|
- # pres_free_surfaceL = 0.0
|
|
|
- pres_free_surfaceL = num.pres0
|
|
|
-
|
|
|
- if occursin("levelset",electrolysis_phase_change_case)
|
|
|
- jump_mass_transfer_rateS = false
|
|
|
- jump_mass_transfer_rateL = true
|
|
|
- else
|
|
|
- jump_mass_transfer_rateS = false
|
|
|
- jump_mass_transfer_rateL = false
|
|
|
- end
|
|
|
-
|
|
|
- mass_transfer_rateS = 0.0
|
|
|
-
|
|
|
- iRe = 1.0 / num.Re
|
|
|
-
|
|
|
- # CFL_sc is not a CFL, it is supposed to be the timestep divided by the cell volume
|
|
|
- CFL_sc = num.timestep_n / num.Δ^2
|
|
|
-
|
|
|
- irho = 1.0
|
|
|
-
|
|
|
- num.mu_cin1 = num.mu1 / num.rho1
|
|
|
- num.mu_cin2 = num.mu2 / num.rho2
|
|
|
-
|
|
|
- if non_dimensionalize==0
|
|
|
- #force L=1 u=1
|
|
|
- Re=num.rho1/num.mu1
|
|
|
- iRe = 1.0/Re
|
|
|
- irho=1.0/num.rho1
|
|
|
- num.visc_coeff=iRe
|
|
|
- else
|
|
|
- num.visc_coeff=iRe
|
|
|
- Re = num.Re
|
|
|
- end
|
|
|
-
|
|
|
- if num.verbosity > 0 && num.solve_Navier_Stokes>0
|
|
|
- printstyled(color=:green, @sprintf "\n Re : %.2e %.2e\n" Re num.visc_coeff)
|
|
|
- printstyled(color=:magenta, @sprintf "\n CFL_sc : %.2e\n" CFL_sc)
|
|
|
- end
|
|
|
-
|
|
|
- # Electrolysis
|
|
|
- num.current_radius = 0.0
|
|
|
-
|
|
|
- # TODO kill_dead_cells! for [:,:,iscal]
|
|
|
- #region init number of moles
|
|
|
- if electrolysis
|
|
|
- if electrolysis_phase_change_case != "none"
|
|
|
- num.current_radius = num.R
|
|
|
-
|
|
|
- printstyled(color=:green, @sprintf "\n radius: %.2e \n" num.current_radius)
|
|
|
-
|
|
|
- p_liq= num.pres0 + mean(veci(phL.pD,grid_p,2)) #TODO here one bubble
|
|
|
- p_g=p_liq + 2 * num.σ / num.current_radius
|
|
|
-
|
|
|
- #TODO using num.temperature0
|
|
|
- if mode_2d==0
|
|
|
- num.nH2 = p_g * 4.0 / 3.0 * pi * num.current_radius ^ 3 / (num.temperature0 * num.Ru)
|
|
|
- elseif mode_2d == 1 #reference thickness for a cylinder
|
|
|
- num.nH2 = p_g * pi * num.current_radius ^ 2 * num.ref_thickness_2d / (num.temperature0 * num.Ru)
|
|
|
- elseif mode_2d==2 #mol/meter
|
|
|
- num.nH2=num.concentration0[num.index_phase_change]* pi * num.current_radius ^ 2
|
|
|
- elseif mode_2d==3 #mol/meter half circle
|
|
|
- num.nH2=1.0/2.0*num.concentration0[num.index_phase_change]* pi * num.current_radius ^ 2
|
|
|
- end
|
|
|
- # num.nH2 = 4.0/3.0 * pi * num.current_radius^3 * num.rho2 / num.MWH2
|
|
|
-
|
|
|
- printstyled(color=:green, @sprintf "\n Mole: %.2e \n" num.nH2)
|
|
|
-
|
|
|
- printstyled(color=:green, @sprintf "\n Mole test: %.2e %.2e\n" num.concentration0[num.index_phase_change]*4.0/3.0*pi*num.current_radius^3 p_g*4.0/3.0*pi*num.current_radius^3/(num.temperature0*num.Ru))
|
|
|
-
|
|
|
- end
|
|
|
- end # if electrolysis
|
|
|
- #endregion init number of moles
|
|
|
-
|
|
|
- #endregion Initialize simulation parameters
|
|
|
-
|
|
|
- local NB_indices;
|
|
|
-
|
|
|
- #region Allocations
|
|
|
-
|
|
|
- if num.solve_solid == 1
|
|
|
- local Cum1S = fzeros(grid_u)
|
|
|
- local Cvm1S = fzeros(grid_v)
|
|
|
- local Mm1_S
|
|
|
- local Mum1_S
|
|
|
- local Mvm1_S
|
|
|
- end
|
|
|
-
|
|
|
-
|
|
|
- local Cum1L = fzeros(grid_u)
|
|
|
- local Cvm1L = fzeros(grid_v)
|
|
|
- local Mm1_L
|
|
|
- local Mum1_L
|
|
|
- local Mvm1_L
|
|
|
-
|
|
|
- if extend_liquid_velocity #num.advection_LS_mode == 16
|
|
|
- # *_ext_vel : variables for the second system of Navier-Stokes equations
|
|
|
- local Cum1L_ext_vel = fzeros(grid_u)
|
|
|
- local Cvm1L_ext_vel = fzeros(grid_v)
|
|
|
- # local Mm1_L_ext_vel
|
|
|
- # local Mum1_L_ext_vel
|
|
|
- # local Mvm1_L_ext_vel
|
|
|
- p_ext_vel = zeros(grid_p)
|
|
|
- pD_ext_vel = fzeros(grid_p)
|
|
|
- phi_ext_vel = zeros(grid_p)
|
|
|
- u_ext_vel = zeros(grid_u)
|
|
|
- v_ext_vel= zeros(grid_v)
|
|
|
- u_predictionD_ext_vel= fzeros(grid_u)
|
|
|
- v_predictionD_ext_vel= fzeros(grid_v)
|
|
|
- uD_ext_vel= fzeros(grid_u)
|
|
|
- vD_ext_vel= fzeros(grid_v)
|
|
|
- u_prediction_ext_vel= zeros(grid_u)
|
|
|
- v_prediction_ext_vel= zeros(grid_v)
|
|
|
- uT_ext_vel = zeros(grid_p)
|
|
|
- pres_grad_x = fzeros(grid_u)
|
|
|
- pres_grad_y = fzeros(grid_v)
|
|
|
- else
|
|
|
- u_ext_vel = nothing
|
|
|
- v_ext_vel= nothing
|
|
|
-
|
|
|
- end
|
|
|
-
|
|
|
- θ_out = zeros(grid_p, 4)
|
|
|
- utmp = copy(grid_p.LS[1].u)
|
|
|
- rhs_LS = fzeros(grid_p)
|
|
|
-
|
|
|
- #vectors used reset at start of functions to limit allocations
|
|
|
- tmp_vec_u = zeros(grid_u)
|
|
|
- tmp_vec_v = zeros(grid_v)
|
|
|
- tmp_vec_u0 = zeros(grid_u)
|
|
|
- tmp_vec_v0 = zeros(grid_v)
|
|
|
-
|
|
|
- tmp_vec_u1 = zeros(grid_u)
|
|
|
- tmp_vec_v1 = zeros(grid_v)
|
|
|
-
|
|
|
- tmp_vec_p = zeros(grid_p)
|
|
|
- tmp_vec_p0 = zeros(grid_p)
|
|
|
- tmp_vec_p1 = zeros(grid_p)
|
|
|
-
|
|
|
- tmp_vec_1D = fnzeros(grid_p,num)
|
|
|
- tmp_vec_1D_2 = fnzeros(grid_p,num)
|
|
|
-
|
|
|
- # tmp_vec_1D_u = fnzeros(grid_p,num)
|
|
|
- tmp_vec_1D_v = fnzeros(grid_v,num)
|
|
|
- tmp_vec_1D_v0 = fnzeros(grid_v,num)
|
|
|
-
|
|
|
- tmp_vec_1D_p = fnzeros(grid_p,num)
|
|
|
- tmp_vec_1D_p0 = fnzeros(grid_p,num)
|
|
|
-
|
|
|
-
|
|
|
- # nb_gaz_acceptors = zeros(grid_p)
|
|
|
-
|
|
|
- #region electrolysis
|
|
|
- #TODO harmonic conductivity
|
|
|
- if electrolysis
|
|
|
- if num.nb_transported_scalars>1
|
|
|
- elec_cond = zeros(grid_p)
|
|
|
- elec_condD = fnzeros(grid_p,num)
|
|
|
-
|
|
|
- if heat
|
|
|
- elec_condD .= compute_ele_cond.(num.Faraday,num.diffusion_coeff[num.index_electrolyte],num.Ru, phL.TD, phL.trans_scalD[:,num.index_electrolyte])
|
|
|
- elec_cond .= reshape(vec1(elec_condD,grid_p),grid_p)
|
|
|
- # elec_cond .= compute_ele_cond.(num.Faraday,num.diffusion_coeff[num.index_electrolyte],num.Ru, phL.T, phL.trans_scal)
|
|
|
- # elec_cond = 2*num.Faraday^2 .*phL.trans_scal[:,:,2].*num.diffusion_coeff[2]./(num.Ru.*phL.T) #phL.T
|
|
|
- else
|
|
|
- elec_condD .= compute_ele_cond.(num.Faraday,num.diffusion_coeff[num.index_electrolyte],num.Ru, num.temperature0, phL.trans_scalD[:,num.index_electrolyte])
|
|
|
- elec_cond .= reshape(vec1(elec_condD,grid_p),grid_p)
|
|
|
- # elec_cond .= compute_ele_cond.(num.Faraday,num.diffusion_coeff[num.index_electrolyte],num.Ru, num.temperature0, phL.trans_scal)
|
|
|
- # elec_cond = 2*num.Faraday^2 .*phL.trans_scal[:,:,2].*num.diffusion_coeff[2]./(num.Ru*num.temperature0)
|
|
|
- end
|
|
|
- else
|
|
|
- elec_cond = ones(grid_p)
|
|
|
- elec_condD = fnones(grid_p,num)
|
|
|
- printstyled(color=:green, @sprintf "\n conductivity one")
|
|
|
- end
|
|
|
-
|
|
|
- # update_electrical_conductivity!(num,grid,elec_cond,elec_condD,heat;phL)
|
|
|
-
|
|
|
-
|
|
|
- # if num.electrolysis_reaction == "Butler_no_concentration"
|
|
|
- # i_butler = zeros(grid_p.ny) #left wall
|
|
|
- # elseif num.electrolysis_reaction == "fixed_current"
|
|
|
- # i_butler = zeros(grid_p.nx) #bottom wall
|
|
|
- # end
|
|
|
-
|
|
|
- printstyled(color=:red, @sprintf "\n Reaction")
|
|
|
-
|
|
|
- i_butler = Float64[] # empty vector, defined in scope
|
|
|
-
|
|
|
- print("\n electrolysis_reaction_symb ",num.electrolysis_reaction_symb)
|
|
|
-
|
|
|
- if num.electrolysis_reaction_symb === :none
|
|
|
-
|
|
|
- else
|
|
|
-
|
|
|
- size_BC_reaction = if num.electrolysis_reaction_symb === :Butler_no_concentration
|
|
|
- grid_p.ny
|
|
|
- elseif num.electrolysis_reaction_symb === :fixed_current
|
|
|
- grid_p.nx
|
|
|
- else
|
|
|
- error("Unknown electrolysis_reaction")
|
|
|
- end
|
|
|
-
|
|
|
- resize!(i_butler, size_BC_reaction)
|
|
|
- fill!(i_butler, 0.0)
|
|
|
-
|
|
|
- end
|
|
|
-
|
|
|
-
|
|
|
- end #electrolysis
|
|
|
- #endregion electrolysis
|
|
|
-
|
|
|
- if levelset
|
|
|
-
|
|
|
- # At every iteration, update_all_ls_data is called twice,
|
|
|
- # once inside run.jl and another one (if there's advection of the levelset) inside set_heat!.
|
|
|
- # The difference between both is a flag as last argument, inside run.jl is
|
|
|
- # implicitly defined as true and inside set_heat is false.
|
|
|
- # If you're calling your version of set_heat several times, then you're calling the version with the flag set to false, but for the convective term it has to be set to true, so that's why
|
|
|
- # The flag=true, the capacities are set for the convection, the flag=false they are set for the other operators
|
|
|
-
|
|
|
- NB_indices = update_all_ls_data(num, grid_p, grid_u, grid_v, BC_int, periodic_x, periodic_y)
|
|
|
-
|
|
|
- # printstyled(color=:red, @sprintf "\n levelset:\n")
|
|
|
- # println(grid_p.LS[1].geoL.dcap[1,1,:])
|
|
|
-
|
|
|
- if save_radius
|
|
|
- n_snaps = iszero(num.max_iterations%num.save_every) ? num.max_iterations÷num.save_every+1 : num.max_iterations÷num.save_every+2
|
|
|
- local radius = zeros(n_snaps)
|
|
|
- radius[1] = find_radius(grid_p, grid_p.LS[1])
|
|
|
- end
|
|
|
- if hill
|
|
|
- local radius = zeros(num.max_iterations+1)
|
|
|
- a = zeros(length(grid_p.LS[1].MIXED))
|
|
|
- for i in eachindex(grid_p.LS[1].MIXED)
|
|
|
- a[i] = grid_p.LS[1].geoL.projection[grid_p.LS[1].MIXED[i]].pos.y
|
|
|
- end
|
|
|
- radius[1] = mean(a)
|
|
|
- end
|
|
|
- elseif !levelset
|
|
|
- grid_p.LS[1].MIXED = [CartesianIndex(-1,-1)]
|
|
|
- grid_u.LS[1].MIXED = [CartesianIndex(-1,-1)]
|
|
|
- grid_v.LS[1].MIXED = [CartesianIndex(-1,-1)]
|
|
|
- end
|
|
|
-
|
|
|
- # if save_length
|
|
|
- # fwd.length[1] = arc_length2(grid_p.LS[1].geoS.projection, grid_p.LS[1].MIXED)
|
|
|
- # end
|
|
|
-
|
|
|
- #endregion
|
|
|
-
|
|
|
-
|
|
|
- #region Initialisation of bulk and interfacial values
|
|
|
- # Initialisation of bulk and interfacial values
|
|
|
-
|
|
|
- #TODO which grid_p.LS
|
|
|
-
|
|
|
- #TODO perio, intfc, ... check init_fields_2!
|
|
|
-
|
|
|
- #No scalar in "solid" phase
|
|
|
- # @views init_fields_multiple_levelsets!(num,phS.trans_scalD[:,iscal],phS.trans_scal[:,:,iscal],HS,BC_trans_scal[iscal],grid_p,num.concentration0[iscal])
|
|
|
- # @views phS.trans_scal[:,:,iscal] .= num.concentration0[iscal]
|
|
|
-
|
|
|
- # TODO reset zero
|
|
|
- tmp_vec_u .= 0.0
|
|
|
-
|
|
|
- if num.solve_solid == 1
|
|
|
- #from LS 1 to centroid grid_u.LS[end].geoS
|
|
|
- get_height!(grid_u.LS[1],grid_u.ind,grid_u.dx,grid_u.dy,grid_u.LS[end].geoS,tmp_vec_u) #here tmp_vec_u solid
|
|
|
-
|
|
|
- init_fields_multiple_levelsets!(num,phS.uD,phS.u,tmp_vec_u,BC_uS,grid_u,num.uD,"uS")
|
|
|
- # init_fields_multiple_levelsets!(num,phS.u_predictionD,phS.u,HSu,BC_uS,grid_u,num.uD)
|
|
|
- end
|
|
|
-
|
|
|
- get_height!(grid_u.LS[1],grid_u.ind,grid_u.dx,grid_u.dy,grid_u.LS[end].geoL,tmp_vec_u) #here tmp_vec_u liquid
|
|
|
-
|
|
|
- init_fields_multiple_levelsets!(num,phL.uD,phL.u,tmp_vec_u,BC_uL,grid_u,num.uD,"uL")
|
|
|
- # init_fields_multiple_levelsets!(num,phL.u_predictionD,phL.u,HLu,BC_uL,grid_u,num.uD)
|
|
|
-
|
|
|
-
|
|
|
- # TODO reset zero
|
|
|
- tmp_vec_v .= 0.0
|
|
|
- if num.solve_solid == 1
|
|
|
- get_height!(grid_v.LS[1],grid_v.ind,grid_v.dx,grid_v.dy,grid_v.LS[end].geoS,tmp_vec_v)
|
|
|
-
|
|
|
- init_fields_multiple_levelsets!(num,phS.vD,phS.v,tmp_vec_v,BC_vS,grid_v,num.vD,"vS")
|
|
|
- # init_fields_multiple_levelsets!(num,phS.v_predictionD,phS.v,HSv,BC_vS,grid_v,num.vD)
|
|
|
- end
|
|
|
-
|
|
|
- get_height!(grid_v.LS[1],grid_v.ind,grid_v.dx,grid_v.dy,grid_v.LS[end].geoL,tmp_vec_v)
|
|
|
-
|
|
|
- init_fields_multiple_levelsets!(num,phL.vD,phL.v,tmp_vec_v,BC_vL,grid_v,num.vD,"vL")
|
|
|
- # init_fields_multiple_levelsets!(num,phL.v_predictionD,phL.v,HLv,BC_vL,grid_v,num.vD)
|
|
|
-
|
|
|
-
|
|
|
- # TODO reset zero
|
|
|
- tmp_vec_p .= 0.0
|
|
|
-
|
|
|
- if num.solve_solid == 1
|
|
|
-
|
|
|
- get_height!(grid_p.LS[1],grid_p.ind,grid_p.dx,grid_p.dy,grid_p.LS[end].geoS,tmp_vec_p) #here tmp_vec_p solid
|
|
|
-
|
|
|
- init_fields_multiple_levelsets!(num,phS.pD,phS.p,tmp_vec_p,BC_pS,grid_p,num.pres_intfc,"pS")
|
|
|
-
|
|
|
- if heat
|
|
|
- init_fields_multiple_levelsets!(num,phS.TD,phS.T,tmp_vec_p,BC_TS,grid_p,num.θd,"TS")
|
|
|
- end
|
|
|
-
|
|
|
- #Electrolysis
|
|
|
- if electrolysis && num.electrical_potential > 0
|
|
|
- init_fields_multiple_levelsets!(num,phS.phi_eleD,phS.phi_ele,tmp_vec_p,BC_phi_ele,grid_p,num.phi_ele0,"phiS")
|
|
|
- end
|
|
|
- end
|
|
|
-
|
|
|
- get_height!(grid_p.LS[1],grid_p.ind,grid_p.dx,grid_p.dy,grid_p.LS[end].geoL,tmp_vec_p) #here tmp_vec_p liquid
|
|
|
-
|
|
|
- init_fields_multiple_levelsets!(num,phL.pD,phL.p,tmp_vec_p,BC_pL,grid_p,num.pres_intfc,"pL")
|
|
|
-
|
|
|
- if heat
|
|
|
- init_fields_multiple_levelsets!(num,phL.TD,phL.T,tmp_vec_p,BC_TL,grid_p,num.θd,"TL")
|
|
|
- end
|
|
|
-
|
|
|
- # Electrolysis
|
|
|
- if electrolysis
|
|
|
-
|
|
|
- printstyled(color=:green, @sprintf "\n Check %s %s %s %s %.2e %.2e %2i\n" heat heat_convection electrolysis electrolysis_convection num.timestep_n num.θd num.nb_transported_scalars)
|
|
|
-
|
|
|
- for iscal=1:num.nb_transported_scalars
|
|
|
- @views phL.trans_scal[:,:,iscal] .= num.concentration0[iscal]
|
|
|
- @views init_fields_multiple_levelsets!(num,phL.trans_scalD[:,iscal],phL.trans_scal[:,:,iscal],
|
|
|
- tmp_vec_p,BC_trans_scal[iscal],grid_p,num.concentration0[iscal],"scalL")
|
|
|
- # tmp_vec_p,BC_trans_scal[iscal],grid_p,num.concentration0[iscal],"testscalL")
|
|
|
-
|
|
|
- end
|
|
|
-
|
|
|
- if num.electrical_potential > 0 #TODO init phi =0 or with Neumann for BC of concentration? for now not done since "phiL" given in arg
|
|
|
- init_fields_multiple_levelsets!(num,phL.phi_eleD,phL.phi_ele,tmp_vec_p,BC_phi_ele,grid_p,num.phi_ele0,"phiL")
|
|
|
- if num.electrical_potential == 3
|
|
|
- vecb(phL.phi_eleD,grid_p) .= 0.0
|
|
|
- end
|
|
|
- end
|
|
|
- end
|
|
|
-
|
|
|
- if electrolysis
|
|
|
-
|
|
|
- printstyled(color=:green, @sprintf "\n Check init_fields_2!\n")
|
|
|
- # print_electrolysis_statistics(num,grid_p,phL)
|
|
|
-
|
|
|
- PDI_status = @ccall "libpdi".PDI_multi_expose("print_variables"::Cstring,
|
|
|
- "nstep"::Cstring, num.current_iter ::Ref{Clonglong}, PDI_OUT::Cint,
|
|
|
- "time"::Cstring, num.time::Ref{Cdouble}, PDI_OUT::Cint,
|
|
|
- "u_1D"::Cstring, phL.uD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "v_1D"::Cstring, phL.vD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "p_1D"::Cstring, phL.pD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "levelset_p"::Cstring, grid_p.LS[num.iLSpdi].u::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "levelset_u"::Cstring, grid_u.LS[num.iLSpdi].u::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "levelset_v"::Cstring, grid_v.LS[num.iLSpdi].u::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # "levelset_p_wall"::Cstring, LStable::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "trans_scal_1DT"::Cstring, phL.trans_scalD'::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "phi_ele_1D"::Cstring, phL.phi_eleD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # "i_current_x"::Cstring, Eus::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # "i_current_y"::Cstring, Evs::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # "velocity_x"::Cstring, us::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # "velocity_y"::Cstring, vs::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # "radius"::Cstring, current_radius::Ref{Cdouble}, PDI_OUT::Cint,
|
|
|
- # "intfc_vtx_num"::Cstring, intfc_vtx_num::Ref{Clonglong}, PDI_OUT::Cint,
|
|
|
- # "intfc_seg_num"::Cstring, intfc_seg_num::Ref{Clonglong}, PDI_OUT::Cint,
|
|
|
- # "intfc_vtx_x"::Cstring, intfc_vtx_x::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # "intfc_vtx_y"::Cstring, intfc_vtx_y::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # "intfc_vtx_field"::Cstring, intfc_vtx_field::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # "intfc_vtx_connectivities"::Cstring, intfc_vtx_connectivities::Ptr{Clonglong}, PDI_OUT::Cint,
|
|
|
- C_NULL::Ptr{Cvoid})::Cint
|
|
|
-
|
|
|
- if num.solve_solid == 1
|
|
|
-
|
|
|
- if (any(isnan, phL.uD) || any(isnan, phL.vD) || any(isnan, phL.TD) ||
|
|
|
- any(isnan, phS.uD) || any(isnan, phS.vD) || any(isnan, phS.TD) ||
|
|
|
- any(isnan, phL.trans_scalD) || any(isnan, phL.phi_eleD) ||
|
|
|
- norm(phL.u) > 1e8 || norm(phL.T) > 1e8 ||
|
|
|
- norm(phS.u) > 1e8 || norm(phS.T) > 1e8
|
|
|
- # norm(phL.trans_scal) > 1e8
|
|
|
- || norm(phL.phi_ele) > 1e8
|
|
|
- # ||
|
|
|
- # any(phL.trans_scal .<0)
|
|
|
- )
|
|
|
- println(@sprintf "\n CRASHED start \n")
|
|
|
-
|
|
|
- print("\n phL.uD: ",any(isnan, phL.uD) , "\n phL.vD: ",any(isnan, phL.vD) , "\n phL.TD: ",any(isnan, phL.TD) ,
|
|
|
- "\n phS.uD: ",any(isnan, phS.uD) , "\n phS.vD: ",any(isnan, phS.vD) , "\n phS.TD: ",any(isnan, phS.TD) ,
|
|
|
- "\n phL.trans_scalD: ",any(isnan, phL.trans_scalD) , "\n phL.phi_eleD: ",any(isnan, phL.phi_eleD) ,
|
|
|
- "\n phL.u: ",norm(phL.u) > 1e8 , "\n phS.u: ",norm(phS.u) > 1e8 , "\n phL.T: ",norm(phL.T) > 1e8 ,
|
|
|
- "\n phS.T: ",norm(phS.T) > 1e8 , "\n phL.trans_scal: ",norm(phL.trans_scal) > 1e8 ,
|
|
|
- "\n phL.phi_ele: ",norm(phL.phi_ele) > 1e8,"\n any(phL.trans_scal .<0): ", any(phL.trans_scal .<0))
|
|
|
-
|
|
|
- num.status = 1
|
|
|
- return num.current_iter
|
|
|
-
|
|
|
- end
|
|
|
- else
|
|
|
-
|
|
|
- if (any(isnan, phL.uD) || any(isnan, phL.vD) || any(isnan, phL.TD) ||
|
|
|
- any(isnan, phL.trans_scalD) || any(isnan, phL.phi_eleD) ||
|
|
|
- norm(phL.u) > 1e8 || norm(phL.T) > 1e8 ||
|
|
|
- # norm(phL.trans_scal) > 1e8 ||
|
|
|
- norm(phL.phi_ele) > 1e8
|
|
|
- # || any(phL.trans_scal .<0)
|
|
|
- )
|
|
|
- println(@sprintf "\n CRASHED start \n")
|
|
|
-
|
|
|
- # println(@sprintf "\n CRASHED after %d iterations \n" num.current_iter)
|
|
|
-
|
|
|
- print("\n phL.uD: ",any(isnan, phL.uD) , "\n phL.vD: ",any(isnan, phL.vD) , "\n phL.TD: ",any(isnan, phL.TD) ,
|
|
|
- "\n phL.trans_scalD: ",any(isnan, phL.trans_scalD) , "\n phL.phi_eleD: ",any(isnan, phL.phi_eleD) ,
|
|
|
- "\n phL.u: ",norm(phL.u) > 1e8, "\n phL.T: ",norm(phL.T) > 1e8 ,
|
|
|
- "\n phL.trans_scal: ",norm(phL.trans_scal) > 1e8 ,
|
|
|
- "\n phL.phi_ele: ",norm(phL.phi_ele) > 1e8,"\n any(phL.trans_scal .<0): ", any(phL.trans_scal .<0))
|
|
|
-
|
|
|
- num.status = 1
|
|
|
- return num.current_iter
|
|
|
-
|
|
|
- end
|
|
|
-
|
|
|
- end
|
|
|
-
|
|
|
- interpolate_staggered_u_v_to_scalar_grid_one_fluid_or_one_phase!(num,grid_p,grid_u,grid_v,phL.u,phL.v,tmp_vec_p,tmp_vec_p0)
|
|
|
-
|
|
|
- PDI_status = @ccall "libpdi".PDI_multi_expose("write_data"::Cstring,
|
|
|
- "nstep"::Cstring, num.current_iter::Ref{Clonglong}, PDI_OUT::Cint,
|
|
|
- "time"::Cstring, num.time::Ref{Cdouble}, PDI_OUT::Cint,
|
|
|
- "timestep"::Cstring, num.timestep_n::Ref{Cdouble}, PDI_OUT::Cint,
|
|
|
- "nx"::Cstring, grid_p.nx::Ref{Clonglong}, PDI_OUT::Cint,
|
|
|
- "ny"::Cstring, grid_p.ny::Ref{Clonglong}, PDI_OUT::Cint,
|
|
|
- "u_1D"::Cstring, phL.uD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "v_1D"::Cstring, phL.vD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "p_1D"::Cstring, phL.pD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "levelset_p"::Cstring, grid_p.LS[num.iLSpdi].u::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "levelset_u"::Cstring, grid_u.LS[num.iLSpdi].u::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "levelset_v"::Cstring, grid_v.LS[num.iLSpdi].u::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "trans_scal_1DT"::Cstring, phL.trans_scalD'::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "phi_ele_1D"::Cstring, phL.phi_eleD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # "i_current_x"::Cstring, tmp_vec_p::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # "i_current_y"::Cstring, tmp_vec_p0::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # "normal_x"::Cstring, normal_x::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # "normal_y"::Cstring, normal_y::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # grid_u.LS[iLS].α
|
|
|
- # "normal_angle"::Cstring, grid_p.LS[num.iLSpdi].α::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "velocity_x"::Cstring, tmp_vec_p::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "velocity_y"::Cstring, tmp_vec_p0::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "radius"::Cstring, num.current_radius::Ref{Cdouble}, PDI_OUT::Cint,
|
|
|
- # "intfc_vtx_num"::Cstring, intfc_vtx_num::Ref{Clonglong}, PDI_OUT::Cint,
|
|
|
- # "intfc_seg_num"::Cstring, intfc_seg_num::Ref{Clonglong}, PDI_OUT::Cint,
|
|
|
- # "intfc_vtx_x"::Cstring, intfc_vtx_x::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # "intfc_vtx_y"::Cstring, intfc_vtx_y::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # "intfc_vtx_field"::Cstring, intfc_vtx_field::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # "intfc_vtx_connectivities"::Cstring, intfc_vtx_connectivities::Ptr{Clonglong}, PDI_OUT::Cint,
|
|
|
- C_NULL::Ptr{Cvoid})::Cint
|
|
|
-
|
|
|
-
|
|
|
-
|
|
|
-
|
|
|
- end
|
|
|
-
|
|
|
- if num.solve_solid == 1
|
|
|
- kill_dead_cells!(phS.T, grid_p, grid_p.LS[end].geoS)
|
|
|
- end
|
|
|
- kill_dead_cells!(phL.T, grid_p, grid_p.LS[end].geoL)
|
|
|
-
|
|
|
-
|
|
|
- #TODO check timestep coefficients num.n-1
|
|
|
-
|
|
|
- #Electrolysis
|
|
|
- # TODO kill_dead_cells! for [:,:,iscal]
|
|
|
- if electrolysis
|
|
|
- for iscal=1:num.nb_transported_scalars
|
|
|
- # @views kill_dead_cells!(phS.trans_scal[:,:,iscal], grid_p, grid_p.LS[end].geoS) #TODO
|
|
|
- # @views kill_dead_cells!(phL.trans_scal[:,:,iscal], grid_p, grid_p.LS[end].geoL)
|
|
|
- # @views kill_dead_cells_val!(phS.trans_scal[:,:,iscal], grid_p, grid_p.LS[end].geoS) #TODO
|
|
|
- # @views kill_dead_cells_val!(phL.trans_scal[:,:,iscal], grid_p, grid_p.LS[end].geoL,num.concentration0[iscal])
|
|
|
- # @views kill_dead_cells_val!(phL.trans_scal[:,:,iscal], grid_p, grid_p.LS[end].geoL,0.0)
|
|
|
-
|
|
|
- if num.kill_dead_cells == 0
|
|
|
- @views kill_dead_cells_val!(phL.trans_scal[:,:,iscal], grid_p, grid_p.LS[end].geoL,0.0)
|
|
|
- else
|
|
|
- @views kill_dead_cells_val!(phL.trans_scal[:,:,iscal], grid_p, grid_p.LS[end].geoL,num.concentration0[iscal])
|
|
|
- end
|
|
|
-
|
|
|
- @views veci(phL.trans_scalD[:,iscal],grid_p,1) .= vec(phL.trans_scal[:,:,iscal])
|
|
|
-
|
|
|
-
|
|
|
- end
|
|
|
- end #if electrolysis
|
|
|
-
|
|
|
- #endregion
|
|
|
-
|
|
|
-
|
|
|
-
|
|
|
-
|
|
|
-
|
|
|
-
|
|
|
- # print("\n vecb_L(elec_condD, grid_p) after kill \n ", vecb_L(phL.trans_scalD[:,2], grid_p) )
|
|
|
-
|
|
|
-
|
|
|
- # if num.nb_transported_scalars>0
|
|
|
- # printstyled(color=:green, @sprintf "\n after kill \n")
|
|
|
- # print_electrolysis_statistics(num,phL)
|
|
|
- # printstyled(color=:green, @sprintf "\n average T %s\n" average!(phL.T, grid_p, grid_p.LS[1].geoL, num))
|
|
|
- # end
|
|
|
-
|
|
|
- #region Allocations
|
|
|
-
|
|
|
- # Set matrices/operators
|
|
|
- if is_Forward_Euler(time_scheme) || is_Crank_Nicolson(time_scheme)
|
|
|
- NB_indices = update_all_ls_data(num, grid_p, grid_u, grid_v, BC_int, periodic_x, periodic_y, false)
|
|
|
-
|
|
|
- # printstyled(color=:red, @sprintf "\n levelset 2:\n")
|
|
|
- # println(grid_p.LS[1].geoL.dcap[1,1,:])
|
|
|
-
|
|
|
- if navier_stokes || heat || electrolysis
|
|
|
- if num.solve_solid == 1
|
|
|
- geoS = [grid_p.LS[iLS].geoS for iLS in 1:num._nLS]
|
|
|
- geo_uS = [grid_u.LS[iLS].geoS for iLS in 1:num._nLS]
|
|
|
- geo_vS = [grid_v.LS[iLS].geoS for iLS in 1:num._nLS]
|
|
|
- Lpm1_S, bc_Lpm1_S, bc_Lpm1_b_S, Lum1_S, bc_Lum1_S, bc_Lum1_b_S, Lvm1_S, bc_Lvm1_S, bc_Lvm1_b_S = set_matrices!(
|
|
|
- num, grid_p, geoS, grid_u, geo_uS, grid_v, geo_vS,
|
|
|
- op.opC_pS, op.opC_uS, op.opC_vS, periodic_x, periodic_y
|
|
|
- )
|
|
|
- end
|
|
|
-
|
|
|
- geoL = [grid_p.LS[iLS].geoL for iLS in 1:num._nLS]
|
|
|
- geo_uL = [grid_u.LS[iLS].geoL for iLS in 1:num._nLS]
|
|
|
- geo_vL = [grid_v.LS[iLS].geoL for iLS in 1:num._nLS]
|
|
|
- Lpm1_L, bc_Lpm1_L, bc_Lpm1_b_L, Lum1_L, bc_Lum1_L, bc_Lum1_b_L, Lvm1_L, bc_Lvm1_L, bc_Lvm1_b_L = set_matrices!(
|
|
|
- num, grid_p, geoL, grid_u, geo_uL, grid_v, geo_vL,
|
|
|
- op.opC_pL, op.opC_uL, op.opC_vL, periodic_x, periodic_y
|
|
|
- )
|
|
|
- end
|
|
|
-
|
|
|
- # Contains volumes
|
|
|
- if num.solve_solid == 1
|
|
|
- Mm1_S = copy(op.opC_pS.M)
|
|
|
- Mum1_S = copy(op.opC_uS.M)
|
|
|
- Mvm1_S = copy(op.opC_vS.M)
|
|
|
- end
|
|
|
-
|
|
|
- # M defined in set_cutcell_matrices!
|
|
|
- Mm1_L = copy(op.opC_pL.M)
|
|
|
- Mum1_L = copy(op.opC_uL.M)
|
|
|
- Mvm1_L = copy(op.opC_vL.M)
|
|
|
-
|
|
|
- if num.one_fluid_model == 1
|
|
|
- @error("\n Mum1L needs to be redefined for one-fluid ")
|
|
|
-
|
|
|
- end
|
|
|
-
|
|
|
-
|
|
|
- if navier_stokes || heat || electrolysis
|
|
|
-
|
|
|
- #Allocations
|
|
|
- #TODO pre-allocate at start to save up allocations
|
|
|
- #TODO optimize allocations
|
|
|
-
|
|
|
- # #Preallocate for mass flux computations
|
|
|
- # if electrolysis_phase_change_case != "None"
|
|
|
- # mass_transfer_rate_vec1 = fzeros(grid_p)
|
|
|
- # mass_transfer_rate_vecb = fzeros(grid_p)
|
|
|
- # mass_transfer_rate_veci = fzeros(grid_p)
|
|
|
- # mass_transfer_rate = zeros(grid_p)
|
|
|
- # else
|
|
|
- # mass_transfer_rate = 0.0
|
|
|
- # end
|
|
|
- mass_transfer_rate_vec1 = fzeros(grid_p)
|
|
|
- mass_transfer_rate_vecb = fzeros(grid_p)
|
|
|
- mass_transfer_rate_veci = fzeros(grid_p)
|
|
|
- mass_transfer_rate = zeros(grid_p)
|
|
|
-
|
|
|
- mass_transfer_rate_redistributed = zeros(grid_p)
|
|
|
- nb_gaz_acceptors = zeros(Int64,grid_p.ny,grid_p.nx)
|
|
|
- # nb_gaz_acceptors = zeros(grid_p)
|
|
|
-
|
|
|
- #Allocations for scalar grid_p
|
|
|
- ni = grid_p.nx * grid_p.ny
|
|
|
- nb = 2 * grid_p.nx + 2 * grid_p.ny
|
|
|
- nt = (num.nLS + 1) * ni + nb
|
|
|
-
|
|
|
- a1_p = spdiagm(ni,ni,zeros(ni))
|
|
|
-
|
|
|
- Ascal = spzeros(nt, nt)
|
|
|
- Bscal = spzeros(nt, nt)
|
|
|
- rhs_scal = fnzeros(grid_p, num)
|
|
|
- F_residual = fnzeros(grid_p, num)
|
|
|
-
|
|
|
-
|
|
|
- all_CUTCT = zeros(grid_p.ny * grid_p.nx, num.nb_transported_scalars)
|
|
|
-
|
|
|
- if num.one_fluid_model == 0
|
|
|
- if num.solve_solid == 1
|
|
|
- AϕS = spzeros(nt, nt)
|
|
|
- end
|
|
|
- AϕL = spzeros(nt, nt)
|
|
|
-
|
|
|
- else
|
|
|
- nt_pressure = ni + nb
|
|
|
- AϕL = spzeros(nt_pressure, nt_pressure)
|
|
|
- rhs_phi = zeros(nt_pressure)
|
|
|
- end
|
|
|
-
|
|
|
- if electrolysis
|
|
|
- # Aphi_eleL = spzeros(nt, nt)
|
|
|
-
|
|
|
- # coeffDu = zeros(grid_u)
|
|
|
- # coeffDv = zeros(grid_v)
|
|
|
- coeffDx_interface = zeros(grid_u)
|
|
|
- coeffDy_interface = zeros(grid_v)
|
|
|
- end
|
|
|
-
|
|
|
- if num.solve_solid == 1
|
|
|
- ATS = spzeros(nt, nt)
|
|
|
- BTS = spzeros(nt, nt)
|
|
|
- end
|
|
|
- ATL = spzeros(nt, nt)
|
|
|
- BTL = spzeros(nt, nt)
|
|
|
-
|
|
|
- ni = grid_u.nx * grid_u.ny
|
|
|
- nb = 2 * grid_u.nx + 2 * grid_u.ny
|
|
|
- nt = (num.nLS + 1) * ni + nb
|
|
|
- if num.solve_solid == 1
|
|
|
- AuS = spzeros(nt, nt)
|
|
|
- BuS = spzeros(nt, nt)
|
|
|
- end
|
|
|
- AuL = spzeros(nt, nt)
|
|
|
- BuL = spzeros(nt, nt)
|
|
|
-
|
|
|
- ni = grid_v.nx * grid_v.ny
|
|
|
- nb = 2 * grid_v.nx + 2 * grid_v.ny
|
|
|
- nt = (num.nLS + 1) * ni + nb
|
|
|
- if ns_solid_phase
|
|
|
- AvS = spzeros(nt, nt)
|
|
|
- BvS = spzeros(nt, nt)
|
|
|
- end
|
|
|
- AvL = spzeros(nt, nt)
|
|
|
- BvL = spzeros(nt, nt)
|
|
|
-
|
|
|
- #region Allocate for Navier case
|
|
|
- ni = grid_u.nx * grid_u.ny + grid_v.nx * grid_v.ny
|
|
|
- nb = 2 * grid_u.nx + 2 * grid_u.ny + 2 * grid_v.nx + 2 * grid_v.ny
|
|
|
- nt = (num.nLS - num.nNavier + 1) * ni + num.nNavier * grid_p.nx * grid_p.ny + nb
|
|
|
- #dev was done with nNavier == ?
|
|
|
-
|
|
|
- #when no Navier: nt = (num.nLS + 1) * ni + nb
|
|
|
-
|
|
|
-
|
|
|
- # if ((num.one_fluid_model == 1 && num.solve_Navier_Stokes_liquid_phase == 1) && num.pressure_velocity_coupling != 0)
|
|
|
- # @error("\nCoupled pressure velocity + one-fluid model error")
|
|
|
- # return
|
|
|
- # end
|
|
|
-
|
|
|
-
|
|
|
- if (num.one_fluid_model == 1 && num.solve_Navier_Stokes_liquid_phase == 1)
|
|
|
- rho_one_fluid = zeros(grid_p)
|
|
|
- mu_one_fluid = zeros(grid_p)
|
|
|
- volume_fraction = zeros(grid_p)
|
|
|
- interface_length = zeros(grid_p)
|
|
|
- levelset_one_fluid = zeros(grid_p)
|
|
|
-
|
|
|
- rho_one_fluid_u = zeros(grid_u)
|
|
|
- # mu_one_fluid_u = zeros(grid_u)
|
|
|
-
|
|
|
- rho_one_fluid_v = zeros(grid_v)
|
|
|
- # mu_one_fluid_v = zeros(grid_v)
|
|
|
-
|
|
|
- volumic_surface_tension_u = zeros(grid_u)
|
|
|
- volumic_surface_tension_v = zeros(grid_v)
|
|
|
-
|
|
|
- # Pre-allocate arrays
|
|
|
- levelset_1D = fnzeros(grid_p, num)
|
|
|
- levelset_heavyside_2D = zeros(grid_p)
|
|
|
-
|
|
|
- convection_u = fzeros(grid_u)
|
|
|
- convection_v = fzeros(grid_v)
|
|
|
-
|
|
|
- viscosity_coeff_for_du_dx = zeros(grid_u.ny, grid_u.nx+1)
|
|
|
- viscosity_coeff_for_du_dy = zeros(grid_u.ny+1, grid_u.nx)
|
|
|
- viscosity_coeff_for_dv_dx = zeros(grid_v.ny, grid_v.nx+1)
|
|
|
- viscosity_coeff_for_dv_dy = zeros(grid_v.ny+1, grid_v.nx)
|
|
|
-
|
|
|
- velocity_and_BC_convection_u_x = zeros(grid_u) #Du_x
|
|
|
- velocity_and_BC_convection_u_y = zeros(grid_u)
|
|
|
- velocity_and_BC_convection_v_x = zeros(grid_v)
|
|
|
- velocity_and_BC_convection_v_y = zeros(grid_v)
|
|
|
-
|
|
|
- else
|
|
|
- volume_fraction = nothing
|
|
|
- interface_length = nothing
|
|
|
- end
|
|
|
-
|
|
|
- if (num.pressure_velocity_coupling == 0 && num.one_fluid_model == 0)
|
|
|
-
|
|
|
- if ns_solid_phase
|
|
|
- AuvS = spzeros(nt, nt)
|
|
|
- BuvS = spzeros(nt, nt)
|
|
|
- end
|
|
|
-
|
|
|
- AuvL = spzeros(nt, nt)
|
|
|
- BuvL = spzeros(nt, nt)
|
|
|
-
|
|
|
- else
|
|
|
- # We solve for u, v with borders
|
|
|
- # dev was done with nNavier == ?
|
|
|
-
|
|
|
- ni_p = grid_p.nx * grid_p.ny
|
|
|
- nb_p = 2 * grid_p.nx + 2 * grid_p.ny
|
|
|
-
|
|
|
- ni_u = grid_u.nx * grid_u.ny
|
|
|
- nb_u = 2 * grid_u.nx + 2 * grid_u.ny
|
|
|
-
|
|
|
- ni_v = grid_v.nx * grid_v.ny
|
|
|
- nb_v = 2 * grid_v.nx + 2 * grid_v.ny
|
|
|
-
|
|
|
- ni_uv = ni_u + ni_v
|
|
|
- nb_uv = nb_u + nb_v
|
|
|
-
|
|
|
- if num.pressure_velocity_coupling == 0
|
|
|
- if (num.one_fluid_model == 1 && num.solve_Navier_Stokes_liquid_phase == 1)
|
|
|
- nt = ni_uv + nb_uv
|
|
|
- ncol_A = ni_uv + nb_uv
|
|
|
- else
|
|
|
- nt = (num.nLS - num.nNavier + 1) * ni_uv + num.nNavier * ni_p + nb_uv
|
|
|
- ncol_A = (num.nLS - num.nNavier + 1) * ni_uv + num.nNavier * ni_p + nb_uv
|
|
|
- end
|
|
|
-
|
|
|
- AuvL = spzeros(ncol_A, nt)
|
|
|
- BuvL = spzeros(ncol_A, nt)
|
|
|
- rhs_uv = zeros(ncol_A)
|
|
|
- # AϕL = spzeros(nt, nt)
|
|
|
-
|
|
|
-
|
|
|
-
|
|
|
- elseif num.pressure_velocity_coupling == 1
|
|
|
- nt = (num.nLS - num.nNavier + 1) * ni_uv + num.nNavier * ni_p + nb_uv + (num.nLS + 1) * ni_p + nb_p
|
|
|
- ncol_A = (num.nLS - num.nNavier + 1) * ni_uv + num.nNavier * ni_p + nb_uv + ni_p
|
|
|
- # so 1 * ni + 1 * ni_p +nb + ni_p + nb
|
|
|
- # u v Navier, pression
|
|
|
-
|
|
|
- # AuvL = spzeros(nt, nt)
|
|
|
- # BuvL = spzeros(nt, nt)
|
|
|
-
|
|
|
- AuvL = spzeros(ncol_A, nt)
|
|
|
- BuvL = spzeros(ncol_A, nt)
|
|
|
- rhs_uv = zeros(ncol_A)
|
|
|
-
|
|
|
- elseif num.pressure_velocity_coupling == 2
|
|
|
- nt = (num.nLS - num.nNavier + 1) * ni_uv + num.nNavier * ni_p + nb_uv + (num.nLS + 1) * ni_p + nb_p
|
|
|
- ncol_A = nt
|
|
|
-
|
|
|
- # so 1 * ni + 1 * ni_p +nb + ni_p + nb
|
|
|
- # u v Navier, pression
|
|
|
-
|
|
|
- # AuvL = spzeros(nt, nt)
|
|
|
- # BuvL = spzeros(nt, nt)
|
|
|
-
|
|
|
- AuvL = spzeros(ncol_A, nt)
|
|
|
- BuvL = spzeros(ncol_A, nt)
|
|
|
- rhs_uv = zeros(ncol_A)
|
|
|
-
|
|
|
- elseif num.pressure_velocity_coupling == 3 #no BC for pressure
|
|
|
-
|
|
|
-
|
|
|
-
|
|
|
- if (num.one_fluid_model == 1 && num.solve_Navier_Stokes_liquid_phase == 1)
|
|
|
- nt = ni_uv + nb_uv + ni_p
|
|
|
- ncol_A = ni_uv + nb_uv + ni_p
|
|
|
- else
|
|
|
- nt = (num.nLS - num.nNavier + 1) * ni_uv + num.nNavier * ni_p + nb_uv + ni_p
|
|
|
- ncol_A = nt
|
|
|
- end
|
|
|
-
|
|
|
- AuvL = spzeros(ncol_A, nt)
|
|
|
- BuvL = spzeros(ncol_A, nt)
|
|
|
- rhs_uv = zeros(ncol_A)
|
|
|
- # AϕL = spzeros(nt, nt)
|
|
|
-
|
|
|
-
|
|
|
- elseif num.pressure_velocity_coupling == 4 # BC for pressure on interfaces (bubble)
|
|
|
-
|
|
|
- n_phase = 2
|
|
|
-
|
|
|
- nt = n_phase * ((num.nLS - num.nNavier + 1) * ni_uv + num.nNavier * ni_p + (num.nLS + 1) * ni_p ) + nb_uv
|
|
|
-
|
|
|
- ncol_A = nt
|
|
|
-
|
|
|
- # so 1 * ni + 1 * ni_p +nb + ni_p + nb
|
|
|
- # u v Navier, pression
|
|
|
-
|
|
|
- # AuvL = spzeros(nt, nt)
|
|
|
- # BuvL = spzeros(nt, nt)
|
|
|
-
|
|
|
- AuvL = spzeros(ncol_A, nt)
|
|
|
- BuvL = spzeros(ncol_A, nt)
|
|
|
- rhs_uv = zeros(ncol_A)
|
|
|
-
|
|
|
- # elseif num.pressure
|
|
|
-
|
|
|
- end
|
|
|
-
|
|
|
- end
|
|
|
-
|
|
|
- #endregion Allocate for Navier case
|
|
|
-
|
|
|
- #endregion Allocations
|
|
|
-
|
|
|
- #TODO remove allocations
|
|
|
-
|
|
|
-
|
|
|
- # Adapt cell volume W for gradients
|
|
|
- # cf 4/3 factor for Laplacian
|
|
|
- if num.laplacian == 1
|
|
|
- AvLcopy = copy(AvL)
|
|
|
- Lvm1_L,bc_Lvm1_L,bc_Lvm1_b_L=compute_divergence!(num,
|
|
|
- # grid_p,
|
|
|
- # grid_u,
|
|
|
- grid_v,
|
|
|
- op.opC_vL,
|
|
|
- AvLcopy,
|
|
|
- # rhs_scal,
|
|
|
- # tmp_vec_p, #a0
|
|
|
- tmp_vec_1D_v0,
|
|
|
- tmp_vec_1D_v,
|
|
|
- Lvm1_L,
|
|
|
- bc_Lvm1_L,
|
|
|
- bc_Lvm1_b_L
|
|
|
- # tmp_vec_u0,
|
|
|
- # tmp_vec_v0,
|
|
|
- # tmp_vec_1D,
|
|
|
- # ls_advection
|
|
|
- )
|
|
|
-
|
|
|
- ApLcopy = copy(Ascal)
|
|
|
- Lpm1_L,bc_Lpm1_L,bc_Lpm1_b_L=compute_divergence!(num,
|
|
|
- # grid_p,
|
|
|
- # grid_u,
|
|
|
- grid_p,
|
|
|
- op.opC_pL,
|
|
|
- ApLcopy,
|
|
|
- # rhs_scal,
|
|
|
- # tmp_vec_p, #a0
|
|
|
- tmp_vec_1D_p0,
|
|
|
- tmp_vec_1D_p,
|
|
|
- Lpm1_L,
|
|
|
- bc_Lpm1_L,
|
|
|
- bc_Lpm1_b_L
|
|
|
- # tmp_vec_u0,
|
|
|
- # tmp_vec_v0,
|
|
|
- # tmp_vec_1D,
|
|
|
- # ls_advection
|
|
|
- )
|
|
|
-
|
|
|
- end
|
|
|
-
|
|
|
- if num.pressure_velocity_coupling == 0
|
|
|
- #TODO why this call without interface initialization ?
|
|
|
- if !navier
|
|
|
- if (num.solve_solid == 1) && ns_solid_phase
|
|
|
- _ = FE_set_momentum(
|
|
|
- num, grid_u, op.opC_uS,
|
|
|
- AuS, BuS,
|
|
|
- iRe.*Lum1_S, iRe.*bc_Lum1_S, iRe.*bc_Lum1_b_S, Mum1_S, BC_uS,
|
|
|
- true
|
|
|
- )
|
|
|
-
|
|
|
- _ = FE_set_momentum(
|
|
|
- num, grid_v, op.opC_vS,
|
|
|
- AvS, BvS,
|
|
|
- iRe.*Lvm1_S, iRe.*bc_Lvm1_S, iRe.*bc_Lvm1_b_S, Mvm1_S, BC_vS,
|
|
|
- true
|
|
|
- )
|
|
|
- end
|
|
|
- else
|
|
|
- _ = FE_set_momentum_coupled(
|
|
|
- BC_int, num, grid_p, grid_u, grid_v,
|
|
|
- op.opC_pS, op.opC_uS, op.opC_vS,
|
|
|
- AuvS, BuvS,
|
|
|
- iRe.*Lum1_S, iRe.*bc_Lum1_S, iRe.*bc_Lum1_b_S, Mum1_S, BC_uS,
|
|
|
- iRe.*Lvm1_S, iRe.*bc_Lvm1_S, iRe.*bc_Lvm1_b_S, Mvm1_S, BC_vS,
|
|
|
- true
|
|
|
- )
|
|
|
- end
|
|
|
- elseif num.pressure_velocity_coupling > 1
|
|
|
- if num.pressure_velocity_coupling == 4
|
|
|
- # Coupled resolution of u and v, going to two-phases for pressure
|
|
|
- _ = FE_set_momentum_coupled_two_phases(
|
|
|
- BC_int, num, grid_p, grid_u, grid_v,
|
|
|
- op,
|
|
|
- AuvL, BuvL,rhs_uv,
|
|
|
- iRe.*Lum1_L, iRe.*bc_Lum1_L, iRe.*bc_Lum1_b_L, Mum1_L, BC_uL,
|
|
|
- iRe.*Lvm1_L, iRe.*bc_Lvm1_L, iRe.*bc_Lvm1_b_L, Mvm1_L, BC_vL,
|
|
|
- iRe.*Lum1_S, iRe.*bc_Lum1_S, iRe.*bc_Lum1_b_S, Mum1_S, BC_uS,
|
|
|
- iRe.*Lvm1_S, iRe.*bc_Lvm1_S, iRe.*bc_Lvm1_b_S, Mvm1_S, BC_vS,
|
|
|
- true,BC_pL,phL,phS
|
|
|
- )
|
|
|
- else
|
|
|
- if (num.one_fluid_model == 1 && num.solve_Navier_Stokes_liquid_phase == 1)
|
|
|
- # _ = FE_set_momentum_coupled2_one_fluid(
|
|
|
- # BC_int, num, grid_p, grid_u, grid_v,
|
|
|
- # op.opC_pL, op.opC_uL, op.opC_vL,
|
|
|
- # AuvL, BuvL,rhs_uv,
|
|
|
- # diffusion_bulk_u, diffusion_LS_u, diffusion_border_u, Mum1, BC_uL,
|
|
|
- # diffusion_bulk_v, diffusion_LS_v, diffusion_border_v, Mvm1, BC_vL,
|
|
|
- # cross_term_diffusion_bulk_d_dv_dx_dy,cross_term_diffusion_bulk_d_du_dy_dx,
|
|
|
- # cross_term_diffusion_bulk_d_dv_dx_dy_border,cross_term_diffusion_bulk_d_du_dy_dx_border,
|
|
|
- # rho_one_fluid_u,rho_one_fluid_v,
|
|
|
- # true,
|
|
|
- # BC_pL,phL
|
|
|
- # )
|
|
|
- else
|
|
|
-
|
|
|
- # Coupled resolution of u and v
|
|
|
- _ = FE_set_momentum_coupled2(
|
|
|
- BC_int, num, grid_p, grid_u, grid_v,
|
|
|
- op.opC_pL, op.opC_uL, op.opC_vL,
|
|
|
- AuvL, BuvL,rhs_uv,
|
|
|
- iRe.*Lum1_L, iRe.*bc_Lum1_L, iRe.*bc_Lum1_b_L, Mum1_L, BC_uL,
|
|
|
- iRe.*Lvm1_L, iRe.*bc_Lvm1_L, iRe.*bc_Lvm1_b_L, Mvm1_L, BC_vL,
|
|
|
- true,BC_pL,phL
|
|
|
- )
|
|
|
- end
|
|
|
- end
|
|
|
-
|
|
|
-
|
|
|
- end
|
|
|
-
|
|
|
-
|
|
|
- #TODO remove alloc a0_p
|
|
|
- a0_p = []
|
|
|
- for i in 1:num.nLS
|
|
|
- push!(a0_p, zeros(grid_p))
|
|
|
- end
|
|
|
-
|
|
|
-
|
|
|
- if !advection && (num.solve_solid == 1)
|
|
|
- #call to set_heat! is there to set up the matrices for the heat equation.
|
|
|
- #If the level-set is not advected, then after this call there is no need to update these matrices anymore
|
|
|
-
|
|
|
- _ = set_poisson(
|
|
|
- BC_int, num, grid_p, a0_p, op.opC_pS, op.opC_uS, op.opC_vS,
|
|
|
- AϕS, Lpm1_S, bc_Lpm1_S, bc_Lpm1_b_S, BC_pS,
|
|
|
- true
|
|
|
- )
|
|
|
-
|
|
|
- set_heat!(
|
|
|
- BC_int[1], num, grid_p, op.opC_TS, grid_p.LS[1].geoS, phS, num.θd, BC_TS, grid_p.LS[1].MIXED, grid_p.LS[1].geoS.projection,
|
|
|
- ATS, BTS,rhs_scal,
|
|
|
- op.opS, grid_u, grid_u.LS[1].geoS, grid_v, grid_v.LS[1].geoS,
|
|
|
- periodic_x, periodic_y, heat_convection, true, BC_int
|
|
|
- )
|
|
|
- end
|
|
|
-
|
|
|
- if test_laplacian
|
|
|
- num.exact_laplacian = test_laplacian_pressure(num,grid_v,phL.vD,op.opC_pL, Lvm1_L, bc_Lvm1_L, bc_Lvm1_b_L)
|
|
|
- return
|
|
|
- end
|
|
|
-
|
|
|
- # Segregated resolution for u and v since the viscosity is assumed to be constant in a phase
|
|
|
- if num.pressure_velocity_coupling == 0
|
|
|
- if !navier
|
|
|
- _ = FE_set_momentum(
|
|
|
- num, grid_u, op.opC_uL,
|
|
|
- AuL, BuL,
|
|
|
- iRe.*Lum1_L, iRe.*bc_Lum1_L, iRe.*bc_Lum1_b_L, Mum1_L, BC_uL,
|
|
|
- true
|
|
|
- )
|
|
|
- _ = FE_set_momentum(
|
|
|
- num, grid_v, op.opC_vL,
|
|
|
- AvL, BvL,
|
|
|
- iRe.*Lvm1_L, iRe.*bc_Lvm1_L, iRe.*bc_Lvm1_b_L, Mvm1_L, BC_vL,
|
|
|
- true
|
|
|
- )
|
|
|
- else
|
|
|
- # Coupled resolution of u and v if navier BC is activated
|
|
|
- _ = FE_set_momentum_coupled(
|
|
|
- BC_int, num, grid_p, grid_u, grid_v,
|
|
|
- op.opC_pL, op.opC_uL, op.opC_vL,
|
|
|
- AuvL, BuvL,
|
|
|
- iRe.*Lum1_L, iRe.*bc_Lum1_L, iRe.*bc_Lum1_b_L, Mum1_L, BC_uL,
|
|
|
- iRe.*Lvm1_L, iRe.*bc_Lvm1_L, iRe.*bc_Lvm1_b_L, Mvm1_L, BC_vL,
|
|
|
- true
|
|
|
- )
|
|
|
- end
|
|
|
-
|
|
|
- elseif num.pressure_velocity_coupling > 1
|
|
|
-
|
|
|
-
|
|
|
- if (num.one_fluid_model == 1 && num.solve_Navier_Stokes_liquid_phase == 1)
|
|
|
- # _ = FE_set_momentum_coupled2_one_fluid(
|
|
|
- # BC_int, num, grid_p, grid_u, grid_v,
|
|
|
- # op.opC_pL, op.opC_uL, op.opC_vL,
|
|
|
- # AuvL, BuvL,rhs_uv,
|
|
|
- # diffusion_bulk_u, diffusion_LS_u, diffusion_border_u, Mum1, BC_uL,
|
|
|
- # diffusion_bulk_v, diffusion_LS_v, diffusion_border_v, Mvm1, BC_vL,
|
|
|
- # cross_term_diffusion_bulk_d_dv_dx_dy,cross_term_diffusion_bulk_d_du_dy_dx,
|
|
|
- # cross_term_diffusion_bulk_d_dv_dx_dy_border,cross_term_diffusion_bulk_d_du_dy_dx_border,
|
|
|
- # rho_one_fluid_u,rho_one_fluid_v,
|
|
|
- # true,
|
|
|
- # BC_pL,phL
|
|
|
- # )
|
|
|
- else
|
|
|
- # Coupled resolution of u and v
|
|
|
- _ = FE_set_momentum_coupled2(
|
|
|
- BC_int, num, grid_p, grid_u, grid_v,
|
|
|
- op.opC_pL, op.opC_uL, op.opC_vL,
|
|
|
- AuvL, BuvL,rhs_uv,
|
|
|
- iRe.*Lum1_L, iRe.*bc_Lum1_L, iRe.*bc_Lum1_b_L, Mum1_L, BC_uL,
|
|
|
- iRe.*Lvm1_L, iRe.*bc_Lvm1_L, iRe.*bc_Lvm1_b_L, Mvm1_L, BC_vL,
|
|
|
- true,BC_pL,phL
|
|
|
- )
|
|
|
- end #one-fluid
|
|
|
-
|
|
|
- end
|
|
|
-
|
|
|
- a0_p = []
|
|
|
- for i in 1:num.nLS
|
|
|
- push!(a0_p, zeros(grid_p))
|
|
|
- end
|
|
|
-
|
|
|
- if num.one_fluid_model == 0
|
|
|
- _ = set_poisson(
|
|
|
- BC_int, num, grid_p, a0_p, op.opC_pL, op.opC_uL, op.opC_vL,
|
|
|
- AϕL, Lpm1_L, bc_Lpm1_L, bc_Lpm1_b_L, BC_pL,
|
|
|
- true
|
|
|
- )
|
|
|
- else
|
|
|
- _ = set_poisson_one_fluid(
|
|
|
- BC_int, num, grid_p, a0_p, op.opC_pL, op.opC_uL, op.opC_vL,
|
|
|
- AϕL, Lpm1_L, bc_Lpm1_L, bc_Lpm1_b_L, BC_pL,
|
|
|
- true
|
|
|
- )
|
|
|
- end
|
|
|
-
|
|
|
- if num.nLS>1 #monophasic
|
|
|
- #call to set_heat! is there to set up the matrices for the heat equation.
|
|
|
- #If the level-set is not advected, then after this call there is no need to update these matrices anymore
|
|
|
- set_heat!(
|
|
|
- BC_int[1], num, grid_p, op.opC_TL, grid_p.LS[1].geoL, phL, num.θd, BC_TL, grid_p.LS[1].MIXED, grid_p.LS[1].geoL.projection,
|
|
|
- ATL, BTL,rhs_scal,
|
|
|
- op.opL, grid_u, grid_u.LS[1].geoL, grid_v, grid_v.LS[1].geoL,
|
|
|
- periodic_x, periodic_y, heat_convection, true, BC_int
|
|
|
- )
|
|
|
- end
|
|
|
-
|
|
|
- end
|
|
|
- else
|
|
|
- error("Unknown time scheme. Available options are ForwardEuler and CrankNicolson")
|
|
|
- end
|
|
|
-
|
|
|
- if heat_convection || electrolysis_convection
|
|
|
- NB_indices = update_all_ls_data(num, grid_p, grid_u, grid_v, BC_int, periodic_x, periodic_y)
|
|
|
- # println(grid_p.LS[1].geoL.dcap[1,1,:])
|
|
|
- end
|
|
|
-
|
|
|
- # V0S = volume(grid_p.LS[end].geoS)
|
|
|
- # V0L = volume(grid_p.LS[end].geoL)
|
|
|
-
|
|
|
- if num.debug == "allocations_start"
|
|
|
- print("\n STOP allocations")
|
|
|
- return
|
|
|
- end
|
|
|
-
|
|
|
- #region restart
|
|
|
- # - file: decl_hdf5_test_02_r${rank}.h5
|
|
|
- # when: $input=1
|
|
|
- # read: [ reals, values ]
|
|
|
-
|
|
|
- PDI_status = @ccall "libpdi".PDI_multi_expose("restart"::Cstring,
|
|
|
- "nstep"::Cstring, num.current_iter::Ref{Clonglong}, PDI_OUT::Cint,
|
|
|
- "time"::Cstring, num.time::Ref{Cdouble}, PDI_OUT::Cint,
|
|
|
- "timestep"::Cstring, num.timestep_n::Ref{Cdouble}, PDI_OUT::Cint,
|
|
|
- "nx"::Cstring, grid_p.nx::Ref{Clonglong}, PDI_OUT::Cint,
|
|
|
- "ny"::Cstring, grid_p.ny::Ref{Clonglong}, PDI_OUT::Cint,
|
|
|
- "u_1D"::Cstring, phL.uD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "v_1D"::Cstring, phL.vD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "p_1D"::Cstring, phL.pD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "levelset_p"::Cstring, grid_p.LS[num.iLSpdi].u::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "levelset_u"::Cstring, grid_u.LS[num.iLSpdi].u::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "levelset_v"::Cstring, grid_v.LS[num.iLSpdi].u::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "trans_scal_1DT"::Cstring, phL.trans_scalD'::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "phi_ele_1D"::Cstring, phL.phi_eleD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "velocity_x"::Cstring, tmp_vec_p::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "velocity_y"::Cstring, tmp_vec_p0::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "radius"::Cstring, num.current_radius::Ref{Cdouble}, PDI_OUT::Cint,
|
|
|
- C_NULL::Ptr{Cvoid})::Cint
|
|
|
-
|
|
|
- #endregion restart
|
|
|
-
|
|
|
- interpolate_staggered_u_v_to_scalar_grid_one_fluid_or_one_phase!(num,grid_p,grid_u,grid_v,phL.u,phL.v,tmp_vec_p,tmp_vec_p0)
|
|
|
-
|
|
|
- PDI_status = @ccall "libpdi".PDI_multi_expose("write_data"::Cstring,
|
|
|
- "nstep"::Cstring, num.current_iter::Ref{Clonglong}, PDI_OUT::Cint,
|
|
|
- "time"::Cstring, num.time::Ref{Cdouble}, PDI_OUT::Cint,
|
|
|
- "timestep"::Cstring, num.timestep_n::Ref{Cdouble}, PDI_OUT::Cint,
|
|
|
- "nx"::Cstring, grid_p.nx::Ref{Clonglong}, PDI_OUT::Cint,
|
|
|
- "ny"::Cstring, grid_p.ny::Ref{Clonglong}, PDI_OUT::Cint,
|
|
|
- "u_1D"::Cstring, phL.uD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "v_1D"::Cstring, phL.vD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "p_1D"::Cstring, phL.pD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "levelset_p"::Cstring, grid_p.LS[num.iLSpdi].u::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "levelset_u"::Cstring, grid_u.LS[num.iLSpdi].u::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "levelset_v"::Cstring, grid_v.LS[num.iLSpdi].u::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "trans_scal_1DT"::Cstring, phL.trans_scalD'::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "phi_ele_1D"::Cstring, phL.phi_eleD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "velocity_x"::Cstring, tmp_vec_p::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "velocity_y"::Cstring, tmp_vec_p0::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "radius"::Cstring, num.current_radius::Ref{Cdouble}, PDI_OUT::Cint,
|
|
|
- C_NULL::Ptr{Cvoid})::Cint
|
|
|
-
|
|
|
- #compute numerical radius
|
|
|
- PDI_status = @ccall "libpdi".PDI_multi_expose("compute_radius"::Cstring,
|
|
|
- "levelset_p"::Cstring, grid_p.LS[num.iLSpdi].u::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "mesh_p_x"::Cstring, grid_p.x::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "mesh_p_y"::Cstring, grid_p.y::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "radius_vec"::Cstring, radius_pdi::Ptr{Cdouble}, PDI_INOUT::Cint,
|
|
|
- C_NULL::Ptr{Cvoid})::Cint
|
|
|
- num.current_radius = radius_pdi[1]
|
|
|
-
|
|
|
-
|
|
|
- print("\n radius pdi ", radius_pdi[1])
|
|
|
-
|
|
|
- # slice = levelset_p[nx//2,:]
|
|
|
- # x_1D = mesh_p_y[nx//2,:]
|
|
|
-
|
|
|
- # try
|
|
|
- # radius_vertical = compute_radius_from_levelset_slice(grid_p.LS[num.iLSpdi].u[:,div(grid_p.nx,2)],
|
|
|
- # grid_p.y[:,div(grid_p.nx,2)])
|
|
|
- # # catch error
|
|
|
- # # volume_cell = grid_p.LS[num.iLSpdi].geoS.cap[:,:,5] #bubble
|
|
|
- # volume_cell = grid_p.LS[num.iLSpdi].geoL.cap[:,:,5] #drop
|
|
|
-
|
|
|
- # center_of_mass_x, center_of_mass_y = calculate_centroid(grid_p.x, grid_p.y, volume_cell)
|
|
|
-
|
|
|
-
|
|
|
- # indices_bubble_mass_center = find_slice_coord_bubble_mass_center(center_of_mass_x,center_of_mass_y,
|
|
|
- # num,grid_p)
|
|
|
- # radius_horizontal = compute_radius_from_levelset_slice(grid_p.LS[num.iLSpdi].u[indices_bubble_mass_center[1],:],
|
|
|
- # grid_p.y[indices_bubble_mass_center[1],:])
|
|
|
-
|
|
|
- # if isnothing(radius_vertical)
|
|
|
- # print("\n error radius vertical")
|
|
|
- # if isnothing(radius_horizontal)
|
|
|
- # print("\n error radius horizontal and vertical")
|
|
|
- # else
|
|
|
- # num.current_radius = radius_horizontal
|
|
|
- # end
|
|
|
-
|
|
|
- # else
|
|
|
- # if isnothing(radius_horizontal)
|
|
|
- # print("\n error radius horizontal")
|
|
|
- # else
|
|
|
- # num.current_radius = max(radius_horizontal, radius_vertical)
|
|
|
- # end
|
|
|
- # end
|
|
|
- # # end
|
|
|
-
|
|
|
- # num.current_radius = maximum(compute_radii_from_slices(
|
|
|
- # grid_p.LS[num.iLSpdi].u,
|
|
|
- # x_grid::AbstractMatrix,
|
|
|
- # y_grid::AbstractMatrix,
|
|
|
- # slices::Vector{Tuple{Union{Int, Colon}, Union{Int, Colon}}};
|
|
|
- # center_x,center_y)
|
|
|
-
|
|
|
- if num.sphere_post_processing == 1
|
|
|
- compute_bubble_drop_radius(num, grid_p)
|
|
|
- end
|
|
|
-
|
|
|
- if (num.one_fluid_model == 1 && num.solve_Navier_Stokes_liquid_phase == 1 && num.solve_Navier_Stokes == 1)
|
|
|
- # rise_velocity_y =0.0
|
|
|
- # PDI_status = @ccall "libpdi".PDI_multi_expose("post_processing_rising_bubble_first_share"::Cstring,
|
|
|
- # "nstep"::Cstring, num.current_iter ::Ref{Clonglong}, PDI_OUT::Cint,
|
|
|
- # # "rho_one_fluid"::Cstring, rho_one_fluid::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # # "rho_one_fluid_u"::Cstring, rho_one_fluid_u::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # # "rho_one_fluid_v"::Cstring, rho_one_fluid_v::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # # "mu_one_fluid"::Cstring, mu_one_fluid::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # "rise_velocity_y"::Cstring, rise_velocity_y::Ref{Cdouble}, PDI_OUT::Cint,
|
|
|
- # "timestep"::Cstring, num.timestep_n::Ref{Cdouble}, PDI_OUT::Cint,
|
|
|
- # # "volume_fraction"::Cstring, volume_fraction::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # # "volume_cell"::Cstring, grid_p.LS[end].geoS.dcap[:,:,5]::Ptr{Cdouble}, PDI_OUT::Cint, #geoS for bubble phase
|
|
|
- # # "mesh_p_x"::Cstring, grid_p.x::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # # "mesh_p_y"::Cstring, grid_p.y::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # # "dcap_1"::Cstring, grid_p.LS[num.iLSpdi].geoS.dcap[:,:,1]::Ptr{Cdouble}, PDI_OUT::Cint, #geoS for bubble phase
|
|
|
- # # "dcap_2"::Cstring, grid_p.LS[num.iLSpdi].geoS.dcap[:,:,2]::Ptr{Cdouble}, PDI_OUT::Cint, #geoS for bubble phase
|
|
|
- # # "dcap_3"::Cstring, grid_p.LS[num.iLSpdi].geoS.dcap[:,:,3]::Ptr{Cdouble}, PDI_OUT::Cint, #geoS for bubble phase
|
|
|
- # # "dcap_4"::Cstring, grid_p.LS[num.iLSpdi].geoS.dcap[:,:,4]::Ptr{Cdouble}, PDI_OUT::Cint, #geoS for bubble phase
|
|
|
- # C_NULL::Ptr{Cvoid})::Cint
|
|
|
-
|
|
|
- update_one_fluid_density_viscosity(num,grid_p,grid_u,grid_v,volume_fraction,levelset_one_fluid,rho_one_fluid,
|
|
|
- rho_one_fluid_u,rho_one_fluid_v,mu_one_fluid,tmp_vec_p0)
|
|
|
-
|
|
|
- # PDI_status = @ccall "libpdi".PDI_multi_expose("print_timestep"::Cstring,
|
|
|
- # "nstep"::Cstring, num.current_iter ::Ref{Clonglong}, PDI_OUT::Cint,
|
|
|
- # "time"::Cstring, num.time::Ref{Cdouble}, PDI_OUT::Cint,
|
|
|
- # "timestep"::Cstring, num.timestep_n::Ref{Cdouble}, PDI_OUT::Cint,
|
|
|
- # C_NULL::Ptr{Cvoid})::Cint
|
|
|
-
|
|
|
- # #region initial values
|
|
|
- # PDI_status = @ccall "libpdi".PDI_multi_expose("post_processing_rising_bubble"::Cstring,
|
|
|
- # #endregion initial values
|
|
|
-
|
|
|
- end
|
|
|
-
|
|
|
-
|
|
|
- if num.one_fluid_model == 1 && num.solve_Navier_Stokes == 1
|
|
|
- conservation = compute_conservation_mass(num,phL, grid_p ,grid_u, grid_v, rho_one_fluid)
|
|
|
- else
|
|
|
- conservation = 0 # TODO
|
|
|
- end
|
|
|
-
|
|
|
- # Compute divergence of velocity
|
|
|
- Duv = op.opC_pL.AxT * vec1(phL.uD,grid_u) .+ op.opC_pL.Gx_b * vecb(phL.uD,grid_u) .+
|
|
|
- op.opC_pL.AyT * vec1(phL.vD,grid_v) .+ op.opC_pL.Gy_b * vecb(phL.vD,grid_v)
|
|
|
- for iLS in 1:num.nLS
|
|
|
- if !is_navier(BC_int[iLS]) && !is_navier_cl(BC_int[iLS]) #otherwise normal velocity null if no blowing
|
|
|
- Duv .+= op.opC_pL.Gx[iLS] * veci(phL.uD,grid_u,iLS+1) .+
|
|
|
- op.opC_pL.Gy[iLS] * veci(phL.vD,grid_v,iLS+1)
|
|
|
- end
|
|
|
- end
|
|
|
-
|
|
|
-
|
|
|
- if num.solve_Navier_Stokes == 1 && num.solve_Navier_Stokes == 1
|
|
|
- PDI_status = @ccall "libpdi".PDI_multi_expose("print_conservation"::Cstring,
|
|
|
- "nstep"::Cstring, num.current_iter::Ref{Clonglong}, PDI_OUT::Cint,
|
|
|
- "conservation"::Cstring, conservation::Ref{Cdouble}, PDI_OUT::Cint,
|
|
|
- "velocity_divergence"::Cstring, Duv::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # "p_1D"::Cstring, phL.pD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- C_NULL::Ptr{Cvoid})::Cint
|
|
|
-
|
|
|
- #save initialised electrical potential and current (iter 0)
|
|
|
- if num.electrical_potential>0
|
|
|
- compute_grad_phi_ele!(num, grid_p, grid_u, grid_v, grid_u.LS[end], grid_v.LS[end], phL,
|
|
|
- op.opC_pL, elec_cond,tmp_vec_u,tmp_vec_v,tmp_vec_p,tmp_vec_p0,tmp_vec_p1) #TODO current
|
|
|
- end
|
|
|
-
|
|
|
- # PDI_status = @ccall "libpdi".PDI_multi_expose("check_conservation"::Cstring,
|
|
|
- # "nstep"::Cstring, num.current_iter::Ref{Clonglong}, PDI_OUT::Cint,
|
|
|
- # "conservation"::Cstring, conservation::Ref{Cdouble}, PDI_OUT::Cint,
|
|
|
- # "velocity_divergence"::Cstring, Duv::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # "u_1D"::Cstring, phL.uD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # "v_1D"::Cstring, phL.vD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # "p_1D"::Cstring, phL.pD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # C_NULL::Ptr{Cvoid})::Cint
|
|
|
-
|
|
|
- end
|
|
|
-
|
|
|
-
|
|
|
- simulation_finished = false
|
|
|
-
|
|
|
- #TODO variable time steps
|
|
|
-
|
|
|
- #region time loop
|
|
|
- while (num.current_iter < num.max_iterations + 1) && (num.time < num.end_time) && (num.stop_simulation == 0)
|
|
|
-
|
|
|
- #update time
|
|
|
- num.time += num.timestep_n
|
|
|
- num.current_iter += 1
|
|
|
-
|
|
|
-
|
|
|
- #region start iter
|
|
|
-
|
|
|
- #region Adapt timestep
|
|
|
-
|
|
|
-
|
|
|
- num.timestep_n = adapt_timestep!(num, phL, phS, grid_u, grid_v,adapt_timestep_mode)
|
|
|
-
|
|
|
- # printstyled(color=:green, @sprintf "\n num.CFL : %.2e dt : %.2e num.timestep_n : %.2e\n" num.CFL num.timestep_n num.timestep_n)
|
|
|
- # print("\n num.stop_simulation start loop ",num.stop_simulation)
|
|
|
-
|
|
|
- #endregion adapt time
|
|
|
-
|
|
|
- if grid_p.LS[1].geoL.dcap[1,1,:] == 0.0
|
|
|
- printstyled(color=:red, @sprintf "\n Error operator null \n")
|
|
|
- return
|
|
|
- end
|
|
|
-
|
|
|
- # Print information at the start of temporal iteration and check definition of operators
|
|
|
- # cf update_all_ls_data (true VS false)
|
|
|
- PDI_status = @ccall "libpdi".PDI_multi_expose("print_start_temporal_iteration"::Cstring,
|
|
|
- "nstep"::Cstring, num.current_iter::Ref{Clonglong}, PDI_OUT::Cint,
|
|
|
- "time"::Cstring, num.time::Ref{Cdouble}, PDI_OUT::Cint,
|
|
|
- # "levelset_p"::Cstring, grid_p.LS[num.index_levelset_pdi].u::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "dcap"::Cstring, grid_p.LS[num.index_levelset_pdi].geoL.dcap[:,:,:]::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # grid_p.LS[1].geoL.dcap[1,1,:]
|
|
|
- C_NULL::Ptr{Cvoid})::Cint
|
|
|
-
|
|
|
-
|
|
|
-
|
|
|
-
|
|
|
- # if num.io_pdi>0
|
|
|
-
|
|
|
- # #or permutedims(grid_p.LS[num.iLSpdi].geoL.dcap, (3, 2, 1)) (3, 1, 2)
|
|
|
- # try
|
|
|
- # # num.iLSpdi = 1 # TODO all grid_p.LS
|
|
|
- # PDI_status = @ccall "libpdi".PDI_multi_expose("write_capacities"::Cstring,
|
|
|
- # # "dcap"::Cstring, permutedims(grid_p.LS[num.iLSpdi].geoL.dcap, (3, 2, 1))::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # "dcap_1"::Cstring, grid_p.LS[num.iLSpdi].geoL.dcap[:,:,1]::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # "dcap_2"::Cstring, grid_p.LS[num.iLSpdi].geoL.dcap[:,:,2]::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # "dcap_3"::Cstring, grid_p.LS[num.iLSpdi].geoL.dcap[:,:,3]::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # "dcap_4"::Cstring, grid_p.LS[num.iLSpdi].geoL.dcap[:,:,4]::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
-
|
|
|
- # C_NULL::Ptr{Cvoid})::Cint
|
|
|
- # catch error
|
|
|
- # printstyled(color=:red, @sprintf "\n PDI error \n")
|
|
|
- # print(error)
|
|
|
- # printstyled(color=:red, @sprintf "\n PDI error \n")
|
|
|
- # end
|
|
|
- # end #if io_pdi
|
|
|
-
|
|
|
- #endregion start iter
|
|
|
-
|
|
|
- # PDI (IO)
|
|
|
- if electrolysis
|
|
|
-
|
|
|
- #TODO not necessary to expose everything now for ex only grid_p.LS ? and the rest later
|
|
|
-
|
|
|
- intfc_vtx_x,intfc_vtx_y,intfc_vtx_field,intfc_vtx_connectivities,intfc_vtx_num, intfc_seg_num = convert_interfacial_D_to_segments(num,grid_p,phL.TD,1,2)
|
|
|
-
|
|
|
- #TODO when to write elec dat, ...
|
|
|
-
|
|
|
- if num.io_pdi>0
|
|
|
-
|
|
|
- try
|
|
|
- # printstyled(color=:red, @sprintf "\n PDI test \n" )
|
|
|
-
|
|
|
-
|
|
|
-
|
|
|
- # phi_array=phL.phi_ele #do not transpose since python row major
|
|
|
-
|
|
|
-
|
|
|
- if num.electrical_potential > 0
|
|
|
-
|
|
|
- # print("\n type of elec_cond ", typeof(elec_cond)," \n")
|
|
|
- try
|
|
|
- compute_grad_phi_ele!(num, grid_p, grid_u, grid_v,grid_u.LS[end], grid_v.LS[end],
|
|
|
- phL,
|
|
|
- # phS,
|
|
|
- op.opC_pL,
|
|
|
- # op.opC_pS,
|
|
|
- elec_cond,tmp_vec_u,tmp_vec_v,tmp_vec_p,tmp_vec_p0,tmp_vec_p1) #TODO current
|
|
|
- catch
|
|
|
- @error("compute_grad_phi_ele!")
|
|
|
- break
|
|
|
- end
|
|
|
- end
|
|
|
- # #store in us, vs instead of Eus, Evs
|
|
|
- # interpolate_grid_liquid!(grid_p,grid_u,grid_v,phL.Eu, phL.Ev,tmp_vec_p,tmp_vec_p0)
|
|
|
-
|
|
|
- # @ccall "libpdi".PDI_multi_expose("write_data_elec"::Cstring,
|
|
|
- # "i_current_x"::Cstring, tmp_vec_p::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # "i_current_y"::Cstring, tmp_vec_p0::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # "i_current_mag"::Cstring, phL.i_current_mag::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # "phi_ele_1D"::Cstring, phL.phi_eleD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # C_NULL::Ptr{Cvoid})::Cint
|
|
|
-
|
|
|
- interpolate_staggered_u_v_to_scalar_grid_one_fluid_or_one_phase!(num,grid_p,grid_u,grid_v,phL.u,phL.v,tmp_vec_p,tmp_vec_p0)
|
|
|
-
|
|
|
- # print("\n before write \n ")
|
|
|
-
|
|
|
- # num.iLSpdi = 1 # all grid_p.LS iLS = 1 # or all grid_p.LS ?
|
|
|
-
|
|
|
- # Exposing data to PDI for IO
|
|
|
- # if writing "D" array (bulk, interface, border), add "_1D" to the name
|
|
|
-
|
|
|
- # printstyled(color=:magenta, @sprintf "\n PDI write_data_start_loop %.5i \n" num.current_iter)
|
|
|
-
|
|
|
- PDI_status = @ccall "libpdi".PDI_multi_expose("write_data_start_loop"::Cstring,
|
|
|
- "nstep"::Cstring, num.current_iter::Ref{Clonglong}, PDI_OUT::Cint,
|
|
|
- "time"::Cstring, num.time::Ref{Cdouble}, PDI_OUT::Cint,
|
|
|
- "u_1D"::Cstring, phL.uD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "v_1D"::Cstring, phL.vD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "p_1D"::Cstring, phL.pD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "levelset_p"::Cstring, grid_p.LS[num.iLSpdi].u::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "levelset_u"::Cstring, grid_u.LS[num.iLSpdi].u::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "levelset_v"::Cstring, grid_v.LS[num.iLSpdi].u::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "trans_scal_1DT"::Cstring, phL.trans_scalD'::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "velocity_x"::Cstring, tmp_vec_p::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "velocity_y"::Cstring, tmp_vec_p0::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "radius"::Cstring, num.current_radius::Ref{Cdouble}, PDI_OUT::Cint,
|
|
|
- "intfc_vtx_num"::Cstring, intfc_vtx_num::Ref{Clonglong}, PDI_OUT::Cint,
|
|
|
- "intfc_seg_num"::Cstring, intfc_seg_num::Ref{Clonglong}, PDI_OUT::Cint,
|
|
|
- "intfc_vtx_x"::Cstring, intfc_vtx_x::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "intfc_vtx_y"::Cstring, intfc_vtx_y::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "intfc_vtx_field"::Cstring, intfc_vtx_field::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "intfc_vtx_connectivities"::Cstring, intfc_vtx_connectivities::Ptr{Clonglong}, PDI_OUT::Cint,
|
|
|
- C_NULL::Ptr{Cvoid})::Cint
|
|
|
-
|
|
|
- catch error
|
|
|
- printstyled(color=:red, @sprintf "\n PDI error \n")
|
|
|
- print(error)
|
|
|
- printstyled(color=:red, @sprintf "\n PDI error \n")
|
|
|
- end
|
|
|
-
|
|
|
- end #if io_pdi
|
|
|
- end #if electrolysis
|
|
|
-
|
|
|
-
|
|
|
- if !stefan
|
|
|
- grid_p.V .= speed #*ones(grid_p.ny, grid_p.nx)
|
|
|
- end
|
|
|
-
|
|
|
- #region Heat equation
|
|
|
- # Solve heat equation
|
|
|
- if heat
|
|
|
- if heat_solid_phase
|
|
|
- kill_dead_cells!(phS.T, grid_p, grid_p.LS[1].geoS)
|
|
|
- veci(phS.TD,grid_p,1) .= vec(phS.T)
|
|
|
- set_heat!(
|
|
|
- BC_int[1], num, grid_p, op.opC_TS, grid_p.LS[1].geoS, phS, num.θd, BC_TS, grid_p.LS[1].MIXED, grid_p.LS[1].geoS.projection,
|
|
|
- ATS, BTS,rhs_scal,
|
|
|
- op.opS, grid_u, grid_u.LS[1].geoS, grid_v, grid_v.LS[1].geoS,
|
|
|
- periodic_x, periodic_y, heat_convection, advection, BC_int
|
|
|
- )
|
|
|
- mul!(rhs_scal, BTS, phS.TD, 1.0, 1.0)
|
|
|
-
|
|
|
- phS.TD .= ATS \ rhs_scal
|
|
|
- phS.T .= reshape(veci(phS.TD,grid_p,1), grid_p)
|
|
|
- end
|
|
|
- if heat_liquid_phase
|
|
|
- kill_dead_cells!(phL.T, grid_p, grid_p.LS[1].geoL)
|
|
|
- veci(phL.TD,grid_p,1) .= vec(phL.T)
|
|
|
- set_heat!(
|
|
|
- BC_int[1], num, grid_p, op.opC_TL, grid_p.LS[1].geoL, phL, num.θd, BC_TL, grid_p.LS[1].MIXED, grid_p.LS[1].geoL.projection,
|
|
|
- ATL, BTL,rhs_scal,
|
|
|
- op.opL, grid_u, grid_u.LS[1].geoL, grid_v, grid_v.LS[1].geoL,
|
|
|
- periodic_x, periodic_y, heat_convection, advection, BC_int
|
|
|
- )
|
|
|
- mul!(rhs_scal, BTL, phL.TD, 1.0, 1.0)
|
|
|
-
|
|
|
- phL.TD .= ATL \ rhs_scal
|
|
|
- phL.T .= reshape(veci(phL.TD,grid_p,1), grid_p)
|
|
|
-
|
|
|
- end
|
|
|
- end # #if heat
|
|
|
- #endregion heat equation
|
|
|
-
|
|
|
-
|
|
|
- #region Electrolysis
|
|
|
- #TODO scalar without electrical potential
|
|
|
- if electrolysis && (num.solve_potential == 1 ) && (num.solve_species == 1)
|
|
|
- if electrolysis_liquid_phase
|
|
|
-
|
|
|
- #region Electrolysis: Poisson
|
|
|
- if num.electrical_potential > 0
|
|
|
- solve_poisson_loop!(num, grid_p, grid_u, grid_v, op, Ascal, rhs_scal,F_residual,
|
|
|
- tmp_vec_p,tmp_vec_p0,tmp_vec_p1, a1_p, BC_phi_ele, phL, phS,elec_cond,
|
|
|
- elec_condD, tmp_vec_u, tmp_vec_v, i_butler, ls_advection, heat)
|
|
|
- end
|
|
|
- #endregion Poisson
|
|
|
-
|
|
|
-
|
|
|
-
|
|
|
- #TODO check BC not overwritten by different scalars
|
|
|
- #TODO check ls advection true when num.n scalars
|
|
|
-
|
|
|
- # print("\n vecb_L",vecb_L(phL.phi_eleD, grid_p))
|
|
|
- # print("\n phL.phi_ele[:,1]",phL.phi_ele[:,1])
|
|
|
-
|
|
|
- #TODO phL.phi_ele[:,1] or vecb_L(phL.phi_eleD, grid_p)
|
|
|
- # we suppose phi(x=0)=... cf Khalighi
|
|
|
- #but here BC
|
|
|
-
|
|
|
- # New start scalar loop
|
|
|
-
|
|
|
- #region velocity
|
|
|
- # TODO
|
|
|
- interpolate_interface_velocity!(phL,grid_u,grid_v)
|
|
|
-
|
|
|
- #endregion velocity
|
|
|
-
|
|
|
- #region Impose velocity
|
|
|
- if imposed_velocity == "zero"
|
|
|
- phL.u .= 0.0
|
|
|
- phL.v .= 0.0
|
|
|
- phL.uD .= 0.0
|
|
|
- phL.vD .= 0.0
|
|
|
- phL.p .= 0.0
|
|
|
- phL.pD .= 0.0
|
|
|
-
|
|
|
-
|
|
|
- elseif imposed_velocity == "constant"
|
|
|
-
|
|
|
- #Required to modify whole uD vD
|
|
|
- phL.u .= 0.0
|
|
|
- phL.v .= BC_vL.bottom.val
|
|
|
- phL.uD .= 0.0
|
|
|
- phL.vD .= BC_vL.bottom.val
|
|
|
-
|
|
|
- phL.pD .= 0.0
|
|
|
-
|
|
|
- if ((num.current_iter-1)%show_every == 0)
|
|
|
- printstyled(color=:red, @sprintf "\n Imposed velocity v min %.2e max %.2e\n" minimum(phL.vD) maximum(phL.vD))
|
|
|
- printstyled(color=:red, @sprintf "\n Imposed velocity u min %.2e max %.2e\n" minimum(phL.uD) maximum(phL.uD))
|
|
|
- end
|
|
|
-
|
|
|
- elseif imposed_velocity == "Poiseuille_bottom_top"
|
|
|
-
|
|
|
- vPoiseuille = Poiseuille_fmax.(grid_v.x,num.v_inlet,num.L0)
|
|
|
-
|
|
|
-
|
|
|
- #Required to modify whole uD vD
|
|
|
- phL.u .= 0.0
|
|
|
- phL.v .= vPoiseuille
|
|
|
-
|
|
|
- phL.uD .= 0.0 #Dirichlet and Neumann
|
|
|
- # phL.vD .= BC_vL.bottom.val
|
|
|
- # vecb...TODO
|
|
|
- # ...
|
|
|
-
|
|
|
- elseif imposed_velocity == "radial"
|
|
|
- impose_radial_velocity(phS,phL,num)
|
|
|
- end #imposed_velocity
|
|
|
- #endregion Impose velocity
|
|
|
-
|
|
|
- #region Update current
|
|
|
-
|
|
|
- # if num.electrolysis_reaction == "Butler_no_concentration"
|
|
|
-
|
|
|
- # update_electrical_current_from_Butler_Volmer!(num,grid_p,heat,phL.phi_eleD,i_butler;phL.T)
|
|
|
-
|
|
|
- # # #TODO dev multiple levelsets
|
|
|
- # # if heat
|
|
|
- # # i_butler = butler_volmer_no_concentration.(num.alpha_a,num.alpha_c,num.Faraday,num.i0,vecb_L(phL.phi_eleD, grid_p),
|
|
|
- # # num.phi_ele1,num.Ru,phL.T)
|
|
|
- # # else
|
|
|
- # # if num.nLS == 1
|
|
|
- # # i_butler = butler_volmer_no_concentration.(num.alpha_a,num.alpha_c,num.Faraday,num.i0,vecb_L(phL.phi_eleD, grid_p),
|
|
|
- # # num.phi_ele1,num.Ru,num.temperature0)
|
|
|
- # # # else
|
|
|
- # # #imposed by LS 2
|
|
|
- # # # iLS_elec = 2
|
|
|
- # # # i_butler = butler_volmer_no_concentration.(num.alpha_a,num.alpha_c,num.Faraday,num.i0,veci(phL.phi_eleD, grid_p,iLS_elec+1),
|
|
|
- # # # num.phi_ele1,num.Ru,num.temperature0)
|
|
|
- # # end
|
|
|
- # # end
|
|
|
- # end
|
|
|
-
|
|
|
- # if num.electrolysis_reaction_symb in (:Butler_no_concentration, :fixed_current)
|
|
|
- if num.electrolysis_reaction_symb === :Butler_no_concentration
|
|
|
- # if num.verbosity>0
|
|
|
- # print("\n electrode_definition_function ",electrode_definition_function)
|
|
|
- # end
|
|
|
-
|
|
|
- update_electrical_current_from_Butler_Volmer_func!(num,grid_p,heat,phL.phi_eleD,i_butler,electrode_definition_function;phL.T)
|
|
|
- #if fixed do not update
|
|
|
- end
|
|
|
-
|
|
|
- #endregion Update current
|
|
|
-
|
|
|
-
|
|
|
- #region Scalar transport: update boundary conditions
|
|
|
- if num.nb_transported_scalars>0
|
|
|
-
|
|
|
- # printstyled(color=:magenta, @sprintf "\n num.nb_transported_scalars %.5i " num.nb_transported_scalars)
|
|
|
-
|
|
|
- for iscal=1:num.nb_transported_scalars
|
|
|
- # @views kill_dead_cells_val!(phL.trans_scal[:,:,iscal], grid_p, grid_p.LS[1].geoL,0.0)
|
|
|
- # @views kill_dead_cells_val!(phL.trans_scal[:,:,iscal], grid_p, grid_p.LS[1].geoL,num.concentration0[iscal])
|
|
|
-
|
|
|
- if num.kill_dead_cells == 0
|
|
|
- @views kill_dead_cells_val!(phL.trans_scal[:,:,iscal], grid_p, grid_p.LS[end].geoL,0.0)
|
|
|
- else
|
|
|
- @views kill_dead_cells_val!(phL.trans_scal[:,:,iscal], grid_p, grid_p.LS[end].geoL,num.concentration0[iscal])
|
|
|
- end
|
|
|
-
|
|
|
- @views veci(phL.trans_scalD[:,iscal],grid_p,1) .= vec(phL.trans_scal[:,:,iscal])
|
|
|
-
|
|
|
-
|
|
|
- if num.nLS == 1
|
|
|
- #BC for LS 2 in scalar transport : done in scalar loop
|
|
|
- # if iscal==1 || iscal==2
|
|
|
- # inv_stoechiometric_coeff = -1.0/2.0 #H2 and KOH
|
|
|
- # elseif iscal == 3
|
|
|
- # inv_stoechiometric_coeff = 1.0 #H2O consummed
|
|
|
- # end
|
|
|
- inv_stoechiometric_coeff = num.inv_stoechiometric_coeff[iscal]
|
|
|
-
|
|
|
- if num.electrolysis_reaction_symb === :Butler_no_concentration
|
|
|
-
|
|
|
- # BC at left wall
|
|
|
- # -(-/i_butler) because i=-lambda grad phi and BC at left: -e_x
|
|
|
- BC_trans_scal[iscal].left.val = i_butler./(num.Faraday*num.diffusion_coeff[iscal])*inv_stoechiometric_coeff
|
|
|
-
|
|
|
- # debug
|
|
|
- # for testn in 1:grid_p.ny
|
|
|
- # printstyled(color=:green, @sprintf "\n jtmp : %.5i j : %.5i border %.5e\n" testn grid_p.ny-testn+1 vecb_L(phL.trans_scalD[:,iscal], grid_p)[testn])
|
|
|
- # end
|
|
|
-
|
|
|
- # elseif num.electrolysis_reaction_symb === :fixed_current && num.nLS == 1
|
|
|
- #already initialised in convergence.jl
|
|
|
- # print("\n scalar ", iscal, " ", BC_trans_scal[iscal].bottom.val)
|
|
|
- # BC_trans_scal[iscal].bottom.val = i_butler./(num.Faraday*num.diffusion_coeff[iscal])*inv_stoechiometric_coeff
|
|
|
-
|
|
|
- end
|
|
|
-
|
|
|
- end
|
|
|
-
|
|
|
- # if num.scalar_transport_implementation > 0
|
|
|
-
|
|
|
- # printstyled(color=:red, @sprintf "\n scalar_transport_implementation")
|
|
|
-
|
|
|
- # if num.time>num.nucleation_time
|
|
|
- # BC_trans_scal[1].int = Dirichlet(val = num.concentration0[1])
|
|
|
- # else
|
|
|
- # BC_trans_scal[1].int = Neumann()
|
|
|
- # print("\n BC_trans_scal[1].int ", BC_trans_scal[1].int )
|
|
|
-
|
|
|
- # end
|
|
|
- # end
|
|
|
-
|
|
|
-
|
|
|
-
|
|
|
-
|
|
|
-
|
|
|
- end
|
|
|
-
|
|
|
- #TODO convection_Cdivu BC divergence
|
|
|
- #TODO check num.nb_transported_scalars>1
|
|
|
-
|
|
|
- if ((num.current_iter-1)%show_every == 0)
|
|
|
- # printstyled(color=:cyan, @sprintf "\n before scalar transport \n")
|
|
|
- # print_electrolysis_statistics(num,grid_p,phL)
|
|
|
-
|
|
|
- PDI_status = @ccall "libpdi".PDI_multi_expose("print_variables"::Cstring,
|
|
|
- "nstep"::Cstring, num.current_iter::Ref{Clonglong}, PDI_OUT::Cint,
|
|
|
- "time"::Cstring, num.time::Ref{Cdouble}, PDI_OUT::Cint,
|
|
|
- "u_1D"::Cstring, phL.uD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "v_1D"::Cstring, phL.vD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "p_1D"::Cstring, phL.pD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "levelset_p"::Cstring, grid_p.LS[num.iLSpdi].u::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "levelset_u"::Cstring, grid_u.LS[num.iLSpdi].u::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "levelset_v"::Cstring, grid_v.LS[num.iLSpdi].u::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "trans_scal_1DT"::Cstring, phL.trans_scalD'::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "phi_ele_1D"::Cstring, phL.phi_eleD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- C_NULL::Ptr{Cvoid})::Cint
|
|
|
- end
|
|
|
-
|
|
|
- # printstyled(color=:red, @sprintf "\n levelset: before scalar_transport\n")
|
|
|
- # println(grid_p.LS[1].geoL.dcap[1,1,:])
|
|
|
-
|
|
|
-
|
|
|
- #TODO better alloc
|
|
|
- # geo = [grid_p.LS[iLS].geoL for iLS in 1:num._nLS]
|
|
|
- # geo_u = [grid_u.LS[iLS].geoL for iLS in 1:num._nLS]
|
|
|
- # geo_v = [grid_v.LS[iLS].geoL for iLS in 1:num._nLS]
|
|
|
-
|
|
|
- # laps = set_matrices!(
|
|
|
- # num, grid_p, geo, grid_u, geo_u, grid_v, geo_v,
|
|
|
- # op.opC_pL, op.opC_uL, op.opC_vL,
|
|
|
- # periodic_x, periodic_y)
|
|
|
-
|
|
|
- if electrolysis_advection
|
|
|
-
|
|
|
- # if imposed_velocity == "zero" && num.current_iter ==16
|
|
|
- # num.ϵwall = 0.0
|
|
|
- # print("\n changed eps wall")
|
|
|
- # end
|
|
|
- #printstyled(color=:cyan, @sprintf "\n epsilon %.2e %.2e \n" num.ϵ num.ϵwall )
|
|
|
-
|
|
|
- update_all_ls_data(num, grid_p, grid_u, grid_v, BC_int, periodic_x, periodic_y)
|
|
|
- end
|
|
|
-
|
|
|
-
|
|
|
- # printstyled(color=:red, @sprintf "\n return before debug mem\n")
|
|
|
- # return
|
|
|
-
|
|
|
- # if num.nLS ==1 #TODO dev multiple levelsets
|
|
|
-
|
|
|
-
|
|
|
- #PDI (IO)
|
|
|
-
|
|
|
- if num.io_pdi>0
|
|
|
-
|
|
|
- #or permutedims(grid_p.LS[num.iLSpdi].geoL.dcap, (3, 2, 1)) (3, 1, 2)
|
|
|
- try
|
|
|
- # num.iLSpdi = 1 # TODO all grid_p.LS
|
|
|
- PDI_status = @ccall "libpdi".PDI_multi_expose("write_capacities"::Cstring,
|
|
|
- # "dcap"::Cstring, permutedims(grid_p.LS[num.iLSpdi].geoL.dcap, (3, 2, 1))::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "dcap_1"::Cstring, grid_p.LS[num.iLSpdi].geoL.dcap[:,:,1]::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "dcap_2"::Cstring, grid_p.LS[num.iLSpdi].geoL.dcap[:,:,2]::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "dcap_3"::Cstring, grid_p.LS[num.iLSpdi].geoL.dcap[:,:,3]::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "dcap_4"::Cstring, grid_p.LS[num.iLSpdi].geoL.dcap[:,:,4]::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- C_NULL::Ptr{Cvoid})::Cint
|
|
|
- catch error
|
|
|
- printstyled(color=:red, @sprintf "\n PDI error \n")
|
|
|
- print(error)
|
|
|
- printstyled(color=:red, @sprintf "\n PDI error \n")
|
|
|
- end
|
|
|
- end #if io_pdi
|
|
|
-
|
|
|
-
|
|
|
-
|
|
|
- # # #TODO better alloc
|
|
|
- # geo = [grid_p.LS[iLS].geoL for iLS in 1:num._nLS]
|
|
|
- # geo_u = [grid_u.LS[iLS].geoL for iLS in 1:num._nLS]
|
|
|
- # geo_v = [grid_v.LS[iLS].geoL for iLS in 1:num._nLS]
|
|
|
-
|
|
|
- # laps = set_matrices!(
|
|
|
- # num, grid_p, geo, grid_u, geo_u, grid_v, geo_v,
|
|
|
- # op.opC_pL, op.opC_uL, op.opC_vL,
|
|
|
- # periodic_x, periodic_y)
|
|
|
-
|
|
|
- #interface term in rhs comes from op.opC_TL
|
|
|
- #TODO variable CL
|
|
|
-
|
|
|
- PDI_status = @ccall "libpdi".PDI_multi_expose("write_iso"::Cstring,
|
|
|
- "nstep"::Cstring, num.current_iter::Ref{Clonglong}, PDI_OUT::Cint,
|
|
|
- "levelset_iso"::Cstring, grid_p.LS[num.iLSpdi].iso::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- C_NULL::Ptr{Cvoid})::Cint
|
|
|
-
|
|
|
- # printstyled(color=:magenta, @sprintf "\n Before scalar transport\n")
|
|
|
- # printstyled(color=:cyan, @sprintf "\n Before scalar transport\n")
|
|
|
- # printstyled(color=:green, @sprintf "\n Before scalar transport\n")
|
|
|
-
|
|
|
-
|
|
|
- scalar_transport!(num, grid_p, grid_u, grid_v,
|
|
|
- op.opC_TL, #op
|
|
|
- op.opL, #op_conv
|
|
|
- phL,
|
|
|
- BC_trans_scal,
|
|
|
- BC_int,
|
|
|
- Ascal,
|
|
|
- Bscal,
|
|
|
- all_CUTCT,
|
|
|
- rhs_scal,
|
|
|
- tmp_vec_p, #used to store a0 for rhs of interfacial value
|
|
|
- tmp_vec_u,
|
|
|
- tmp_vec_v,
|
|
|
- mass_transfer_rate,
|
|
|
- periodic_x,
|
|
|
- periodic_y,
|
|
|
- electrolysis_convection,
|
|
|
- ls_advection)
|
|
|
-
|
|
|
- print("\n num.stop_simulation after scalar ",num.stop_simulation)
|
|
|
-
|
|
|
- mask_1D = fnzeros(grid_p,num)
|
|
|
- compute_mask_1D!(num,grid_p,mask_1D)
|
|
|
-
|
|
|
- PDI_status = @ccall "libpdi".PDI_multi_expose("check_concentrations"::Cstring,
|
|
|
- "nstep"::Cstring, num.current_iter ::Ref{Clonglong}, PDI_OUT::Cint,
|
|
|
- "time"::Cstring, num.time::Ref{Cdouble}, PDI_OUT::Cint,
|
|
|
- "trans_scal_1DT"::Cstring, phL.trans_scalD'::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "mask_1D"::Cstring, mask_1D::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- C_NULL::Ptr{Cvoid})::Cint
|
|
|
-
|
|
|
- PDI_status = @ccall "libpdi".PDI_multi_expose("write_mask_one_phase"::Cstring,
|
|
|
- "nstep"::Cstring, num.current_iter ::Ref{Clonglong}, PDI_OUT::Cint,
|
|
|
- "mask_1D"::Cstring, mask_1D::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- C_NULL::Ptr{Cvoid})::Cint
|
|
|
-
|
|
|
- # if fix
|
|
|
- # print("\n BC_trans_scal ",BC_trans_scal[1])
|
|
|
-
|
|
|
- # print("\n BC_trans_scal ",BC_trans_scal[1].bottom.val)
|
|
|
-
|
|
|
-
|
|
|
- concentration_boundary_layer_width,averaged_electrode_concentration = compute_concentration_boundary_layer_width(num,grid_p,num.diffusion_coeff[1],
|
|
|
- num.current,num.saturation_concentration_H2,phL.trans_scalD[:,1],mask_1D,electrode_definition_function)
|
|
|
-
|
|
|
- # compute_concentration_boundary_layer_width(num,diffusion_coefficient,current,saturation_concentration,concentration_1D,mask_1D,func)
|
|
|
-
|
|
|
-
|
|
|
- sherwood = 2*num.R / concentration_boundary_layer_width #num.R : initial radius
|
|
|
-
|
|
|
- #TODO test rewrite diffusion + convection at interface with dot m
|
|
|
-
|
|
|
- sherwood_bubble = num.ambiant_pressure / (num.Ru * num.temperature0) * 2 * num.current_radius * (num.current_radius-num.previous_radius) /num.timestep_n / (num.diffusion_coeff[1] * (averaged_electrode_concentration - num.saturation_concentration_H2))
|
|
|
- # sherwood = compute_sherwood()
|
|
|
-
|
|
|
-
|
|
|
- PDI_status = @ccall "libpdi".PDI_multi_expose("write_sherwood"::Cstring,
|
|
|
- "nstep"::Cstring, num.current_iter ::Ref{Clonglong}, PDI_OUT::Cint,
|
|
|
- "concentration_boundary_layer_width"::Cstring, concentration_boundary_layer_width::Ref{Cdouble}, PDI_OUT::Cint,
|
|
|
- "sherwood"::Cstring, sherwood::Ref{Cdouble}, PDI_OUT::Cint,
|
|
|
- "sherwood_bubble"::Cstring, sherwood_bubble::Ref{Cdouble}, PDI_OUT::Cint,
|
|
|
- C_NULL::Ptr{Cvoid})::Cint
|
|
|
-
|
|
|
- # PDI_status = @ccall "libpdi".PDI_multi_expose("check_concentrations"::Cstring,
|
|
|
-
|
|
|
- # scalar_transport!(BC_trans_scal, num, grid_p, , grid_p.LS[1].geoL, phL, num.concentration0,
|
|
|
- # grid_p.LS[1].MIXED, grid_p.LS[1].geoL.projection, op.opL, grid_u, grid_u.LS[1].geoL, grid_v, grid_v.LS[1].geoL,
|
|
|
- # periodic_x, periodic_y, electrolysis_convection, ls_advection, BC_int, num.diffusion_coeff,Ascal,Bscal,all_CUTCT,rhs_scal)
|
|
|
-
|
|
|
- # else
|
|
|
- # printstyled(color=:red, @sprintf "\n TODO multiple LS \n" )
|
|
|
-
|
|
|
- # end
|
|
|
-
|
|
|
- # scalar_transport_2!(BC_trans_scal, num, grid_p, op.opC_TL, grid_p.LS[1].geoL, phL, num.concentration0,
|
|
|
- # grid_p.LS[1].MIXED, grid_p.LS[1].geoL.projection, op.opL, grid_u, grid_u.LS[1].geoL, grid_v, grid_v.LS[1].geoL,
|
|
|
- # periodic_x, periodic_y, electrolysis_convection, true, BC_int, num.diffusion_coeff)
|
|
|
-
|
|
|
- if ((num.current_iter-1)%show_every == 0)
|
|
|
- # printstyled(color=:cyan, @sprintf "\n after scalar transport \n")
|
|
|
- # print_electrolysis_statistics(num,grid_p,phL)
|
|
|
-
|
|
|
- PDI_status = @ccall "libpdi".PDI_multi_expose("print_variables"::Cstring,
|
|
|
- "nstep"::Cstring, num.current_iter ::Ref{Clonglong}, PDI_OUT::Cint,
|
|
|
- "time"::Cstring, num.time::Ref{Cdouble}, PDI_OUT::Cint,
|
|
|
- "u_1D"::Cstring, phL.uD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "v_1D"::Cstring, phL.vD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "p_1D"::Cstring, phL.pD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "levelset_p"::Cstring, grid_p.LS[num.iLSpdi].u::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "levelset_u"::Cstring, grid_u.LS[num.iLSpdi].u::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "levelset_v"::Cstring, grid_v.LS[num.iLSpdi].u::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "trans_scal_1DT"::Cstring, phL.trans_scalD'::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "phi_ele_1D"::Cstring, phL.phi_eleD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- C_NULL::Ptr{Cvoid})::Cint
|
|
|
-
|
|
|
-
|
|
|
- # scal_error=0.0
|
|
|
- # for iscal in 1:num.nb_transported_scalars
|
|
|
- # # printstyled(color=:cyan, @sprintf "\n after scalar transport %.2i %.2e \n" iscal maximum(abs.(phL.trans_scal[:,:,iscal])))
|
|
|
- # # printstyled(color=:cyan, @sprintf "\n after scalar transport %.2i %.2e \n" iscal maximum(abs.(phL.trans_scalD[:,iscal])))
|
|
|
- # # scal_error_bulk = maximum(abs.(phL.trans_scal[:,:,iscal].-num.concentration0[iscal])./num.concentration0[iscal])
|
|
|
- # # scal_error_border = maximum(abs.(vecb(phL.trans_scalD[:,iscal],grid_p).-num.concentration0[iscal])./num.concentration0[iscal])
|
|
|
- # # scal_error = max(scal_error_bulk,scal_error_border,scal_error)
|
|
|
- # end
|
|
|
-
|
|
|
- # printstyled(color=:cyan, @sprintf "\n error after scalar transport %.2e num.CFL %.2e\n" scal_error num.v_inlet*num.timestep_n/grid_p.dx[1,1])
|
|
|
-
|
|
|
-
|
|
|
- end
|
|
|
-
|
|
|
- if imposed_velocity != "none" && num.debug== "scalar_testing"
|
|
|
- scal_error=0.0
|
|
|
- for iscal in 1:num.nb_transported_scalars
|
|
|
-
|
|
|
- # print("\n maximum ",maximum(phL.trans_scalD[:,iscal]), )
|
|
|
- # printstyled(color=:cyan, @sprintf "\n error after scalar transport max %.2e min %.2e\n" maximum(phL.trans_scalD[:,iscal]) minimum(phL.trans_scalD[:,iscal]))
|
|
|
-
|
|
|
- scal_error_bulk = maximum(abs.(phL.trans_scal[:,:,iscal].-num.concentration0[iscal])./num.concentration0[iscal])
|
|
|
- scal_error_border = maximum(abs.(vecb(phL.trans_scalD[:,iscal],grid_p).-num.concentration0[iscal])./num.concentration0[iscal])
|
|
|
- scal_error = max(scal_error_bulk,scal_error_border,scal_error)
|
|
|
-
|
|
|
- end
|
|
|
-
|
|
|
- printstyled(color=:cyan, @sprintf "\n error after scalar transport %.2e num.CFL %.2e\n" scal_error num.v_inlet*num.timestep_n/grid_p.dx[1,1])
|
|
|
-
|
|
|
- # printstyled(color=:red, @sprintf "\n Poiseuille \n")
|
|
|
-
|
|
|
- # # Check the velocity field before the scalar transport
|
|
|
- # test_Poiseuille(num,phL.vD,grid_v)
|
|
|
-
|
|
|
- # printstyled(color=:cyan, @sprintf "\n pressure min %.2e max %.2e\n" minimum(phL.p[1,:]) maximum(phL.p[1,:]))
|
|
|
-
|
|
|
- # printstyled(color=:cyan, @sprintf "\n pressure min %.2e max %.2e\n" minimum(phL.p[end,:]) maximum(phL.p[end,:]))
|
|
|
-
|
|
|
- # printstyled(color=:cyan, @sprintf "\n pressure min %.2e max %.2e\n" BC_pL.bottom.val BC_pL.top.val )
|
|
|
-
|
|
|
- # compute_grad_p!(num,grid_p, grid_u, grid_v, phL.pD, op.opC_pL, op.opC_uL, op.opC_vL)
|
|
|
-
|
|
|
- end
|
|
|
-
|
|
|
- end #num.nb_transported_scalars>0
|
|
|
-
|
|
|
- #endregion Scalar transport
|
|
|
-
|
|
|
- # if imposed_velocity =="none"
|
|
|
- # printstyled(color=:red, @sprintf "\n after scalar transport \n")
|
|
|
-
|
|
|
- # # Check the velocity field before the scalar transport
|
|
|
- # test_Poiseuille(num,phL,grid_v)
|
|
|
-
|
|
|
- # end
|
|
|
-
|
|
|
-
|
|
|
- end #if electrolysis_liquid_phase
|
|
|
- end #if electrolysis
|
|
|
- #endregion Electrolysis
|
|
|
-
|
|
|
-
|
|
|
-
|
|
|
-
|
|
|
- #PDI (IO)
|
|
|
- if num.io_pdi>0
|
|
|
-
|
|
|
- try
|
|
|
- # printstyled(color=:red, @sprintf "\n PDI test \n" )
|
|
|
-
|
|
|
-
|
|
|
- # phi_array=phL.phi_ele #do not transpose since python row major
|
|
|
-
|
|
|
- # Compute electrical current, interpolate velocity on scalar grid_p
|
|
|
- # if num.electrical_potential>0 compute_grad_phi_ele!(num, grid_p, grid_u, grid_v, phL, phS, op.opC_pL, op.opC_pS) #TODO current
|
|
|
-
|
|
|
- # if electrolysis && num.nb_transported_scalars>1
|
|
|
- # if heat
|
|
|
- # elec_cond = 2*num.Faraday^2 .*phL.trans_scal[:,:,2].*num.diffusion_coeff[2]./(num.Ru.*phL.T) #phL.T
|
|
|
- # else
|
|
|
- # elec_cond = 2*num.Faraday^2 .*phL.trans_scal[:,:,2].*num.diffusion_coeff[2]./(num.Ru*num.temperature0)
|
|
|
- # end
|
|
|
- # else
|
|
|
- # elec_cond = ones(grid_p)
|
|
|
- # printstyled(color=:green, @sprintf "\n conductivity one")
|
|
|
-
|
|
|
- # end
|
|
|
-
|
|
|
- # phL.i_current_mag .*= elec_cond # i=-κ∇ϕ here magnitude
|
|
|
-
|
|
|
-
|
|
|
- #store in us, vs instead of Eus, Evs
|
|
|
- # interpolate_grid_liquid!(grid_p,grid_u,grid_v,phL.Eu, phL.Ev,tmp_vec_p,tmp_vec_p0)
|
|
|
-
|
|
|
- # printstyled(color=:red, @sprintf "\n test i current\n" )
|
|
|
- # print("\n phi ",phL.phi_ele[64,127]," ",phL.phi_ele[64,128]," ",(phL.phi_ele[64,128]-phL.phi_ele[64,127])/grid_p.dx[64,128]," ",(phL.phi_ele[64,128]-phL.phi_ele[64,127])/grid_p.dx[64,127]*elec_cond[64,127], " cond ", elec_cond[64,127]," \n")
|
|
|
-
|
|
|
-
|
|
|
- # tmp_vec_p .*= elec_cond
|
|
|
- # tmp_vec_p0 .*= elec_cond
|
|
|
-
|
|
|
-
|
|
|
- # @ccall "libpdi".PDI_multi_expose("write_data_elec"::Cstring,
|
|
|
- # "i_current_x"::Cstring, tmp_vec_p::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # "i_current_y"::Cstring, tmp_vec_p0::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # "i_current_mag"::Cstring, phL.i_current_mag::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # "phi_ele_1D"::Cstring, phL.phi_eleD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # C_NULL::Ptr{Cvoid})::Cint
|
|
|
-
|
|
|
- interpolate_staggered_u_v_to_scalar_grid_one_fluid_or_one_phase!(num,grid_p,grid_u,grid_v,phL.u,phL.v,tmp_vec_p,tmp_vec_p0)
|
|
|
-
|
|
|
- # num.iLSpdi = 1 # TODO all grid_p.LS
|
|
|
-
|
|
|
- # Exposing data to PDI for IO
|
|
|
- # if writing "D" array (bulk, interface, border), add "_1D" to the name
|
|
|
-
|
|
|
- PDI_status = @ccall "libpdi".PDI_multi_expose("write_data"::Cstring,
|
|
|
- "nstep"::Cstring, num.current_iter::Ref{Clonglong}, PDI_OUT::Cint,
|
|
|
- "time"::Cstring, num.time::Ref{Cdouble}, PDI_OUT::Cint,
|
|
|
- "timestep"::Cstring, num.timestep_n::Ref{Cdouble}, PDI_OUT::Cint,
|
|
|
- "nx"::Cstring, grid_p.nx::Ref{Clonglong}, PDI_OUT::Cint,
|
|
|
- "ny"::Cstring, grid_p.ny::Ref{Clonglong}, PDI_OUT::Cint,
|
|
|
- "u_1D"::Cstring, phL.uD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "v_1D"::Cstring, phL.vD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "p_1D"::Cstring, phL.pD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "levelset_p"::Cstring, grid_p.LS[num.iLSpdi].u::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "levelset_u"::Cstring, grid_u.LS[num.iLSpdi].u::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "levelset_v"::Cstring, grid_v.LS[num.iLSpdi].u::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "trans_scal_1DT"::Cstring, phL.trans_scalD'::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "phi_ele_1D"::Cstring, phL.phi_eleD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "velocity_x"::Cstring, tmp_vec_p::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "velocity_y"::Cstring, tmp_vec_p0::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "radius"::Cstring, num.current_radius::Ref{Cdouble}, PDI_OUT::Cint,
|
|
|
- C_NULL::Ptr{Cvoid})::Cint
|
|
|
-
|
|
|
- # if num.nb_transported_scalars>0
|
|
|
- # PDI_status = @ccall "libpdi".PDI_multi_expose("write_data_species"::Cstring,
|
|
|
- # "nstep"::Cstring, num.current_iter::Ref{Clonglong}, PDI_OUT::Cint,
|
|
|
- # "trans_scal_1DT"::Cstring, phL.trans_scalD'::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # C_NULL::Ptr{Cvoid})::Cint
|
|
|
- # end
|
|
|
-
|
|
|
-
|
|
|
-
|
|
|
- #TODO debug with volume fraction
|
|
|
-
|
|
|
- # A = zeros(gv.ny, gv.nx+1)
|
|
|
- # for jplot in 1:gv.ny
|
|
|
- # for iplot in 1:gv.nx+1
|
|
|
- # II = CartesianIndex(jplot, iplot) #(id_y, id_x)
|
|
|
- # pII = lexicographic(II, grid_p.ny + 1)
|
|
|
- # A[jplot,iplot] = 1 ./ op.opC_vL.iMx.diag[pII]
|
|
|
- # end
|
|
|
- # end
|
|
|
-
|
|
|
-
|
|
|
- # print("\n after write \n ")
|
|
|
-
|
|
|
- # printstyled(color=:red, @sprintf "\n PDI test end\n" )
|
|
|
-
|
|
|
- catch error
|
|
|
- printstyled(color=:red, @sprintf "\n PDI error \n")
|
|
|
- print(error)
|
|
|
- printstyled(color=:red, @sprintf "\n PDI error \n")
|
|
|
- end
|
|
|
- end #if io_pdi
|
|
|
-
|
|
|
- #region compute mass transfer
|
|
|
- if num.solve_Navier_Stokes_liquid_phase == 1 && num.phase_change_method >0
|
|
|
- num.previous_radius = num.current_radius
|
|
|
- # num.status =
|
|
|
- compute_mass_transfer_rate_main!(num, grid_p, grid_u, grid_v, op, phL, phS, BC_int, electrolysis, electrolysis_phase_change_case,
|
|
|
- periodic_x, periodic_y, λ, Vmean, num.iLSpdi, mode_2d, show_every,
|
|
|
- mass_transfer_rate,mass_transfer_rate_vec1,
|
|
|
- mass_transfer_rate_vecb,mass_transfer_rate_veci, mass_transfer_rate_redistributed, tmp_vec_p, tmp_vec_p0, tmp_vec_p1,
|
|
|
- nb_gaz_acceptors, volume_fraction, interface_length)
|
|
|
- #endregion compute mass transfer
|
|
|
- end
|
|
|
-
|
|
|
- #region Navier-Stokes
|
|
|
- if num.solve_Navier_Stokes_liquid_phase == 1
|
|
|
- if num.time < num.nucleation_time
|
|
|
- printstyled(color=:red, @sprintf "\n Navier-Stokes not solved")
|
|
|
-
|
|
|
- navier_stokes = false
|
|
|
- else
|
|
|
- printstyled(color=:red, @sprintf "\n Navier-Stokes solved")
|
|
|
- navier_stokes = true
|
|
|
- end
|
|
|
-
|
|
|
- if navier_stokes
|
|
|
-
|
|
|
-
|
|
|
- # Adapt cell volume W for gradients
|
|
|
- # cf 4/3 factor for Laplacian
|
|
|
- if num.laplacian == 1
|
|
|
- AvLcopy = copy(AvL) #does not work if reset A after operations in compute_divergence!
|
|
|
- Lvm1_L,bc_Lvm1_L,bc_Lvm1_b_L=compute_divergence!(num,
|
|
|
- # grid_p,
|
|
|
- # grid_u,
|
|
|
- grid_v,
|
|
|
- op.opC_vL,
|
|
|
- AvLcopy,
|
|
|
- # rhs_scal,
|
|
|
- # tmp_vec_p, #a0
|
|
|
- tmp_vec_1D_v0,
|
|
|
- tmp_vec_1D_v,
|
|
|
- Lvm1_L,
|
|
|
- bc_Lvm1_L,
|
|
|
- bc_Lvm1_b_L
|
|
|
- # tmp_vec_u0,
|
|
|
- # tmp_vec_v0,
|
|
|
- # tmp_vec_1D,
|
|
|
- # ls_advection
|
|
|
- )
|
|
|
- print("\nbefore pressure")
|
|
|
- II = CartesianIndex(div(grid_v.ny,2),1)
|
|
|
- pII = lexicographic(II,grid_v.ny)
|
|
|
- print("\nLv[pII,:] ",Lvm1_L[pII,:])
|
|
|
- print("\bc_Lvm1_b[pII,:] ",bc_Lvm1_b_L[pII,:])
|
|
|
-
|
|
|
- ApLcopy = copy(Ascal)
|
|
|
- Lpm1_L,bc_Lpm1_L,bc_Lpm1_b_L=compute_divergence!(num,
|
|
|
- # grid_p,
|
|
|
- # grid_u,
|
|
|
- grid_p,
|
|
|
- op.opC_pL,
|
|
|
- ApLcopy,
|
|
|
- # rhs_scal,
|
|
|
- # tmp_vec_p, #a0
|
|
|
- tmp_vec_1D_p0,
|
|
|
- tmp_vec_1D_p,
|
|
|
- Lpm1_L,
|
|
|
- bc_Lpm1_L,
|
|
|
- bc_Lpm1_b_L
|
|
|
- # tmp_vec_u0,
|
|
|
- # tmp_vec_v0,
|
|
|
- # tmp_vec_1D,
|
|
|
- # ls_advection
|
|
|
- )
|
|
|
-
|
|
|
- end
|
|
|
-
|
|
|
- if num.pressure_velocity_coupling == 3
|
|
|
-
|
|
|
- # num.iLSpdi = 1 # TODO all grid_p.LS
|
|
|
- PDI_status = @ccall "libpdi".PDI_multi_expose("print_capacities"::Cstring,
|
|
|
- # "dcap"::Cstring, permutedims(grid_p.LS[num.iLSpdi].geoL.dcap, (3, 2, 1))::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "dcap_1"::Cstring, grid_u.LS[num.iLSpdi].geoL.dcap[:,:,1]::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "dcap_2"::Cstring, grid_u.LS[num.iLSpdi].geoL.dcap[:,:,2]::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "dcap_3"::Cstring, grid_u.LS[num.iLSpdi].geoL.dcap[:,:,3]::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "dcap_4"::Cstring, grid_u.LS[num.iLSpdi].geoL.dcap[:,:,4]::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- C_NULL::Ptr{Cvoid})::Cint
|
|
|
-
|
|
|
- PDI_status = @ccall "libpdi".PDI_multi_expose("print_capacities"::Cstring,
|
|
|
- # "dcap"::Cstring, permutedims(grid_p.LS[num.iLSpdi].geoL.dcap, (3, 2, 1))::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "dcap_1"::Cstring, grid_v.LS[num.iLSpdi].geoL.dcap[:,:,1]::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "dcap_2"::Cstring, grid_v.LS[num.iLSpdi].geoL.dcap[:,:,2]::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "dcap_3"::Cstring, grid_v.LS[num.iLSpdi].geoL.dcap[:,:,3]::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "dcap_4"::Cstring, grid_v.LS[num.iLSpdi].geoL.dcap[:,:,4]::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- C_NULL::Ptr{Cvoid})::Cint
|
|
|
-
|
|
|
- # print("\n cap_1 ",grid_u.LS[num.iLSpdi].geoL.dcap[1,:,1])
|
|
|
- # print("\n cap_1 ",grid_u.LS[num.iLSpdi].geoL.dcap[2,:,1])
|
|
|
-
|
|
|
- # print("\n cap_1 ",grid_v.LS[num.iLSpdi].geoL.dcap[1,:,1])
|
|
|
- # print("\n cap_1 ",grid_v.LS[num.iLSpdi].geoL.dcap[2,:,1])
|
|
|
-
|
|
|
- # print("\n cap_2 ",grid_v.LS[num.iLSpdi].geoL.dcap[1,:,2])
|
|
|
- # print("\n cap_2 ",grid_v.LS[num.iLSpdi].geoL.dcap[2,:,2])
|
|
|
-
|
|
|
- # print("\n cap_3 ",grid_v.LS[num.iLSpdi].geoL.dcap[1,:,3])
|
|
|
- # print("\n cap_3 ",grid_v.LS[num.iLSpdi].geoL.dcap[2,:,3])
|
|
|
-
|
|
|
- # for u grid_p
|
|
|
- # cap 1 same
|
|
|
- # cap_1 [0.0, 1.0e-5, 1.0e-5, 1.0e-5, 1.0e-5, 1.0e-5, 1.0e-5, 1.0e-5, 1.0e-5, 1.0e-5, 1.0e-5]
|
|
|
-
|
|
|
-
|
|
|
-
|
|
|
-
|
|
|
- end
|
|
|
-
|
|
|
-
|
|
|
- # if !advection
|
|
|
- # @time no_slip_condition!(num, grid_p, grid_u, grid_u.LS[1], grid_v, grid_v.LS[1], periodic_x, periodic_y)
|
|
|
- # # grid_u.V .= num.Δ / (1 * num.timestep_n)
|
|
|
- # # grid_v.V .= 0.0
|
|
|
- # end
|
|
|
-
|
|
|
- # Pressure-velocity coupling
|
|
|
-
|
|
|
- if (num.one_fluid_model == 1 && num.solve_Navier_Stokes_liquid_phase == 1)
|
|
|
-
|
|
|
- # interpolate_staggered_u_v_to_scalar_grid_one_fluid_or_one_phase!(num,grid_p,grid_u,grid_v,phL.u,phL.v,tmp_vec_p,tmp_vec_p0)
|
|
|
-
|
|
|
- tmp_vec_p .= 0.0
|
|
|
- tmp_vec_p0 .= 0.0
|
|
|
-
|
|
|
- for j = 1:grid_p.ny
|
|
|
- for i = 1:grid_p.nx
|
|
|
- tmp_vec_p[j,i] =(phL.u[j,i]+phL.u[j,i+1])/2
|
|
|
- tmp_vec_p0[j,i]=(phL.v[j,i]+phL.v[j+1,i])/2
|
|
|
- end
|
|
|
- end
|
|
|
-
|
|
|
-
|
|
|
- update_one_fluid_density_viscosity(num,grid_p,grid_u,grid_v,volume_fraction,levelset_one_fluid,rho_one_fluid,
|
|
|
- rho_one_fluid_u,rho_one_fluid_v,mu_one_fluid,tmp_vec_p0)
|
|
|
-
|
|
|
- total_interface_length = compute_interface_length!(num, grid_p, 1, interface_length)
|
|
|
- # MIXED =
|
|
|
-
|
|
|
- nb_levelsets = num.nLS
|
|
|
- num.nLS = 0 #1 #deactivate cut-cell for one-fluid model
|
|
|
-
|
|
|
- #region deactivate LS
|
|
|
- empty_capacities = vcat(zeros(7), zeros(4))
|
|
|
- full_capacities = vcat(ones(7), 0.5.*ones(4))
|
|
|
-
|
|
|
- for grid_iter in [grid_p,grid_u,grid_v]
|
|
|
-
|
|
|
- grid_iter.LS[1].u .= 1.0
|
|
|
- grid_iter.LS[end].u .= 1.0
|
|
|
-
|
|
|
-
|
|
|
- # for j in 1:grid_iter.ny
|
|
|
- # for i in 1:grid_iter.nx
|
|
|
- # II = CartesianIndex(j,i)
|
|
|
- # # grid_iter.LS[end].geoL.cap[II,:] .= full_capacities
|
|
|
- # # grid_iter.LS[end].geoS.cap[II,:] .= empty_capacities
|
|
|
- # end
|
|
|
- # end
|
|
|
- end
|
|
|
-
|
|
|
- # grid_u.LS[end].geoL.cap[:,:,:] .= full_capacities
|
|
|
- # grid_u.LS[end].geoS.cap[:,:,:] .= empty_capacities
|
|
|
-
|
|
|
- # grid_v.LS[end].geoL.cap[:,:,:] .= full_capacities
|
|
|
- # grid_v.LS[end].geoS.cap[:,:,:] .= empty_capacities
|
|
|
-
|
|
|
- #endregion deactivate LS
|
|
|
-
|
|
|
- NB_indices = update_all_ls_data(num, grid_p, grid_u, grid_v, BC_int, periodic_x, periodic_y,true,true)
|
|
|
-
|
|
|
-
|
|
|
- geoL = [grid_p.LS[iLS].geoL for iLS in 1:num._nLS]
|
|
|
- geo_uL = [grid_u.LS[iLS].geoL for iLS in 1:num._nLS]
|
|
|
- geo_vL = [grid_v.LS[iLS].geoL for iLS in 1:num._nLS]
|
|
|
-
|
|
|
- #region reset centroids
|
|
|
- for grid_iter in [grid_p,grid_u,grid_v]
|
|
|
- for II in grid_iter.ind.inside
|
|
|
- grid_iter.LS[1].geoL.centroid[II] = Point(0.0,0.0)
|
|
|
- end
|
|
|
- end
|
|
|
- #endregion reset centroids
|
|
|
- #TODO reactivate centroids ?
|
|
|
- # TODO update density
|
|
|
-
|
|
|
- # if num.pressure_velocity_coupling == 0
|
|
|
-
|
|
|
- #region update LS
|
|
|
-
|
|
|
- #endregion update LS
|
|
|
-
|
|
|
-
|
|
|
-
|
|
|
-
|
|
|
-
|
|
|
-
|
|
|
- #pbm mass_transfer_rateL
|
|
|
- PDI_status = @ccall "libpdi".PDI_multi_expose("check_mass_transfer_rate_NS"::Cstring,
|
|
|
- # "conservation"::Cstring, conservation::Ref{Cdouble}, PDI_OUT::Cint,
|
|
|
- "mass_transfer_rate"::Cstring, mass_transfer_rate::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # "p_1D"::Cstring, phL.pD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- C_NULL::Ptr{Cvoid})::Cint
|
|
|
-
|
|
|
-
|
|
|
-
|
|
|
- #region update LS for convection (bool=true)+ one fluid
|
|
|
-
|
|
|
- # At every iteration, update_all_ls_data is called twice, once inside run.jl
|
|
|
- # and another one (if there's advection of the levelset) inside set_heat!.
|
|
|
- # The difference between both is a flag as last argument, inside run.jl is implicitly defined
|
|
|
- # as true and inside set_heat! is false. If you're calling your version of set_heat! several times,
|
|
|
- # then you're calling the version with the flag set to false, but for the convective term it has to be set to true.
|
|
|
- # The flag=true, the capacities are set for the convection, the flag=false they are set for the other operators
|
|
|
-
|
|
|
- if advection
|
|
|
- # update_all_ls_data(num, grid_p, grid_u, grid_v, BC_int, periodic_x, periodic_y, true)
|
|
|
- update_all_ls_data(num, grid_p, grid_u, grid_v, BC_int, periodic_x, periodic_y, true,true)
|
|
|
- if num.convection == 0
|
|
|
- # op.opL is op_conv
|
|
|
- set_convection_preallocated!(num, grid_p, geoL[end], grid_u, grid_u.LS, grid_v, grid_v.LS, phL.u, phL.v, op.opL,
|
|
|
- phL, BC_uL, BC_vL,op.opC_pL, op.opC_uL, op.opC_vL,
|
|
|
- velocity_and_BC_convection_u_x ,
|
|
|
- velocity_and_BC_convection_u_y ,
|
|
|
- velocity_and_BC_convection_v_x ,
|
|
|
- velocity_and_BC_convection_v_y)
|
|
|
-
|
|
|
- # Mm1_L .= op.opC_pL.M
|
|
|
- # Mum1_L .= op.opC_uL.M
|
|
|
- # Mvm1_L .= op.opC_vL.M
|
|
|
- else
|
|
|
-
|
|
|
-
|
|
|
- end
|
|
|
-
|
|
|
-
|
|
|
- end
|
|
|
- #endregion update LS for convection (bool=true)+ one fluid
|
|
|
-
|
|
|
-
|
|
|
- #region update LS for other operators than convection (bool=false)+ one fluid
|
|
|
- # update_all_ls_data(num, grid_p, grid_u, grid_v, BC_int, periodic_x, periodic_y, false)
|
|
|
- update_all_ls_data(num, grid_p, grid_u, grid_v, BC_int, periodic_x, periodic_y, false,true)
|
|
|
-
|
|
|
- #endregion
|
|
|
-
|
|
|
-
|
|
|
-
|
|
|
- if (num.one_fluid_model == 1 && num.solve_Navier_Stokes_liquid_phase == 1) && (num.solve_Navier_Stokes == 1)
|
|
|
- if num.activate_interface == 0
|
|
|
- volumic_surface_tension_u .= 0.0
|
|
|
- volumic_surface_tension_v .= 0.0
|
|
|
- else
|
|
|
- if num.surface_tension == 0
|
|
|
- compute_surface_tension_VOF!(num,grid_p, grid_u, grid_v, op.opC_pL, op.opC_uL, op.opC_vL,
|
|
|
- volume_fraction,levelset_one_fluid,volumic_surface_tension_u,volumic_surface_tension_v,tmp_vec_p,tmp_vec_p0)
|
|
|
- elseif num.surface_tension == 1
|
|
|
- # compute_surface_tension_LS!(num,grid_p, grid_u, grid_v, opC_p, opC_u, opC_v,
|
|
|
- # volume_fraction,levelset_one_fluid,volumic_surface_tension_u,volumic_surface_tension_v,tmp_vec_p,tmp_vec_p0)
|
|
|
- compute_surface_tension_LS!(num,grid_p, grid_u, grid_v, op.opC_pL, op.opC_uL, op.opC_vL,
|
|
|
- volume_fraction,levelset_one_fluid,volumic_surface_tension_u,volumic_surface_tension_v,tmp_vec_p,tmp_vec_p0,
|
|
|
- levelset_1D, levelset_heavyside_2D,
|
|
|
- tmp_vec_u0,tmp_vec_v0, #normal_and_dirac_u, normal_and_dirac_v,
|
|
|
- tmp_vec_u1,tmp_vec_v1,#normal_u, normal_v,
|
|
|
- tmp_vec_u,tmp_vec_v,#curvature_u, curvature_v
|
|
|
- )
|
|
|
- # elseif num. TODO HF
|
|
|
- end
|
|
|
-
|
|
|
- end
|
|
|
- end
|
|
|
-
|
|
|
- PDI_status = @ccall "libpdi".PDI_multi_expose("write_one_fluid_surface_tension_concise"::Cstring,
|
|
|
- "nstep"::Cstring, num.current_iter ::Ref{Clonglong}, PDI_OUT::Cint,
|
|
|
- "volumic_surface_tension_u"::Cstring, volumic_surface_tension_u::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "volumic_surface_tension_v"::Cstring, volumic_surface_tension_v::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- C_NULL::Ptr{Cvoid})::Cint
|
|
|
-
|
|
|
-
|
|
|
- #BC TODO
|
|
|
-
|
|
|
- # @error("\n check temporal terms")
|
|
|
-
|
|
|
- # Mum1_L is put in B matrix that multiplies v
|
|
|
- # print("\n advection ", ns_advection, " adv ",advection)
|
|
|
- one_fluid_NS_ls_advection = true #update matrix
|
|
|
-
|
|
|
-
|
|
|
- Lpm1_L, bc_Lpm1_L, bc_Lpm1_b_L, Lum1_L, bc_Lum1_L, bc_Lum1_b_L,
|
|
|
- Lvm1_L, bc_Lvm1_L, bc_Lvm1_b_L, Mm1_L, Mum1_L, Mvm1_L, Cum1L, Cvm1L = solve_one_fluid_NS!(
|
|
|
- time_scheme, BC_int,
|
|
|
- num, grid_p, geoL, grid_u, geo_uL, grid_v, geo_vL, phL,
|
|
|
- BC_uL, BC_vL, BC_pL,
|
|
|
- op.opC_pL, op.opC_uL, op.opC_vL, op.opL,
|
|
|
- # op.opC_TL,
|
|
|
- AuL, BuL, AvL, BvL, AϕL, AuvL, BuvL,rhs_uv,
|
|
|
- Lpm1_L, bc_Lpm1_L, bc_Lpm1_b_L, Lum1_L, bc_Lum1_L, bc_Lum1_b_L, Lvm1_L, bc_Lvm1_L, bc_Lvm1_b_L,
|
|
|
- Cum1L, Cvm1L, Mum1_L, Mvm1_L,
|
|
|
- periodic_x, periodic_y, ns_advection, one_fluid_NS_ls_advection, num.current_iter, Ra, navier,
|
|
|
- volume_fraction,
|
|
|
- levelset_one_fluid,
|
|
|
- rho_one_fluid,
|
|
|
- mu_one_fluid,
|
|
|
- rho_one_fluid_u,
|
|
|
- # mu_one_fluid_u,
|
|
|
- rho_one_fluid_v,
|
|
|
- # mu_one_fluid_v,
|
|
|
- volumic_surface_tension_u,
|
|
|
- volumic_surface_tension_v,
|
|
|
- convection_u,convection_v,
|
|
|
- viscosity_coeff_for_du_dx ,
|
|
|
- viscosity_coeff_for_du_dy ,
|
|
|
- viscosity_coeff_for_dv_dx ,
|
|
|
- viscosity_coeff_for_dv_dy,
|
|
|
- # velocity_and_BC_convection_u_x ,
|
|
|
- # velocity_and_BC_convection_u_y ,
|
|
|
- # velocity_and_BC_convection_v_x ,
|
|
|
- # velocity_and_BC_convection_v_y ,
|
|
|
- tmp_vec_p,
|
|
|
- tmp_vec_p0,
|
|
|
- rhs_phi,
|
|
|
- pres_free_surfaceL,jump_mass_transfer_rateL,mass_transfer_rate )
|
|
|
-
|
|
|
- #region ciprianoMulticomponentDropletEvaporation2024
|
|
|
- if extend_liquid_velocity
|
|
|
-
|
|
|
- # num.phase_change_currently_activated == 1
|
|
|
- phase_change_currently_activated = 0
|
|
|
-
|
|
|
- # no need to recompute Lpm1_L, bc_Lpm1_L, bc_Lpm1_b_L, Lum1_L, bc_Lum1_L, bc_Lum1_b_L,
|
|
|
- # Lvm1_L, bc_Lvm1_L, bc_Lvm1_b_L
|
|
|
-
|
|
|
- # Mm1_L_ext_vel, Mum1_L_ext_vel, Mvm1_L_ext_vel same as Mm1_L, Mum1_L, Mvm1_L
|
|
|
- # because :
|
|
|
- # Mm1_L = copy(op.opC_pL.M)
|
|
|
- # Mum1_L = copy(op.opC_uL.M)
|
|
|
- # Mvm1_L = copy(op.opC_vL.M)
|
|
|
-
|
|
|
- Lpm1_L, bc_Lpm1_L, bc_Lpm1_b_L, Lum1_L, bc_Lum1_L, bc_Lum1_b_L,
|
|
|
- Lvm1_L, bc_Lvm1_L, bc_Lvm1_b_L,
|
|
|
- Mm1_L, Mum1_L, Mvm1_L,
|
|
|
- Cum1L_ext_vel, Cvm1L_ext_vel = solve_one_fluid_NS_no_phase!(
|
|
|
- time_scheme, BC_int,
|
|
|
- num, grid_p, geoL, grid_u, geo_uL, grid_v, geo_vL,
|
|
|
- # phL,
|
|
|
- p_ext_vel, pD_ext_vel, phi_ext_vel, u_ext_vel, v_ext_vel, u_predictionD_ext_vel, v_predictionD_ext_vel, uD_ext_vel, vD_ext_vel, u_prediction_ext_vel, v_prediction_ext_vel, uT_ext_vel,
|
|
|
- pres_grad_x, pres_grad_y,
|
|
|
- phase_change_currently_activated,
|
|
|
- BC_uL, BC_vL, BC_pL,
|
|
|
- op.opC_pL, op.opC_uL, op.opC_vL, op.opL,
|
|
|
- # op.opC_TL,
|
|
|
- AuL, BuL, AvL, BvL, AϕL, AuvL, BuvL,rhs_uv,
|
|
|
- Lpm1_L, bc_Lpm1_L, bc_Lpm1_b_L,
|
|
|
- Lum1_L, bc_Lum1_L, bc_Lum1_b_L,
|
|
|
- Lvm1_L, bc_Lvm1_L, bc_Lvm1_b_L,
|
|
|
- Cum1L_ext_vel, Cvm1L_ext_vel,
|
|
|
- Mum1_L, Mvm1_L,
|
|
|
- periodic_x, periodic_y, ns_advection, advection, num.current_iter, Ra, navier,
|
|
|
- volume_fraction,
|
|
|
- levelset_one_fluid,
|
|
|
- rho_one_fluid,
|
|
|
- mu_one_fluid,
|
|
|
- rho_one_fluid_u,
|
|
|
- rho_one_fluid_v,
|
|
|
- volumic_surface_tension_u,
|
|
|
- volumic_surface_tension_v,
|
|
|
- convection_u,convection_v,
|
|
|
- viscosity_coeff_for_du_dx ,
|
|
|
- viscosity_coeff_for_du_dy ,
|
|
|
- viscosity_coeff_for_dv_dx ,
|
|
|
- viscosity_coeff_for_dv_dy,
|
|
|
- tmp_vec_p,
|
|
|
- tmp_vec_p0,
|
|
|
- rhs_phi,
|
|
|
- pres_free_surfaceL,jump_mass_transfer_rateL,mass_transfer_rate )
|
|
|
-
|
|
|
- PDI_status = @ccall "libpdi".PDI_multi_expose("velocity_extension"::Cstring,
|
|
|
- "u_ext_vel"::Cstring, u_ext_vel::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "v_ext_vel"::Cstring, v_ext_vel::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- C_NULL::Ptr{Cvoid})::Cint
|
|
|
-
|
|
|
- end
|
|
|
- #endregion ciprianoMulticomponentDropletEvaporation2024
|
|
|
-
|
|
|
-
|
|
|
- print("\n num.stop_simulation after NS ",num.stop_simulation)
|
|
|
-
|
|
|
-
|
|
|
- #region reactivate cut-cell for one-fluid model
|
|
|
- num.nLS = nb_levelsets
|
|
|
-
|
|
|
- grid_p.LS[1].u .= levelset_one_fluid
|
|
|
- grid_p.LS[end].u .= levelset_one_fluid
|
|
|
-
|
|
|
- NB_indices = update_all_ls_data(num, grid_p, grid_u, grid_v, BC_int, periodic_x, periodic_y)
|
|
|
-
|
|
|
- # TODO check after activation NS after nucleation
|
|
|
-
|
|
|
- # geoL = [grid_p.LS[iLS].geoL for iLS in 1:num._nLS]
|
|
|
- # geo_uL = [grid_u.LS[iLS].geoL for iLS in 1:num._nLS]
|
|
|
- # geo_vL = [grid_v.LS[iLS].geoL for iLS in 1:num._nLS]
|
|
|
-
|
|
|
- # reactivate centroids
|
|
|
- # printstyled(color=:red, @sprintf "\n reactivate centroids \n")
|
|
|
-
|
|
|
- # x_centroid = grid_p.x .+ getproperty.(grid_p.LS[1].geoL.centroid, :x) .* grid_p.dx #geoS
|
|
|
- # display(x_centroid)
|
|
|
-
|
|
|
- # y_centroid = grid_p.y .+ getproperty.(grid_p.LS[1].geoL.centroid, :y) .* grid_p.dy #geoS
|
|
|
- # display(y_centroid)
|
|
|
-
|
|
|
- # display(levelset_one_fluid)
|
|
|
- #endregion reactivate cut-cell for one-fluid model
|
|
|
-
|
|
|
-
|
|
|
- else
|
|
|
-
|
|
|
- if num.pressure_velocity_coupling == 0
|
|
|
- if ns_solid_phase
|
|
|
- geoS = [grid_p.LS[iLS].geoS for iLS in 1:num._nLS]
|
|
|
- geo_uS = [grid_u.LS[iLS].geoS for iLS in 1:num._nLS]
|
|
|
- geo_vS = [grid_v.LS[iLS].geoS for iLS in 1:num._nLS]
|
|
|
- Lpm1_S, bc_Lpm1_S, bc_Lpm1_b_S, Lum1_S, bc_Lum1_S, bc_Lum1_b_S, Lvm1_S, bc_Lvm1_S, bc_Lvm1_b_S,Mm1_S, Mum1_S, Mvm1_S, Cum1S, Cvm1S = pressure_projection!(
|
|
|
- time_scheme, BC_int,
|
|
|
- num, grid_p, geoS, grid_u, geo_uS, grid_v, geo_vS, phS,
|
|
|
- BC_uS, BC_vS, BC_pS,
|
|
|
- op.opC_pS, op.opC_uS, op.opC_vS, op.opS,
|
|
|
- AuS, BuS, AvS, BvS, AϕS, AuvS, BuvS,
|
|
|
- Lpm1_S, bc_Lpm1_S, bc_Lpm1_b_S, Lum1_S, bc_Lum1_S, bc_Lum1_b_S, Lvm1_S, bc_Lvm1_S, bc_Lvm1_b_S,
|
|
|
- Cum1S, Cvm1S, Mum1_S, Mvm1_S,
|
|
|
- periodic_x, periodic_y, ns_advection, advection, num.current_iter, Ra, navier,pres_free_surfaceS,jump_mass_transfer_rateS,mass_transfer_rateS
|
|
|
- )
|
|
|
- end
|
|
|
- if ns_liquid_phase
|
|
|
- geoL = [grid_p.LS[iLS].geoL for iLS in 1:num._nLS]
|
|
|
- geo_uL = [grid_u.LS[iLS].geoL for iLS in 1:num._nLS]
|
|
|
- geo_vL = [grid_v.LS[iLS].geoL for iLS in 1:num._nLS]
|
|
|
-
|
|
|
- # Mum1_L is put in B matrix that multiplies v
|
|
|
-
|
|
|
- Lpm1_L, bc_Lpm1_L, bc_Lpm1_b_L, Lum1_L, bc_Lum1_L, bc_Lum1_b_L, Lvm1_L, bc_Lvm1_L, bc_Lvm1_b_L, Mm1_L, Mum1_L, Mvm1_L, Cum1L, Cvm1L = pressure_projection!(
|
|
|
- time_scheme, BC_int,
|
|
|
- num, grid_p, geoL, grid_u, geo_uL, grid_v, geo_vL, phL,
|
|
|
- BC_uL, BC_vL, BC_pL,
|
|
|
- op.opC_pL, op.opC_uL, op.opC_vL, op.opL,
|
|
|
- AuL, BuL, AvL, BvL, AϕL, AuvL, BuvL,
|
|
|
- Lpm1_L, bc_Lpm1_L, bc_Lpm1_b_L, Lum1_L, bc_Lum1_L, bc_Lum1_b_L, Lvm1_L, bc_Lvm1_L, bc_Lvm1_b_L,
|
|
|
- Cum1L, Cvm1L, Mum1_L, Mvm1_L,
|
|
|
- periodic_x, periodic_y, ns_advection, advection, num.current_iter, Ra, navier,pres_free_surfaceL,jump_mass_transfer_rateL,mass_transfer_rate
|
|
|
- )
|
|
|
- # if num.current_iter == 1
|
|
|
- # phL.u .= -0.5 .* grid_u.y .+ getproperty.(grid_u.LS[1].geoL.centroid, :y) .* grid_u.dy
|
|
|
- # phL.v .= 0.5 .* grid_v.x .+ getproperty.(grid_v.LS[1].geoL.centroid, :x) .* grid_v.dx
|
|
|
- # phL.u[grid_u.LS[1].SOLID] .= 0.0
|
|
|
- # phL.v[grid_v.LS[1].SOLID] .= 0.0
|
|
|
- # end
|
|
|
- # linear_advection!(
|
|
|
- # num, grid_p, grid_p.LS[1].geoL, grid_u, grid_u.LS[1].geoL, grid_v, grid_v.LS[1].geoL, phL,
|
|
|
- # BC_uL, BC_vL, op.opL
|
|
|
- # )
|
|
|
- end
|
|
|
-
|
|
|
- elseif num.pressure_velocity_coupling > 1
|
|
|
-
|
|
|
- if ns_liquid_phase
|
|
|
- geoL = [grid_p.LS[iLS].geoL for iLS in 1:num._nLS]
|
|
|
- geo_uL = [grid_u.LS[iLS].geoL for iLS in 1:num._nLS]
|
|
|
- geo_vL = [grid_v.LS[iLS].geoL for iLS in 1:num._nLS]
|
|
|
-
|
|
|
- # Mum1_L is put in B matrix that multiplies v
|
|
|
-
|
|
|
- Lpm1_L, bc_Lpm1_L, bc_Lpm1_b_L, Lum1_L, bc_Lum1_L, bc_Lum1_b_L, Lvm1_L, bc_Lvm1_L, bc_Lvm1_b_L, Mm1_L, Mum1_L, Mvm1_L, Cum1L, Cvm1L = coupled_pressure_velocity!(
|
|
|
- time_scheme, BC_int,
|
|
|
- num, grid_p, geoL, grid_u, geo_uL, grid_v, geo_vL, phL,
|
|
|
- BC_uL, BC_vL, BC_pL,
|
|
|
- op.opC_pL, op.opC_uL, op.opC_vL, op.opL,
|
|
|
- AuL, BuL, AvL, BvL, AϕL, AuvL, BuvL,rhs_uv,
|
|
|
- Lpm1_L, bc_Lpm1_L, bc_Lpm1_b_L, Lum1_L, bc_Lum1_L, bc_Lum1_b_L, Lvm1_L, bc_Lvm1_L, bc_Lvm1_b_L,
|
|
|
- Cum1L, Cvm1L, Mum1_L, Mvm1_L,
|
|
|
- periodic_x, periodic_y, ns_advection, advection, num.current_iter, Ra, navier,pres_free_surfaceL,jump_mass_transfer_rateL,mass_transfer_rate
|
|
|
- )
|
|
|
- # if num.current_iter == 1
|
|
|
- # phL.u .= -0.5 .* grid_u.y .+ getproperty.(grid_u.LS[1].geoL.centroid, :y) .* grid_u.dy
|
|
|
- # phL.v .= 0.5 .* grid_v.x .+ getproperty.(grid_v.LS[1].geoL.centroid, :x) .* grid_v.dx
|
|
|
- # phL.u[grid_u.LS[1].SOLID] .= 0.0
|
|
|
- # phL.v[grid_v.LS[1].SOLID] .= 0.0
|
|
|
- # end
|
|
|
- # linear_advection!(
|
|
|
- # num, grid_p, grid_p.LS[1].geoL, grid_u, grid_u.LS[1].geoL, grid_v, grid_v.LS[1].geoL, phL,
|
|
|
- # BC_uL, BC_vL, op.opL
|
|
|
- # )
|
|
|
- end
|
|
|
-
|
|
|
- end #if num.pressure_velocity_coupling
|
|
|
-
|
|
|
- end # if (num.one_fluid_model == 1 && num.solve_Navier_Stokes_liquid_phase == 1)
|
|
|
-
|
|
|
- end # if navier_stokes
|
|
|
-
|
|
|
-
|
|
|
- #region conservation, divergence checks
|
|
|
- # TODO check global mass conservation and divergence free
|
|
|
- not_divergence_free = true
|
|
|
- #need one-fluid capa
|
|
|
- conservation = compute_conservation_mass(num,phL, grid_p ,grid_u, grid_v, rho_one_fluid)
|
|
|
-
|
|
|
- # Compute divergence of velocity
|
|
|
- Duv = op.opC_pL.AxT * vec1(phL.uD,grid_u) .+ op.opC_pL.Gx_b * vecb(phL.uD,grid_u) .+
|
|
|
- op.opC_pL.AyT * vec1(phL.vD,grid_v) .+ op.opC_pL.Gy_b * vecb(phL.vD,grid_v)
|
|
|
- for iLS in 1:num.nLS
|
|
|
- if !is_navier(BC_int[iLS]) && !is_navier_cl(BC_int[iLS]) #otherwise normal velocity null if no blowing
|
|
|
- Duv .+= op.opC_pL.Gx[iLS] * veci(phL.uD,grid_u,iLS+1) .+
|
|
|
- op.opC_pL.Gy[iLS] * veci(phL.vD,grid_v,iLS+1)
|
|
|
- end
|
|
|
- end
|
|
|
- #TODO divergence when Navier ?
|
|
|
-
|
|
|
- PDI_status = @ccall "libpdi".PDI_multi_expose("print_conservation"::Cstring,
|
|
|
- "nstep"::Cstring, num.current_iter::Ref{Clonglong}, PDI_OUT::Cint,
|
|
|
- "conservation"::Cstring, conservation::Ref{Cdouble}, PDI_OUT::Cint,
|
|
|
- "velocity_divergence"::Cstring, Duv::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # "p_1D"::Cstring, phL.pD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- C_NULL::Ptr{Cvoid})::Cint
|
|
|
-
|
|
|
- if maximum(Duv)< num.epsilon_divergence
|
|
|
- not_divergence_free = false
|
|
|
- end
|
|
|
-
|
|
|
- if abs(conservation) > num.epsilon_conservation || not_divergence_free
|
|
|
-
|
|
|
- PDI_status = @ccall "libpdi".PDI_multi_expose("conservation_error"::Cstring,
|
|
|
- "u_1D"::Cstring, phL.uD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "v_1D"::Cstring, phL.vD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "p_1D"::Cstring, phL.pD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- C_NULL::Ptr{Cvoid})::Cint
|
|
|
- end
|
|
|
- #endregion conservation, divergence checks
|
|
|
-
|
|
|
- # PDI_status = @ccall "libpdi".PDI_multi_expose("check_pressure_velocity"::Cstring,
|
|
|
- # "u_1D"::Cstring, phL.uD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # "v_1D"::Cstring, phL.vD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # "p_1D"::Cstring, phL.pD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # C_NULL::Ptr{Cvoid})::Cint
|
|
|
-
|
|
|
- end
|
|
|
- #endregion Navier-Stokes
|
|
|
-
|
|
|
-
|
|
|
- #region old block phase change
|
|
|
- #cf old_phase_change_block.jl
|
|
|
- #endregion old block phase change
|
|
|
-
|
|
|
- # printstyled(color=:red, @sprintf "\n after phase change radius: %.2e \n" num.current_radius)
|
|
|
-
|
|
|
- if verbose && adaptative_t
|
|
|
- println("num.timestep_n = $num.timestep_n")
|
|
|
- end
|
|
|
-
|
|
|
- # printstyled(color=:red, @sprintf "\n advection")
|
|
|
- if num.verbosity>0
|
|
|
- print("\n num.advection_LS_mode ",num.advection_LS_mode," advection ",advection)
|
|
|
- end
|
|
|
-
|
|
|
-
|
|
|
- #region Advection
|
|
|
- if advection || electrolysis_advection
|
|
|
-
|
|
|
- if num.io_pdi>0
|
|
|
-
|
|
|
- PDI_status = @ccall "libpdi".PDI_multi_expose("write_before_LS_adv"::Cstring,
|
|
|
- "normal_velocity_intfc"::Cstring, grid_p.V::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- C_NULL::Ptr{Cvoid})::Cint
|
|
|
-
|
|
|
- end # if num.io_pdi>0
|
|
|
-
|
|
|
- printstyled(color=:red, @sprintf "\n select advection")
|
|
|
-
|
|
|
-
|
|
|
- select_advection!(num, grid_p, BC_int, BC_u, grid_u, grid_v, CFL_sc, periodic_x, periodic_y,
|
|
|
- θ_out, rhs_LS, utmp, electrolysis_phase_change_case, mass_transfer_rate, levelset_1D,
|
|
|
- volume_fraction,
|
|
|
- tmp_vec_p,tmp_vec_p0,
|
|
|
- tmp_vec_u,tmp_vec_v,tmp_vec_u0,tmp_vec_v0,tmp_vec_u1,tmp_vec_v1,
|
|
|
- op,
|
|
|
- phL, u_ext_vel, v_ext_vel)
|
|
|
-
|
|
|
-
|
|
|
-
|
|
|
- # printstyled(color=:red, @sprintf "\n after advection_LS_mode radius: %.2e \n" num.current_radius)
|
|
|
-
|
|
|
- #region reinitialize Levelset
|
|
|
-
|
|
|
- #TODO document
|
|
|
- if num.levelset_reinitialize == 0
|
|
|
-
|
|
|
-
|
|
|
- if analytical
|
|
|
- u[grid_p.ind.b_top[1]] .= sqrt.(grid_p.x[grid_p.ind.b_top[1]] .^ 2 + grid_p.y[grid_p.ind.b_top[1]] .^ 2) .- (num.R + speed*num.current_iter*num.timestep_n);
|
|
|
- u[grid_p.ind.b_bottom[1]] .= sqrt.(grid_p.x[grid_p.ind.b_bottom[1]] .^ 2 + grid_p.y[grid_p.ind.b_bottom[1]] .^ 2) .- (num.R + speed*num.current_iter*num.timestep_n);
|
|
|
- u[grid_p.ind.b_left[1]] .= sqrt.(grid_p.x[grid_p.ind.b_left[1]] .^ 2 + grid_p.y[grid_p.ind.b_left[1]] .^ 2) .- (num.R + speed*num.current_iter*num.timestep_n);
|
|
|
- u[grid_p.ind.b_right[1]] .= sqrt.(grid_p.x[grid_p.ind.b_right[1]] .^ 2 + grid_p.y[grid_p.ind.b_right[1]] .^ 2) .- (num.R + speed*num.current_iter*num.timestep_n);
|
|
|
- elseif num.nb_reinit > 0
|
|
|
- if auto_reinit == 1 && (num.current_iter-1)%num.reinit_every == 0
|
|
|
- for iLS in 1:num.nLS
|
|
|
- if !is_wall(BC_int[iLS])
|
|
|
- ls_rg, rl_rg_v = rg(num, grid_p, grid_p.LS[iLS].u, periodic_x, periodic_y, BC_int)
|
|
|
- println("$(ls_rg)")
|
|
|
- printstyled(color=:green, @sprintf "\n ls_rg : %.2e \n" ls_rg)
|
|
|
- if ls_rg >= num.δreinit || num.current_iter == 1
|
|
|
- print("(ls_rg >= num.δreinit || num.current_iter == 1): yes")
|
|
|
- # println("yes")
|
|
|
- RK2_reinit!(ls_scheme, grid_p, grid_p.ind, iLS, grid_p.LS[iLS].u, num.nb_reinit, periodic_x, periodic_y, BC_u, BC_int)
|
|
|
-
|
|
|
- ls_rg, rl_rg_v = rg(num, grid_p, grid_p.LS[iLS].u, periodic_x, periodic_y, BC_int)
|
|
|
- println("$(ls_rg) ")
|
|
|
- printstyled(color=:green, @sprintf "\n ls_rg : %.2e \n" ls_rg)
|
|
|
- end
|
|
|
- end
|
|
|
- end
|
|
|
- elseif (num.current_iter-1)%num.reinit_every == 0
|
|
|
- for iLS in 1:num.nLS
|
|
|
- if !is_wall(BC_int[iLS])
|
|
|
- RK2_reinit!(ls_scheme, grid_p, grid_p.ind, iLS, grid_p.LS[iLS].u, num.nb_reinit, periodic_x, periodic_y, BC_u, BC_int)
|
|
|
- end
|
|
|
- end
|
|
|
- # elseif num.nLS > 1
|
|
|
- # for iLS in 1:num.nLS
|
|
|
- # if !is_wall(BC_int[iLS])
|
|
|
- # RK2_reinit!(ls_scheme, grid_p, grid_p.ind, iLS, grid_p.LS[iLS].u, 2num.nb_reinit, periodic_x, periodic_y, BC_u, BC_int, true)
|
|
|
- # end
|
|
|
- # end
|
|
|
- end
|
|
|
- end
|
|
|
-
|
|
|
- # Numerical breakup
|
|
|
- if free_surface && breakup ==1
|
|
|
- count, id_break = breakup_n(grid_p.LS[1].u, grid_p.nx, grid_p.ny, grid_p.dx, grid_p.dy, periodic_x, periodic_y, NB_indices, 5e-2)
|
|
|
- println(count)
|
|
|
- if count > count_limit_breakup
|
|
|
- println("BREAK UP!!")
|
|
|
- breakup_f(grid_p, grid_p.LS[1].u, id_break)
|
|
|
- RK2_reinit!(ls_scheme, grid_p, grid_p.ind, 1, grid_p.LS[1].u, num.nb_reinit, periodic_x, periodic_y, BC_u, BC_int)
|
|
|
- end
|
|
|
- end
|
|
|
-
|
|
|
- elseif num.levelset_reinitialize == -1
|
|
|
- printstyled(color=:red, @sprintf "\n No reinit")
|
|
|
-
|
|
|
-
|
|
|
- end #if num.levelset_reinitialize
|
|
|
-
|
|
|
- end
|
|
|
-
|
|
|
- #endregion reinitialize Levelset
|
|
|
-
|
|
|
- # printstyled(color=:red, @sprintf "\n after reinit radius: %.2e \n" num.current_radius)
|
|
|
-
|
|
|
- if verbose
|
|
|
- if (num.current_iter-1)%show_every == 0
|
|
|
- printstyled(color=:green, @sprintf "\n Current iteration : %d (%d%%) | t = %.2e \n" (num.current_iter) 100*(num.current_iter)/num.max_iterations num.time)
|
|
|
-
|
|
|
- #TODO CFL
|
|
|
- # printstyled(color=:green, @sprintf "\n num.CFL : %.2e num.CFL : %.2e num.timestep_n : %.2e\n" num.CFL max(abs.(grid_p.V)..., abs.(phL.u)..., abs.(phL.v)..., abs.(phS.u)..., abs.(phS.v)...)*num.timestep_n/num.Δ num.timestep_n)
|
|
|
-
|
|
|
- if heat && length(grid_p.LS[end].MIXED) != 0
|
|
|
- print(@sprintf "V_mean = %.2e V_max = %.2e V_min = %.2e\n" mean(grid_p.V[grid_p.LS[1].MIXED]) findmax(grid_p.V[grid_p.LS[1].MIXED])[1] findmin(grid_p.V[grid_p.LS[1].MIXED])[1])
|
|
|
- print(@sprintf "κ_mean = %.2e κ_max = %.2e κ_min = %.2e\n" mean(grid_p.LS[1].κ[grid_p.LS[1].MIXED]) findmax(grid_p.LS[1].κ[grid_p.LS[1].MIXED])[1] findmin(grid_p.LS[1].κ[grid_p.LS[1].MIXED])[1])
|
|
|
- elseif advection && length(grid_p.LS[end].MIXED) != 0
|
|
|
- V_mean = mean([mean(grid_u.V[grid_p.LS[1].MIXED]), mean(grid_v.V[grid_p.LS[1].MIXED])])
|
|
|
- V_max = max(findmax(grid_u.V[grid_p.LS[1].MIXED])[1], findmax(grid_v.V[grid_p.LS[1].MIXED])[1])
|
|
|
- V_min = min(findmin(grid_u.V[grid_p.LS[1].MIXED])[1], findmin(grid_v.V[grid_p.LS[1].MIXED])[1])
|
|
|
- # print(@sprintf "Vol_ratio = %.3f%%\n" (volume(grid_p.LS[end].geoL) / V0L * 100))
|
|
|
- print(@sprintf "V_mean = %.2e V_max = %.2e V_min = %.2e\n" V_mean V_max V_min)
|
|
|
- print(@sprintf "κ_mean = %.2e κ_max = %.2e κ_min = %.2e\n" mean(grid_p.LS[1].κ[grid_p.LS[1].MIXED]) findmax(grid_p.LS[1].κ[grid_p.LS[1].MIXED])[1] findmin(grid_p.LS[1].κ[grid_p.LS[1].MIXED])[1])
|
|
|
- end
|
|
|
- if navier_stokes
|
|
|
- if ns_solid_phase
|
|
|
- normuS = norm(phS.u)
|
|
|
- normvS = norm(phS.v)
|
|
|
- normpS = norm(phS.p.*num.timestep_n)
|
|
|
- print("$(@sprintf("norm(uS) %.6e", normuS))\t$(@sprintf("norm(vS) %.6e", normvS))\t$(@sprintf("norm(pS) %.6e", normpS))\n")
|
|
|
- end
|
|
|
- if ns_liquid_phase
|
|
|
- # normuL = norm(phL.u)
|
|
|
- # normvL = norm(phL.v)
|
|
|
- # normpL = norm(phL.p.*num.timestep_n)
|
|
|
- # print("$(@sprintf("norm(uL) %.6e", normuL))\t$(@sprintf("norm(vL) %.6e", normvL))\t$(@sprintf("norm(pL) %.6e", normpL))\n")
|
|
|
- if electrolysis
|
|
|
- # print_electrolysis_statistics(num,grid_p,phL)
|
|
|
- PDI_status = @ccall "libpdi".PDI_multi_expose("print_variables"::Cstring,
|
|
|
- "nstep"::Cstring, num.current_iter ::Ref{Clonglong}, PDI_OUT::Cint,
|
|
|
- "time"::Cstring, num.time::Ref{Cdouble}, PDI_OUT::Cint,
|
|
|
- "u_1D"::Cstring, phL.uD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "v_1D"::Cstring, phL.vD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "p_1D"::Cstring, phL.pD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "levelset_p"::Cstring, grid_p.LS[num.iLSpdi].u::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "levelset_u"::Cstring, grid_u.LS[num.iLSpdi].u::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "levelset_v"::Cstring, grid_v.LS[num.iLSpdi].u::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "trans_scal_1DT"::Cstring, phL.trans_scalD'::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "phi_ele_1D"::Cstring, phL.phi_eleD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- C_NULL::Ptr{Cvoid})::Cint
|
|
|
- end
|
|
|
- end
|
|
|
- end
|
|
|
- end
|
|
|
- end
|
|
|
-
|
|
|
- PDI_status = @ccall "libpdi".PDI_multi_expose("write_after_advection"::Cstring,
|
|
|
- "nstep"::Cstring, num.current_iter::Ref{Clonglong}, PDI_OUT::Cint,
|
|
|
- "time"::Cstring, num.time::Ref{Cdouble}, PDI_OUT::Cint,
|
|
|
- "levelset_p"::Cstring, grid_p.LS[num.iLSpdi].u::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "levelset_u"::Cstring, grid_u.LS[num.iLSpdi].u::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "levelset_v"::Cstring, grid_v.LS[num.iLSpdi].u::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- C_NULL::Ptr{Cvoid})::Cint
|
|
|
-
|
|
|
- # if levelset && (advection || num.current_iter<2 || electrolysis_advection)
|
|
|
- if levelset && (advection || num.current_iter<2)
|
|
|
- try
|
|
|
- NB_indices = update_all_ls_data(num, grid_p, grid_u, grid_v, BC_int, periodic_x, periodic_y)
|
|
|
- catch errorLS
|
|
|
- println(@sprintf "\n CRASHED after %d iterations \n" num.current_iter)
|
|
|
- printstyled(color=:red, @sprintf "\n grid_p.LS not updated \n")
|
|
|
- print(errorLS)
|
|
|
- print(errorLS.task.exception)
|
|
|
- return
|
|
|
- end
|
|
|
- # printstyled(color=:red, @sprintf "\n levelset 4:\n")
|
|
|
- # println(grid_p.LS[1].geoL.dcap[1,1,:])
|
|
|
-
|
|
|
- grid_p.LS[end].geoL.fresh .= false
|
|
|
- grid_u.LS[end].geoL.fresh .= false
|
|
|
- grid_v.LS[end].geoL.fresh .= false
|
|
|
-
|
|
|
- get_fresh_cells!(grid_p, grid_p.LS[end].geoL, Mm1_L, grid_p.ind.all_indices)
|
|
|
- get_fresh_cells!(grid_u, grid_u.LS[end].geoL, Mum1_L, grid_u.ind.all_indices)
|
|
|
- get_fresh_cells!(grid_v, grid_v.LS[end].geoL, Mvm1_L, grid_v.ind.all_indices)
|
|
|
-
|
|
|
- FRESH_L_u = findall(grid_u.LS[end].geoL.fresh)
|
|
|
- FRESH_L_v = findall(grid_v.LS[end].geoL.fresh)
|
|
|
-
|
|
|
- if navier_stokes
|
|
|
-
|
|
|
- init_fresh_cells!(grid_u, veci(phL.uD,grid_u,1), veci(phL.uD,grid_u,1),
|
|
|
- grid_u.LS[end].geoL.projection, FRESH_L_u, periodic_x, periodic_y)
|
|
|
- init_fresh_cells!(grid_v, veci(phL.vD,grid_v,1), veci(phL.vD,grid_v,1),
|
|
|
- grid_v.LS[end].geoL.projection, FRESH_L_v, periodic_x, periodic_y)
|
|
|
- if num.nLS>1 #not monophasic
|
|
|
- init_fresh_cells!(grid_u, veci(phL.uD,grid_u,2), veci(phL.uD,grid_u,1),
|
|
|
- grid_u.LS[end].geoL.projection, FRESH_L_u, periodic_x, periodic_y)
|
|
|
- init_fresh_cells!(grid_v, veci(phL.vD,grid_v,2), veci(phL.vD,grid_v,1),
|
|
|
- grid_v.LS[end].geoL.projection, FRESH_L_v, periodic_x, periodic_y)
|
|
|
- end
|
|
|
- end
|
|
|
-
|
|
|
- if num.solve_solid == 1
|
|
|
- grid_p.LS[end].geoS.fresh .= false
|
|
|
- grid_u.LS[end].geoS.fresh .= false
|
|
|
- grid_v.LS[end].geoS.fresh .= false
|
|
|
-
|
|
|
- get_fresh_cells!(grid_p, grid_p.LS[end].geoS, Mm1_S, grid_p.ind.all_indices)
|
|
|
- get_fresh_cells!(grid_u, grid_u.LS[end].geoS, Mum1_S, grid_u.ind.all_indices)
|
|
|
- get_fresh_cells!(grid_v, grid_v.LS[end].geoS, Mvm1_S, grid_v.ind.all_indices)
|
|
|
-
|
|
|
- FRESH_S_u = findall(grid_u.LS[end].geoS.fresh)
|
|
|
- FRESH_S_v = findall(grid_v.LS[end].geoS.fresh)
|
|
|
-
|
|
|
- if navier_stokes
|
|
|
- init_fresh_cells!(grid_u, veci(phL.uD,grid_u,1), veci(phL.uD,grid_u,1),
|
|
|
- grid_u.LS[end].geoL.projection, FRESH_L_u, periodic_x, periodic_y)
|
|
|
- init_fresh_cells!(grid_v, veci(phL.vD,grid_v,1), veci(phL.vD,grid_v,1),
|
|
|
- grid_v.LS[end].geoL.projection, FRESH_L_v, periodic_x, periodic_y)
|
|
|
- if num.nLS>1 #not monophasic
|
|
|
- init_fresh_cells!(grid_u, veci(phL.uD,grid_u,2), veci(phL.uD,grid_u,1),
|
|
|
- grid_u.LS[end].geoL.projection, FRESH_L_u, periodic_x, periodic_y)
|
|
|
- init_fresh_cells!(grid_v, veci(phL.vD,grid_v,2), veci(phL.vD,grid_v,1),
|
|
|
- grid_v.LS[end].geoL.projection, FRESH_L_v, periodic_x, periodic_y)
|
|
|
- end
|
|
|
- end
|
|
|
- end
|
|
|
-
|
|
|
- if iszero(num.current_iter%num.save_every) || num.current_iter==num.max_iterations
|
|
|
- snap = num.current_iter÷num.save_every+1
|
|
|
- if save_radius
|
|
|
- radius[snap] = find_radius(grid_p, grid_p.LS[1])
|
|
|
- end
|
|
|
- if hill
|
|
|
- a = zeros(length(grid_p.LS[1].MIXED))
|
|
|
- for i in eachindex(grid_p.LS[1].MIXED)
|
|
|
- a[i] = grid_p.LS[1].geoL.projection[grid_p.LS[1].MIXED[i]].pos.y
|
|
|
- end
|
|
|
- radius[snap] = mean(a)
|
|
|
- end
|
|
|
- # if save_length
|
|
|
- # fwd.length[snap] = arc_length2(grid_p.LS[1].geoS.projection, grid_p.LS[1].MIXED)
|
|
|
- # end
|
|
|
- end
|
|
|
- end
|
|
|
-
|
|
|
-
|
|
|
- if num.nLS>1 #not monophasic
|
|
|
- intfc_vtx_x,intfc_vtx_y,intfc_vtx_field,intfc_vtx_connectivities,intfc_vtx_num, intfc_seg_num = convert_interfacial_D_to_segments(num,grid_p,phL.TD,1,2)
|
|
|
-
|
|
|
- barycenter_x_coord = mean(intfc_vtx_x)
|
|
|
-
|
|
|
- PDI_status = @ccall "libpdi".PDI_multi_expose("update_levelset"::Cstring,
|
|
|
- "nstep"::Cstring, num.current_iter ::Ref{Clonglong}, PDI_OUT::Cint,
|
|
|
- "time"::Cstring, num.time::Ref{Cdouble}, PDI_OUT::Cint,
|
|
|
- "levelset_p"::Cstring, grid_p.LS[num.iLSpdi].u::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "levelset_u"::Cstring, grid_u.LS[num.iLSpdi].u::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "levelset_v"::Cstring, grid_v.LS[num.iLSpdi].u::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "intfc_vtx_num"::Cstring, intfc_vtx_num::Ref{Clonglong}, PDI_OUT::Cint,
|
|
|
- "intfc_seg_num"::Cstring, intfc_seg_num::Ref{Clonglong}, PDI_OUT::Cint,
|
|
|
- "intfc_vtx_x"::Cstring, intfc_vtx_x::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "intfc_vtx_y"::Cstring, intfc_vtx_y::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "intfc_vtx_field"::Cstring, intfc_vtx_field::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "intfc_vtx_connectivities"::Cstring, intfc_vtx_connectivities::Ptr{Clonglong}, PDI_OUT::Cint,
|
|
|
- "barycenter_x_coord"::Cstring, barycenter_x_coord::Ref{Cdouble}, PDI_OUT::Cint,
|
|
|
- C_NULL::Ptr{Cvoid})::Cint
|
|
|
- end
|
|
|
-
|
|
|
- #endregion Advection
|
|
|
-
|
|
|
- #region old Navier-Stokes block
|
|
|
- # cf old_Navier_Stokes_block.jl
|
|
|
- #endregion old Navier-Stokes block
|
|
|
-
|
|
|
-
|
|
|
- #region end loop post-processing
|
|
|
-
|
|
|
- if (num.one_fluid_model == 1 && num.solve_Navier_Stokes_liquid_phase == 1)
|
|
|
-
|
|
|
- # interpolate_staggered_u_v_to_scalar_grid_one_fluid_or_one_phase!(num,grid_p,grid_u,grid_v,phL.u,phL.v,tmp_vec_p,tmp_vec_p0)
|
|
|
-
|
|
|
- tmp_vec_p .= 0.0
|
|
|
- tmp_vec_p0 .= 0.0
|
|
|
-
|
|
|
- for j = 1:grid_p.ny
|
|
|
- for i = 1:grid_p.nx
|
|
|
- tmp_vec_p[j,i] =(phL.u[j,i]+phL.u[j,i+1])/2
|
|
|
- tmp_vec_p0[j,i]=(phL.v[j,i]+phL.v[j+1,i])/2
|
|
|
- end
|
|
|
- end
|
|
|
-
|
|
|
- #TODO compute once, write once
|
|
|
- update_one_fluid_density_viscosity(num,grid_p,grid_u,grid_v,volume_fraction,levelset_one_fluid,rho_one_fluid,
|
|
|
- rho_one_fluid_u,rho_one_fluid_v,mu_one_fluid,tmp_vec_p0)
|
|
|
-
|
|
|
- end
|
|
|
-
|
|
|
- #endregion end loop post-processing
|
|
|
-
|
|
|
- # cD, cL, D, L = force_coefficients!(num, grid_p, grid_u, grid_v, op.opL, fwd, phL; step = num.current_iter+1, saveCoeffs = false)
|
|
|
-
|
|
|
- # if iszero(num.current_iter%num.save_every) || num.current_iter==num.max_iterations
|
|
|
- # snap = num.current_iter÷num.save_every+1
|
|
|
- # if num.current_iter==num.max_iterations
|
|
|
- # snap = size(fwd.T,1)
|
|
|
- # end
|
|
|
- # fwd.t[snap] = num.time
|
|
|
- # @views fwd.V[snap,:,:] .= grid_p.V
|
|
|
- # if advection
|
|
|
- # fwdS.Vratio[snap] = volume(grid_p.LS[end].geoS) / V0S
|
|
|
- # fwdL.Vratio[snap] = volume(grid_p.LS[end].geoL) / V0L
|
|
|
- # end
|
|
|
- # end
|
|
|
- # @views fwd.Cd[num.current_iter+1] = cD
|
|
|
- # @views fwd.Cl[num.current_iter+1] = cL
|
|
|
- # # @views fwd.radius[num.current_iter+1] = num.current_radius
|
|
|
-
|
|
|
- # PDI (IO)
|
|
|
- if electrolysis
|
|
|
- if num.io_pdi>0
|
|
|
-
|
|
|
- try
|
|
|
- # printstyled(color=:red, @sprintf "\n PDI test \n" )
|
|
|
-
|
|
|
-
|
|
|
- # phi_array=phL.phi_ele #do not transpose since python row major
|
|
|
-
|
|
|
- # Compute electrical current, interpolate velocity on scalar grid_p
|
|
|
- # if num.electrical_potential>0 compute_grad_phi_ele!(num, grid_p, grid_u, grid_v, phL, phS, op.opC_pL, op.opC_pS) #TODO current
|
|
|
-
|
|
|
- if num.electrical_potential>0
|
|
|
- compute_grad_phi_ele!(num, grid_p, grid_u, grid_v, grid_u.LS[end], grid_v.LS[end], phL,
|
|
|
- op.opC_pL, elec_cond,tmp_vec_u,tmp_vec_v,tmp_vec_p,tmp_vec_p0,tmp_vec_p1) #TODO current
|
|
|
- end
|
|
|
-
|
|
|
- # #store in us, vs instead of Eus, Evs
|
|
|
- # interpolate_grid_liquid!(grid_p,grid_u,grid_v,phL.Eu, phL.Ev,tmp_vec_p,tmp_vec_p0)
|
|
|
-
|
|
|
- # #TODO i_current_mag need cond
|
|
|
-
|
|
|
- # @ccall "libpdi".PDI_multi_expose("write_data_elec"::Cstring,
|
|
|
- # "i_current_x"::Cstring, tmp_vec_p::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # "i_current_y"::Cstring, tmp_vec_p0::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # "i_current_mag"::Cstring, phL.i_current_mag::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # "phi_ele_1D"::Cstring, phL.phi_eleD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- # C_NULL::Ptr{Cvoid})::Cint
|
|
|
-
|
|
|
- # interpolate_grid_liquid!(grid_p,grid_u,grid_v,phL.Eu, phL.Ev,Eus,Evs)
|
|
|
-
|
|
|
- interpolate_staggered_u_v_to_scalar_grid_one_fluid_or_one_phase!(num,grid_p,grid_u,grid_v,phL.u,phL.v,tmp_vec_p,tmp_vec_p0)
|
|
|
-
|
|
|
- # num.iLSpdi = 1 # TODO all grid_p.LS
|
|
|
-
|
|
|
- # Exposing data to PDI for IO
|
|
|
- # if writing "D" array (bulk, interface, border), add "_1D" to the name
|
|
|
-
|
|
|
- PDI_status = @ccall "libpdi".PDI_multi_expose("write_data"::Cstring,
|
|
|
- "nstep"::Cstring, num.current_iter::Ref{Clonglong}, PDI_OUT::Cint,
|
|
|
- "time"::Cstring, num.time::Ref{Cdouble}, PDI_OUT::Cint,
|
|
|
- "timestep"::Cstring, num.timestep_n::Ref{Cdouble}, PDI_OUT::Cint,
|
|
|
- "nx"::Cstring, grid_p.nx::Ref{Clonglong}, PDI_OUT::Cint,
|
|
|
- "ny"::Cstring, grid_p.ny::Ref{Clonglong}, PDI_OUT::Cint,
|
|
|
- "u_1D"::Cstring, phL.uD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "v_1D"::Cstring, phL.vD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "p_1D"::Cstring, phL.pD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "levelset_p"::Cstring, grid_p.LS[num.iLSpdi].u::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "levelset_u"::Cstring, grid_u.LS[num.iLSpdi].u::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "levelset_v"::Cstring, grid_v.LS[num.iLSpdi].u::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "trans_scal_1DT"::Cstring, phL.trans_scalD'::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "phi_ele_1D"::Cstring, phL.phi_eleD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "velocity_x"::Cstring, tmp_vec_p::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "velocity_y"::Cstring, tmp_vec_p0::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "radius"::Cstring, num.current_radius::Ref{Cdouble}, PDI_OUT::Cint,
|
|
|
- C_NULL::Ptr{Cvoid})::Cint
|
|
|
-
|
|
|
-
|
|
|
-
|
|
|
- catch error
|
|
|
- printstyled(color=:red, @sprintf "\n PDI error \n")
|
|
|
- print(error)
|
|
|
- printstyled(color=:red, @sprintf "\n PDI error \n")
|
|
|
- end
|
|
|
-
|
|
|
-
|
|
|
- end #if io_pdi
|
|
|
-
|
|
|
-
|
|
|
-
|
|
|
- if num.status == 1 #due to num.nH2<0...
|
|
|
- return
|
|
|
- end
|
|
|
-
|
|
|
-
|
|
|
- #region test NaN
|
|
|
- if num.solve_solid == 1
|
|
|
-
|
|
|
- if (any(isnan, phL.uD) || any(isnan, phL.vD) || any(isnan, phL.TD) ||
|
|
|
- any(isnan, phS.uD) || any(isnan, phS.vD) || any(isnan, phS.TD) ||
|
|
|
- any(isnan, phL.trans_scalD) || any(isnan, phL.phi_eleD) ||
|
|
|
- norm(phL.u) > 1e8 || norm(phL.T) > 1e8 ||
|
|
|
- norm(phS.u) > 1e8 || norm(phS.T) > 1e8 ||
|
|
|
- # norm(phL.trans_scal) > 1e8 ||
|
|
|
- norm(phL.phi_ele) > 1e8
|
|
|
- # ||
|
|
|
- # any(phL.trans_scal .<0)
|
|
|
- )
|
|
|
- println(@sprintf "\n CRASHED start \n")
|
|
|
-
|
|
|
- print("\n phL.uD: ",any(isnan, phL.uD) , "\n phL.vD: ",any(isnan, phL.vD) , "\n phL.TD: ",any(isnan, phL.TD) ,
|
|
|
- "\n phS.uD: ",any(isnan, phS.uD) , "\n phS.vD: ",any(isnan, phS.vD) , "\n phS.TD: ",any(isnan, phS.TD) ,
|
|
|
- "\n phL.trans_scalD: ",any(isnan, phL.trans_scalD) , "\n phL.phi_eleD: ",any(isnan, phL.phi_eleD) ,
|
|
|
- "\n phL.u: ",norm(phL.u) > 1e8 , "\n phS.u: ",norm(phS.u) > 1e8 , "\n phL.T: ",norm(phL.T) > 1e8 ,
|
|
|
- "\n phS.T: ",norm(phS.T) > 1e8 , "\n phL.trans_scal: ",norm(phL.trans_scal) > 1e8 ,
|
|
|
- "\n phL.phi_ele: ",norm(phL.phi_ele) > 1e8,"\n any(phL.trans_scal .<0): ", any(phL.trans_scal .<0))
|
|
|
-
|
|
|
- num.status = 1
|
|
|
-
|
|
|
- end
|
|
|
- else
|
|
|
-
|
|
|
- if (any(isnan, phL.uD) || any(isnan, phL.vD) || any(isnan, phL.TD) ||
|
|
|
- any(isnan, phL.trans_scalD) || any(isnan, phL.phi_eleD) ||
|
|
|
- norm(phL.u) > 1e8 || norm(phL.T) > 1e8 ||
|
|
|
- # norm(phL.trans_scal) > 1e8 ||
|
|
|
- norm(phL.phi_ele) > 1e8
|
|
|
- # ||
|
|
|
- # any(phL.trans_scal .<0)
|
|
|
- )
|
|
|
- println(@sprintf "\n CRASHED start \n")
|
|
|
-
|
|
|
- # println(@sprintf "\n CRASHED after %d iterations \n" num.current_iter)
|
|
|
-
|
|
|
- print("\n phL.uD: ",any(isnan, phL.uD) , "\n phL.vD: ",any(isnan, phL.vD) , "\n phL.TD: ",any(isnan, phL.TD) ,
|
|
|
- "\n phL.trans_scalD: ",any(isnan, phL.trans_scalD) , "\n phL.phi_eleD: ",any(isnan, phL.phi_eleD) ,
|
|
|
- "\n phL.u: ",norm(phL.u) > 1e8, "\n phL.T: ",norm(phL.T) > 1e8 ,
|
|
|
- "\n phL.trans_scal: ",norm(phL.trans_scal) > 1e8 ,
|
|
|
- "\n phL.phi_ele: ",norm(phL.phi_ele) > 1e8,"\n any(phL.trans_scal .<0): ", any(phL.trans_scal .<0))
|
|
|
-
|
|
|
- num.status = 1
|
|
|
-
|
|
|
- end
|
|
|
-
|
|
|
- end
|
|
|
-
|
|
|
- if num.status == 1
|
|
|
-
|
|
|
- PDI_status = @ccall "libpdi".PDI_multi_expose("print_variables"::Cstring,
|
|
|
- "nstep"::Cstring, num.current_iter ::Ref{Clonglong}, PDI_OUT::Cint,
|
|
|
- "time"::Cstring, num.time::Ref{Cdouble}, PDI_OUT::Cint,
|
|
|
- "u_1D"::Cstring, phL.uD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "v_1D"::Cstring, phL.vD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "p_1D"::Cstring, phL.pD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "levelset_p"::Cstring, grid_p.LS[num.iLSpdi].u::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "levelset_u"::Cstring, grid_u.LS[num.iLSpdi].u::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "levelset_v"::Cstring, grid_v.LS[num.iLSpdi].u::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "trans_scal_1DT"::Cstring, phL.trans_scalD'::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "phi_ele_1D"::Cstring, phL.phi_eleD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- C_NULL::Ptr{Cvoid})::Cint
|
|
|
-
|
|
|
- return num.current_iter
|
|
|
- end
|
|
|
-
|
|
|
-
|
|
|
-
|
|
|
- else #not electrolysis
|
|
|
-
|
|
|
- if num.solve_solid == 1
|
|
|
-
|
|
|
- if (any(isnan, phL.uD) || any(isnan, phL.vD) || any(isnan, phL.TD) ||
|
|
|
- any(isnan, phS.uD) || any(isnan, phS.vD) || any(isnan, phS.TD) ||
|
|
|
- norm(phL.u) > 1e8 || norm(phS.u) > 1e8 || norm(phL.T) > 1e8 || norm(phS.T) > 1e8)
|
|
|
- println(@sprintf "\n CRASHED after %d iterations \n" num.current_iter)
|
|
|
-
|
|
|
- num.status = 1
|
|
|
- return
|
|
|
-
|
|
|
- end
|
|
|
-
|
|
|
- else
|
|
|
-
|
|
|
- if (any(isnan, phL.uD) || any(isnan, phL.vD) || any(isnan, phL.TD) ||
|
|
|
- norm(phL.u) > 1e8 || norm(phL.T) > 1e8 )
|
|
|
- println(@sprintf "\n CRASHED after %d iterations \n" num.current_iter)
|
|
|
-
|
|
|
- num.status = 1
|
|
|
- return
|
|
|
-
|
|
|
- end
|
|
|
-
|
|
|
- end
|
|
|
-
|
|
|
- end #if electrolysis
|
|
|
-
|
|
|
- #endregion test NaN
|
|
|
-
|
|
|
-
|
|
|
- #region update iter number and time
|
|
|
-
|
|
|
- # num.current_iter += 1
|
|
|
-
|
|
|
- if num.time + num.timestep_n > num.end_time
|
|
|
- print("\n num.time + num.timestep_n > num.end_time, break")
|
|
|
- break
|
|
|
- end
|
|
|
-
|
|
|
-
|
|
|
- #endregion update iter number and time
|
|
|
-
|
|
|
- end
|
|
|
- #endregion time loop
|
|
|
-
|
|
|
- #region print end
|
|
|
- if verbose
|
|
|
- try
|
|
|
- printstyled(color=:blue, @sprintf "\n Final iteration : %d (%d%%) | t = %.2e \n" (num.current_iter-1) 100*(num.current_iter-1)/num.max_iterations num.time)
|
|
|
- print("\n num.time ",num.time," stop_simulation ",num.stop_simulation)
|
|
|
- print("\n num.time ",num.time," num.end_time",num.end_time,"num.timestep_n",num.timestep_n)
|
|
|
- print("\n")
|
|
|
- if stefan && advection
|
|
|
- print(@sprintf "V_mean = %.2e V_max = %.2e V_min = %.2e V_stdev = %.5f\n" mean(grid_p.V[grid_p.LS[1].MIXED]) findmax(grid_p.V[grid_p.LS[1].MIXED])[1] findmin(grid_p.V[grid_p.LS[1].MIXED])[1] std(grid_p.V[grid_p.LS[1].MIXED]))
|
|
|
- print(@sprintf "κ_mean = %.2e κ_max = %.2e κ_min = %.2e κ_stdev = %.5f\n" mean(grid_p.LS[1].κ[grid_p.LS[1].MIXED]) findmax(grid_p.LS[1].κ[grid_p.LS[1].MIXED])[1] findmin(grid_p.LS[1].κ[grid_p.LS[1].MIXED])[1] std(grid_p.LS[1].κ[grid_p.LS[1].MIXED]))
|
|
|
- end
|
|
|
- if free_surface && advection
|
|
|
- # print(@sprintf "Vol_ratio = %.3f%%\n" (volume(grid_p.LS[end].geoL) / V0L * 100))
|
|
|
- print(@sprintf "V_mean = %.2e V_max = %.2e V_min = %.2e V_stdev = %.5f\n" mean(grid_p.V[grid_p.LS[1].MIXED]) findmax(grid_p.V[grid_p.LS[1].MIXED])[1] findmin(grid_p.V[grid_p.LS[1].MIXED])[1] std(grid_p.V[grid_p.LS[1].MIXED]))
|
|
|
- print(@sprintf "κ_mean = %.2e κ_max = %.2e κ_min = %.2e κ_stdev = %.5f\n" mean(grid_p.LS[1].κ[grid_p.LS[1].MIXED]) findmax(grid_p.LS[1].κ[grid_p.LS[1].MIXED])[1] findmin(grid_p.LS[1].κ[grid_p.LS[1].MIXED])[1] std(grid_p.LS[1].κ[grid_p.LS[1].MIXED]))
|
|
|
- end
|
|
|
- if navier_stokes
|
|
|
- if ns_solid_phase
|
|
|
- normuS = norm(phS.u)
|
|
|
- normvS = norm(phS.v)
|
|
|
- normpS = norm(phS.p.*num.timestep_n)
|
|
|
- print("$(@sprintf("norm(uS) %.6e", normuS))\t$(@sprintf("norm(vS) %.6e", normvS))\t$(@sprintf("norm(pS) %.6e", normpS))\n")
|
|
|
- end
|
|
|
- if ns_liquid_phase
|
|
|
- normuL = norm(phL.u)
|
|
|
- normvL = norm(phL.v)
|
|
|
- normpL = norm(phL.p.*num.timestep_n)
|
|
|
- print("$(@sprintf("norm(uL) %.6e", normuL))\t$(@sprintf("norm(vL) %.6e", normvL))\t$(@sprintf("norm(pL) %.6e", normpL))\n")
|
|
|
- if electrolysis
|
|
|
- print_electrolysis_statistics(num,grid_p,phL)
|
|
|
- PDI_status = @ccall "libpdi".PDI_multi_expose("print_variables"::Cstring,
|
|
|
- "nstep"::Cstring, num.current_iter ::Ref{Clonglong}, PDI_OUT::Cint,
|
|
|
- "time"::Cstring, num.time::Ref{Cdouble}, PDI_OUT::Cint,
|
|
|
- "u_1D"::Cstring, phL.uD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "v_1D"::Cstring, phL.vD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "p_1D"::Cstring, phL.pD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "levelset_p"::Cstring, grid_p.LS[num.iLSpdi].u::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "levelset_u"::Cstring, grid_u.LS[num.iLSpdi].u::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "levelset_v"::Cstring, grid_v.LS[num.iLSpdi].u::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "trans_scal_1DT"::Cstring, phL.trans_scalD'::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- "phi_ele_1D"::Cstring, phL.phi_eleD::Ptr{Cdouble}, PDI_OUT::Cint,
|
|
|
- C_NULL::Ptr{Cvoid})::Cint
|
|
|
-
|
|
|
- end
|
|
|
- end
|
|
|
- end
|
|
|
- print("\n\n")
|
|
|
- catch
|
|
|
- @show (length(grid_p.LS[end].MIXED))
|
|
|
- end
|
|
|
- end #if verbose
|
|
|
- #endregion print end
|
|
|
-
|
|
|
- if levelset && (save_radius || hill)
|
|
|
- #TODO save radius
|
|
|
- return # radius
|
|
|
- # elseif flapping
|
|
|
- # return xc, yc
|
|
|
- else
|
|
|
- return
|
|
|
- end
|
|
|
-end
|
|
|
-
|
