""" 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