|
|
@@ -0,0 +1,3292 @@
|
|
|
+"""
|
|
|
+
|
|
|
+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
|
|
|
+
|
