# using Revise using Flower # PDI localARGS = ARGS # @show localARGS # print("\n Arguments ", localARGS) # print("\n length(localARGS) ",length(localARGS)) if length(localARGS)>0 yamlfile = localARGS[1] printstyled(color=:magenta, @sprintf "\n YAML file ") print(yamlfile) print("\n") yamlpath = yamlfile yamlnamelist = split(yamlfile, "/") yamlname = yamlnamelist[end] end data = YAML.load_file(yamlpath) # Dictionaries prop_dict = PropertyDict(data) #region study parameters study = PropertyDict(prop_dict.study) #mesh connvergence nb_grid_points = study.meshes n_cases = length(nb_grid_points) print("\n number of points ", nb_grid_points, "\n") #timestep convergence timesteps = study.timesteps #endregion study parameters io = PropertyDict(prop_dict.plot) flower = PropertyDict(prop_dict.flower) mesh = PropertyDict(flower.mesh) sim = PropertyDict(flower.simulation) phys = PropertyDict(flower.physics) macros = PropertyDict(flower.macros) #to parse code from .yml # boundaries_dict = PropertyDict(macros.boundaries_list) study_name = "" #region change one parameter at a time # # print("\n changing parameters",Meta.parseall(study.macro)) # change_one_parameter_at_a_time = PropertyDict(study.change_one_parameter_at_a_time) # # eval(Meta.parseall(study.change_one_parameter_at_a_time.macro)) # # for # # study_name = change_one_parameter_at_a_time # for (key, value) in change_one_parameter_at_a_time # print("\n changing parameter") # print("\n key ",key) # print("\n value ",value) # # eval(value.macro) # print("\n mu ",phys.mu1,phys.mu2) # eval(value["macro"]) # print("\n mu ",phys.mu1,phys.mu2) # end # # Step 2: Modify multiple parameters # # Define a dictionary with the parameters you want to change and their new values # # changes = Dict( # # "parameter1" => "new_value1", # # "parameter2" => "new_value2", # # "parameter3" => "new_value3" # # # Add more parameters as needed # # ) # # changes = eval(study.change_one_parameter_at_a_time) # # Apply the changes to the original YAML data # # for (key, value) in study.change_one_parameter_at_a_time # # print("\n changing parameter") # # print(key) # # eval(value.macro) # # end #endregion change one parameter at a time # print parameters by evaluating Julia code stored in .yml eval(Meta.parseall(macros.print_parameters)) #region attempt at precompiling Flower modules to librairies (.so) # print("\n test juliac") # # @ccall "simple.so".add_julia(2::Cint, 2::Cint)::Cint # # test_juliac = @ccall "simple.so".add_julia(2::Cint, 2::Cint)::Cint # # # Get the current working directory # # current_directory = pwd() # # # Print the current working directory # # println("Current working directory: ", current_directory) # # test_juliac = @ccall "./simple.so".add_julia(2::Cint, 2::Cint)::Cint # # test_juliac = @ccall "/local/home/pr277828/flower/juliac/simple.so".add_julia(2::Cint, 2::Cint)::Cint # test_juliac = @ccall add_julia(2::Cint, 2::Cint)::Cint # print("\n test ",test_juliac) # # juliac_status = @ccall "mylib".load_parameters()::Cint # # juliac_status = @ccall "libmylib.so".load_parameters()::Cint # print("\n test juliac") # # pdipath = "libpdi" # pdipath = "/usr/lib/x86_64-linux-gnu/libpdi.so" #endregion attempt at precompiling Flower modules to librairies (.so) if io.pdi>0 @debug "Before PDI init" yml_file = yamlfile conf = @ccall "libparaconf".PC_parse_path(yml_file::Cstring)::PC_tree_t @debug "after conf" getsubyml = @ccall "libparaconf".PC_get(conf::PC_tree_t,".pdi"::Cstring)::PC_tree_t @debug "after getsubyml" local pdi_status = @ccall "libpdi".PDI_init(getsubyml::PC_tree_t)::Cint @debug "after PDI_init" # Send meta-data to PDI mpi_coords_x = 1 mpi_coords_y = 1 mpi_max_coords_x = 1 mpi_max_coords_y = 1 local nx = 32 local ny = 32 nstep = 0 # nx=gp.nx # ny=gp.ny #TODO check Clonglong ... phys_time = 0.0 #Cdouble # nstep = num.current_iter local PDI_status = @ccall "libpdi".PDI_multi_expose("init_PDI"::Cstring, "mpi_coords_x"::Cstring, mpi_coords_x::Ref{Clonglong}, PDI_OUT::Cint, "mpi_coords_y"::Cstring, mpi_coords_x::Ref{Clonglong}, PDI_OUT::Cint, "mpi_max_coords_x"::Cstring, mpi_max_coords_x::Ref{Clonglong}, PDI_OUT::Cint, "mpi_max_coords_y"::Cstring, mpi_max_coords_y::Ref{Clonglong}, PDI_OUT::Cint, "nx"::Cstring, nx::Ref{Clonglong}, PDI_OUT::Cint, "ny"::Cstring, ny::Ref{Clonglong}, PDI_OUT::Cint, "nb_transported_scalars"::Cstring, phys.nb_transported_scalars::Ref{Clonglong}, PDI_OUT::Cint, "nb_levelsets"::Cstring, phys.nb_levelsets::Ref{Clonglong}, PDI_OUT::Cint, "nstep"::Cstring, nstep::Ref{Clonglong}, PDI_OUT::Cint, # "timestep"::Cstring, nstep::Ref{Clonglong}, PDI_OUT::Cint, # "nb_Navier_slip_BC"::Cstring, num.nNavier::Ref{Clonglong}, PDI_OUT::Cint, C_NULL::Ptr{Cvoid})::Cint print("\n PDI INIT status ",PDI_status) @debug "after PDI_multi_expose" @debug "After full PDI init" # nx_pdi = [nx] # nstep_pdi = [nstep] # local PDI_status = @ccall "libpdi".PDI_multi_expose("test_pycall_write_arr"::Cstring, # "nx_arr"::Cstring, nx_pdi::Ref{Clonglong}, PDI_OUT::Cint, # "nstep_arr"::Cstring, nstep_pdi::Ref{Clonglong}, PDI_OUT::Cint, # C_NULL::Ptr{Cvoid})::Cint # local PDI_status = @ccall "libpdi".PDI_multi_expose("test_pycall_write"::Cstring, # "nx"::Cstring, nx::Ref{Clonglong}, PDI_OUT::Cint, # "nstep"::Cstring, nstep::Ref{Clonglong}, PDI_OUT::Cint, # C_NULL::Ptr{Cvoid})::Cint # local PDI_status = @ccall "libpdi".PDI_multi_expose("test_pycall_write"::Cstring, # "nx"::Cstring, nx::Ref{Clonglong}, PDI_INOUT::Cint, # "nstep"::Cstring, nstep::Ref{Clonglong}, PDI_INOUT::Cint, # C_NULL::Ptr{Cvoid})::Cint # local PDI_status = @ccall "libpdi".PDI_multi_expose("write_pycall"::Cstring, # "nx"::Cstring, nx::Ref{Clonglong}, PDI_OUT::Cint, # "nstep"::Cstring, nstep::Ref{Clonglong}, PDI_OUT::Cint, # C_NULL::Ptr{Cvoid})::Cint end #if io.pdi>0 # arrays to store errors error_list_l1 = zeros(n_cases) error_list_l2 = zeros(n_cases) error_list_linfty = zeros(n_cases) error_list_l1_mixed = zeros(n_cases) error_list_l2_mixed = zeros(n_cases) error_list_linfty_mixed = zeros(n_cases) error_list_l1_full = zeros(n_cases) error_list_l2_full = zeros(n_cases) error_list_linfty_full = zeros(n_cases) cell_volume_list = zeros(n_cases) base_directory = pwd() # print("\n base directory ",base_directory) #region timestep convergence for timestep in timesteps print("\n Timestep ",timestep) timestep_to_string = @sprintf "timestep_%.4e" timestep # print("\n timestep_to_string ",timestep_to_string) timestep_to_string = replace(timestep_to_string, "." => "_") # print("\n timestep_to_string ",timestep_to_string) mkpath(timestep_to_string) cd(timestep_to_string) # print("\n pwd() ",pwd()) #region mesh convergence for (i,study_nb_grid_points) in enumerate(nb_grid_points) mesh_to_string = @sprintf "mesh_%.5i" study_nb_grid_points mesh_ratio = Int((mesh.ymax - mesh.ymin)/ (mesh.xmax - mesh.xmin) ) print("\n mesh ratio ", mesh_ratio ) study_nb_grid_points_x = study_nb_grid_points study_nb_grid_points_y = study_nb_grid_points * mesh_ratio print("\n nx ",study_nb_grid_points_x, " ny ",study_nb_grid_points_y) mkpath(mesh_to_string) cd(mesh_to_string) #delete output.txt: # local PDI_status = @ccall "libpdi".PDI_multi_expose("macro_delete_file"::Cstring,C_NULL::Ptr{Cvoid})::Cint #bug # if study_name !="" # # mkpath(study.change_one_parameter_at_a_time.name) # print("\nstudy ",study.change_one_parameter_at_a_time.name) # end # init regular grid scalar_mesh_x = collect(LinRange(mesh.xmin, mesh.xmax, study_nb_grid_points_x + 1)) scalar_mesh_y = collect(LinRange(mesh.ymin, mesh.ymax, study_nb_grid_points_y + 1)) # # print("\n test juliac") # # juliac_status = @ccall "mylib".load_parameters()::Cint # juliac_status = @ccall "libmylib".load_parameters()::Cint # # print("\n test juliac") # print("\n mu1 mu2 ",phys.mu1," ",typeof(phys.mu1)," ",phys.mu2," ",typeof(phys.mu2)) @debug "Before Numerical" #region test fill struct #aliases from yml to Flower.jl aliases = Dict( :x => :scalar_mesh_x, :y => :scalar_mesh_y, :xcoord => :intfc_x, :ycoord => :intfc_y, :R => :radius, :max_iterations => :max_iter, :save_every => :max_iter, :ϵ => :epsilon, :ϵwall => :epsilon_wall, :nLS => :nb_levelsets, :θd => :temperature0, :β => :beta, :σ => :sigma, :δreinit => :delta_reinit, :n_ext_cl => :n_ext, :timestep_0 => :timestep, :timestep_n => :timestep, :io_pdi => :pdi, # :convection => :convection_mode, #debug ) # fields defined locally, not in yml extra = Dict( :scalar_mesh_x => scalar_mesh_x, :scalar_mesh_y => scalar_mesh_y, # :electrolysis_reaction => Symbol(phys.electrolysis_reaction), # :intfc_x => intfc_x, # :intfc_y => intfc_y, ) default = Numerical{Float64,Int}( x = scalar_mesh_x, y = scalar_mesh_y, timestep_n = timestep, timestep_0 = timestep, electrolysis_reaction_symb = Symbol(phys.electrolysis_reaction), bulk_velocity_symb = Symbol(phys.bulk_velocity), ) # construct default parametric instance with x otherwise L0 and ... not defined in the same way # global num_new = safefill_with_aliases(Numerical{Float64, Int}, sim, phys, io,aliases) # global num_new = safefill_with_aliases_and_extra(Numerical{Float64,Int}, # sim, phys, io, # aliases, # extra) print("\n ns_advection ",(sim.ns_advection ==1)) global num_new = safefill_with_aliases_and_extra_already_init(Numerical{Float64,Int},default, sim, phys, io, aliases, extra) #endregion test fill struct # # global num = Numerical( # # CFL = sim.CFL, # # Re = Re, #in Flower, not real Re # # end_time=phys.end_time, # # x = scalar_mesh_x, # # y = scalar_mesh_y, # # xcoord = phys.intfc_x, # # ycoord = phys.intfc_y, # # case = sim.case, # # R = phys.radius, # # max_iterations = sim.max_iter, # # save_every = sim.max_iter, # # ϵ = sim.epsilon, # # ϵwall = sim.epsilon_wall, # # epsilon_mode = sim.epsilon_mode, # # nLS = phys.nb_levelsets, # # nb_transported_scalars=phys.nb_transported_scalars, # # concentration0=phys.concentration0, # # epsilon_concentration=phys.epsilon_concentration, # # diffusion_coeff=phys.diffusion_coeff, # # temperature0=phys.temperature0, # # i0=phys.i0, # # phi_ele0=phys.phi_ele0, # # phi_ele1=phys.phi_ele1, # # alpha_c=phys.alpha_c, # # alpha_a=phys.alpha_a, # # Ru=phys.Ru, # # Faraday=phys.Faraday, # # MWH2=phys.MWH2, # # θd=phys.temperature0, # # eps=sim.eps, # # mu1=phys.mu1, # # mu2=phys.mu2, # # rho1=phys.rho1, # # rho2=phys.rho2, # # # u_inf = 0.0, # # # v_inf = 0.0, # # pres0=phys.pres0, # # g = phys.g, # # β = phys.beta, # # σ = phys.sigma, # # sigma = phys.sigma, # # reinit_every = sim.reinit_every, # # nb_reinit = sim.nb_reinit, # # δreinit = sim.delta_reinit, # # n_ext_cl = sim.n_ext, # # NB = sim.NB, # # # plot_xscale = io.scale_x, # # timestep_n = timestep, #timestep convergence #sim.timestep_0, # # timestep_0 = timestep, #timestep convergence #sim.timestep_0, # # concentration_check_factor = sim.concentration_check_factor, # # radial_vel_factor = phys.radial_vel_factor, # # debug = sim.debug, # # v_inlet = phys.v_inlet, # # prediction = sim.prediction, # # null_space = sim.null_space, # # io_pdi = io.pdi, # # bulk_conductivity = sim.bulk_conductivity, # # electrical_potential = sim.electrical_potential, # # contact_angle = sim.contact_angle, # # convection_Cdivu = sim.convection_Cdivu, # # convection_mode = sim.convection_mode, # # advection_LS_mode = sim.advection_LS_mode, # # scalar_bc = sim.scalar_bc, # # scalar_scheme = sim.scalar_scheme, # # solver = sim.solver, # # mass_transfer_rate = sim.mass_transfer_rate, # # average_liquid_solid = sim.average_liquid_solid, # # index_phase_change = sim.index_phase_change, # # index_electrolyte = sim.index_electrolyte, # # extend_field = sim.extend_field, # # average_velocity = sim.average_velocity, # # laplacian = sim.laplacian, # # electrical_potential_max_iter = sim.electrical_potential_max_iter, # # electrical_potential_relative_residual = sim.electrical_potential_relative_residual, # # electrical_potential_residual = sim.electrical_potential_residual, # # electrical_potential_nonlinear_solver = sim.electrical_potential_nonlinear_solver, # # electrolysis_reaction = phys.electrolysis_reaction, # # pressure_velocity_coupling = sim.pressure_velocity_coupling, # # pressure_velocity_solver = sim.pressure_velocity_solver, # # solve_solid = sim.solve_solid, # # phase_change_method = sim.phase_change_method, # # one_fluid_model = sim.one_fluid_model, # # smooth_VOF = sim.smooth_VOF, # # surface_tension = sim.surface_tension, # # non_dimensionalize=sim.non_dimensionalize, # # levelset_reinitialize=sim.levelset_reinitialize, # # mu_one_fluid_average = sim.mu_one_fluid_average, # # one_fluid_normal = sim.one_fluid_normal, # # marching_squares_epsilon = sim.marching_squares_epsilon, # # marching_squares_max_iter = sim.marching_squares_max_iter, # # convection = sim.convection_mode, # # nucleation_time = phys.nucleation_time, # # solve_potential = sim.solve_potential, # # solve_species = sim.solve_species, # # kill_dead_cells = sim.kill_dead_cells, # # epsilon_volume_fraction_phase_change = sim.epsilon_volume_fraction_phase_change, # # solve_Navier_Stokes_liquid_phase = sim.solve_Navier_Stokes_liquid_phase, # # mass_transfer_rate_imposed_value = sim.mass_transfer_rate_imposed_value, # # verbosity = sim.verbosity, # # mode_2d = sim.mode_2d, # # mu_cin1 = phys.mu_cin1, # # mu_cin2 = phys.mu_cin2, # # # u_inf = phys.u_inf, # # ) # if (num == num_new) # print("\n New init of num OK") # else # @error("\n init num") # end # function diff_struct(a, b) # @assert typeof(a) == typeof(b) "Types differ" # for name in fieldnames(typeof(a)) # va = getfield(a, name) # vb = getfield(b, name) # if va != vb # println("Field $name differs:") # println(" a.$name = $va") # println(" b.$name = $vb") # end # end # end # diff_struct(num, num_new) # print("\n num x ",num.x) # print("\n num x ",num_new.x) # print("\n num y ",num.y) # print("\n num y ",num_new.y) global num = num_new if num.verbosity >0 print("\n Parameters num ",num) end Broadcast.broadcastable(num::Numerical) = Ref(num) #do not broadcast num @debug "After Numerical" # num.electrolysis_reaction = Symbol(num.electrolysis_reaction) global gp, gu, gv = init_meshes(num) #gp, gu, gv = init_meshes(num) does not work for eval(Meta.parseall(macros.boundaries)) global op, phS, phL = init_fields(num, gp, gu, gv) gp.LS[1].u .= 1.0 #deactivate interface # Init fields eval(Meta.parseall(macros.init_fields)) # Define boundary conditions eval(Meta.parseall(macros.boundaries)) if num.io_pdi>0 # Send meta-data to PDI nx=gp.nx ny=gp.ny try local PDI_status = @ccall "libpdi".PDI_multi_expose("init_PDI"::Cstring, "mpi_coords_x"::Cstring, mpi_coords_x::Ref{Clonglong}, PDI_OUT::Cint, "mpi_coords_y"::Cstring, mpi_coords_x::Ref{Clonglong}, PDI_OUT::Cint, "mpi_max_coords_x"::Cstring, mpi_max_coords_x::Ref{Clonglong}, PDI_OUT::Cint, "mpi_max_coords_y"::Cstring, mpi_max_coords_y::Ref{Clonglong}, PDI_OUT::Cint, "nx"::Cstring, nx::Ref{Clonglong}, PDI_OUT::Cint, "ny"::Cstring, ny::Ref{Clonglong}, PDI_OUT::Cint, "nb_transported_scalars"::Cstring, phys.nb_transported_scalars::Ref{Clonglong}, PDI_OUT::Cint, "nb_levelsets"::Cstring, phys.nb_levelsets::Ref{Clonglong}, PDI_OUT::Cint, "nstep"::Cstring, num.current_iter::Ref{Clonglong}, PDI_OUT::Cint, "nb_Navier_slip_BC"::Cstring, num.nNavier::Ref{Clonglong}, PDI_OUT::Cint, "timestep"::Cstring, timestep::Ref{Cdouble}, PDI_OUT::Cint, C_NULL::Ptr{Cvoid})::Cint catch error print("\n Bug init_PDI \n") print("\n PDI_status ",PDI_status, "\n") print("\n error",error) end @debug "after PDI_multi_expose" @debug "After full PDI init" end #if num.io_pdi>0 # Define interfaces (for bubbles, drops...) eval(Meta.parseall(macros.interface)) # if sim.activate_interface == 1 # gp.LS[1].u .= sqrt.((gp.x .- phys.intfc_x).^2 + (gp.y .- phys.intfc_y).^2) - phys.radius * ones(gp) # #modify velocity field near interface # su = sqrt.((gv.x .- phys.intfc_x).^2 .+ (gv.y .- phys.intfc_y).^2) # R1 = phys.radius + 3.0*num.Δ # bl = 4.0 # for II in gv.ind.all_indices # if su[II] <= R1 # phL.v[II] = 0.0 # # elseif su[II] > R1 # # uL[II] = tanh(bl*(su[II]-R1)) # end # end # elseif sim.activate_interface == -1 # gp.LS[1].u .= sqrt.((gp.x .- phys.intfc_x).^2 + (gp.y .- phys.intfc_y).^2) - phys.radius * ones(gp) # gp.LS[1].u .*= -1.0 # else # gp.LS[1].u .= 1.0 # end # test_LS(gp) # Create segments from interface # x,y,field,connectivities,num_vtx = convert_interfacial_D_to_segments(num,gp,phL.T,1) # print("\n number of interface points ", num_vtx) # # print("\n x",x) # # print("\n x",y) # # print("\n x",field) # print("\n x",connectivities) # print("\n x",num_vtx) # printstyled(color=:green, @sprintf "\n Initialisation0 \n") # print_electrolysis_statistics(num,gp,phL) if sim.time_scheme == "FE" time_scheme = FE else time_scheme = CN end if num.io_pdi>0 try # printstyled(color=:red, @sprintf "\n before pdi \n") # printstyled(color=:red, @sprintf "\n PDI test \n" ) # if num.solve_electrical_potential>0 # phi_array=phL.phi_ele #do not transpose since python row major # compute_grad_phi_ele!(num, grid, grid_u, grid_v, phL, phS, op.opC_pL, op.opC_pS) #TODO current # Eus,Evs = interpolate_grid_liquid(grid,grid_u,grid_v,phL.Eu, phL.Ev) # us,vs = interpolate_grid_liquid(grid,grid_u,grid_v,phL.u,phL.v) # print("\n before write \n ") #TODO "levelset_p_tot"::Cstring, gp.LS[2].u::Ptr{Cdouble}, PDI_OUT::Cint, iLSpdi = 1 # all LS iLS = 1 # or all 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_linftyp %.5i \n" num.current_iter) #print("\n size LS wall ", size( gp.LS[2].u)) LStable = zeros(gp) if phys.nb_levelsets>1 LStable = gp.LS[2].u end # Compute electrical current, interpolate velocity on scalar grid # compute_grad_phi_ele!(num, gp, gu, gv, phL, phS, op.opC_pL, op.opC_pS) #TODO current # compute_grad_phi_ele!(num, grid, grid_u, grid_v, 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 us=zeros(gp) #TODO allocate only once vs=zeros(gp) # #store in us, vs instead of Eus, Evs # interpolate_grid_liquid!(gp,gu,gv,phL.Eu, phL.Ev,us,vs) # #TODO i_current_mag need cond # @ccall "libpdi".PDI_multi_expose("write_data_elec"::Cstring, # "i_current_x"::Cstring, us::Ptr{Cdouble}, PDI_OUT::Cint, # "i_current_y"::Cstring, vs::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!(gp,gu,gv,phL.u,phL.v,us,vs) current_radius = phys.radius PDI_status = @ccall "libpdi".PDI_multi_expose("write_initialization"::Cstring, "nstep"::Cstring, num.current_iter::Ref{Clonglong}, PDI_OUT::Cint, "time"::Cstring, phys_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, gp.LS[iLSpdi].u::Ptr{Cdouble}, PDI_OUT::Cint, "levelset_u"::Cstring, gu.LS[iLSpdi].u::Ptr{Cdouble}, PDI_OUT::Cint, "levelset_v"::Cstring, gv.LS[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 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 velocity_y_bubble = 0.0 .* gp.LS[end].geoS.dcap[:,:,5] volume_fraction = gp.LS[end].geoS.dcap[:,:,5] center_of_mass_x = 0.0 center_of_mass_y = 0.0 circularity = 1.0 rise_velocity_y = 0.0 iLSpdi = 1 velocity_y = velocity_y_bubble area_liq = 0.0 area_gaz = 0.0 # PDI_status = @ccall "libpdi".PDI_multi_expose("write_postprocessing_rising_bubble"::Cstring, # "nstep"::Cstring, num.current_iter ::Ref{Clonglong}, PDI_OUT::Cint, # "center_of_mass_x"::Cstring, center_of_mass_x::Ref{Cdouble}, PDI_OUT::Cint, # "center_of_mass_y"::Cstring, center_of_mass_y::Ref{Cdouble}, PDI_OUT::Cint, # "rise_velocity_y"::Cstring, rise_velocity_y::Ref{Cdouble}, PDI_OUT::Cint, # "circularity"::Cstring, circularity::Ref{Cdouble}, PDI_OUT::Cint, # "velocity_y_bubble"::Cstring, velocity_y_bubble::Ptr{Cdouble}, PDI_OUT::Cint, # "area_gaz"::Cstring, area_gaz::Ref{Cdouble}, PDI_OUT::Cint, # "area_liq"::Cstring, area_liq::Ref{Cdouble}, PDI_OUT::Cint, # C_NULL::Ptr{Cvoid})::Cint # PDI_status = @ccall "libpdi".PDI_multi_expose("post_processing_rising_bubble_new"::Cstring, # "nstep"::Cstring, num.current_iter ::Ref{Clonglong}, PDI_OUT::Cint, # "velocity_y"::Cstring, velocity_y_bubble::Ptr{Cdouble}, PDI_OUT::Cint, # "volume_fraction"::Cstring, volume_fraction::Ptr{Cdouble}, PDI_OUT::Cint, # # "volume_liq_cell"::Cstring, grid_p.LS[end].geoL.dcap[:,:,5]::Ptr{Cdouble}, PDI_OUT::Cint, #geoS for bubble phase # # "volume_cell"::Cstring, grid_p.LS[end].geoS.dcap[:,:,5]::Ptr{Cdouble}, PDI_OUT::Cint, #geoS for bubble phase # "mesh_p_x"::Cstring, gp.x::Ptr{Cdouble}, PDI_OUT::Cint, # "mesh_p_y"::Cstring, gp.y::Ptr{Cdouble}, PDI_OUT::Cint, # "dcap_1"::Cstring, gp.LS[iLSpdi].geoS.dcap[:,:,1]::Ptr{Cdouble}, PDI_OUT::Cint, #geoS for bubble phase # "dcap_2"::Cstring, gp.LS[iLSpdi].geoS.dcap[:,:,2]::Ptr{Cdouble}, PDI_OUT::Cint, #geoS for bubble phase # "dcap_3"::Cstring, gp.LS[iLSpdi].geoS.dcap[:,:,3]::Ptr{Cdouble}, PDI_OUT::Cint, #geoS for bubble phase # "dcap_4"::Cstring, gp.LS[iLSpdi].geoS.dcap[:,:,4]::Ptr{Cdouble}, PDI_OUT::Cint, #geoS for bubble phase # C_NULL::Ptr{Cvoid})::Cint PDI_status = @ccall "libpdi".PDI_multi_expose("post_processing_rising_bubble"::Cstring, "nstep"::Cstring, num.current_iter ::Ref{Clonglong}, PDI_OUT::Cint, "velocity_y"::Cstring, velocity_y::Ptr{Cdouble}, PDI_OUT::Cint, "volume_fraction"::Cstring, volume_fraction::Ptr{Cdouble}, PDI_OUT::Cint, "volume_liq_cell"::Cstring, gp.LS[iLSpdi].geoL.dcap[:,:,5]::Ptr{Cdouble}, PDI_OUT::Cint, #geoS for bubble phase "volume_cell"::Cstring, gp.LS[iLSpdi].geoS.dcap[:,:,5]::Ptr{Cdouble}, PDI_OUT::Cint, #geoS for bubble phase "mesh_p_x"::Cstring, gp.x::Ptr{Cdouble}, PDI_OUT::Cint, "mesh_p_y"::Cstring, gp.y::Ptr{Cdouble}, PDI_OUT::Cint, "dcap_1"::Cstring, gp.LS[iLSpdi].geoS.dcap[:,:,1]::Ptr{Cdouble}, PDI_OUT::Cint, #geoS for bubble phase "dcap_2"::Cstring, gp.LS[iLSpdi].geoS.dcap[:,:,2]::Ptr{Cdouble}, PDI_OUT::Cint, #geoS for bubble phase "dcap_3"::Cstring, gp.LS[iLSpdi].geoS.dcap[:,:,3]::Ptr{Cdouble}, PDI_OUT::Cint, #geoS for bubble phase "dcap_4"::Cstring, gp.LS[iLSpdi].geoS.dcap[:,:,4]::Ptr{Cdouble}, PDI_OUT::Cint, #geoS for bubble phase C_NULL::Ptr{Cvoid})::Cint # if num.io_pdi>0 # iLSpdi = 1 # TODO all grid.LS # PDI_status = @ccall "libpdi".PDI_multi_expose("write_capacities"::Cstring, # "dcap"::Cstring, permutedims(gp.LS[iLSpdi].geoL.dcap, (3, 1, 2))::Ptr{Cdouble}, PDI_OUT::Cint, # C_NULL::Ptr{Cvoid})::Cint # # try # # iLSpdi = 1 # TODO all grid.LS # # PDI_status = @ccall "libpdi".PDI_multi_expose("write_capacities"::Cstring, # # "dcap"::Cstring, gp.LS[iLSpdi].geoL.dcap'::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 # printstyled(color=:red, @sprintf "\n after pdi \n") # printstyled(color=:red, @sprintf "\n before run_forward \n") # print("\n BC_uL ",BC_uL) run_forward!( num, gp, gu, gv, op, phS, phL; periodic_x = (sim.periodic_x == 1), periodic_y = (sim.periodic_y == 1), BC_uL = BC_uL, BC_uS=BC_uS, BC_vL = BC_vL, BC_vS=BC_vS, BC_pL = BC_pL, BC_pS=BC_pS, BC_u = BC_u, BC_int = BC_int, BC_trans_scal=BC_trans_scal, BC_phi_ele = BC_phi_ele, auto_reinit = sim.auto_reinit, time_scheme = time_scheme, electrolysis = true, navier_stokes = true, ns_advection = (sim.ns_advection ==1), ns_liquid_phase = (sim.solve_Navier_Stokes_liquid_phase == 1), verbose = true, show_every = sim.show_every, electrolysis_convection = (sim.electrolysis_convection ==1), electrolysis_liquid_phase = true, electrolysis_phase_change_case = sim.electrolysis_phase_change_case, imposed_velocity = sim.imposed_velocity, adapt_timestep_mode = sim.adapt_timestep_mode,#1, non_dimensionalize=sim.non_dimensionalize, mode_2d = sim.mode_2d, breakup = sim.breakup, ) @debug "After run" # #cut small cells for error # number_small_cells_for_error = 0 # cutoff_for_error_volume = 1e-12 #TODO # for II in gp.ind.all_indices # if gp.LS[1].geoL.cap[II,5] < cutoff_for_error_volume # number_small_cells_for_error += 1 # Tana[II] = 0.0 # T[II] = 0.0 # end # end # printstyled(color=:green, @sprintf "\n number_small_cells_for_error %.3i \n" number_small_cells_for_error) if study.compute_errors != "None" #"Poiseuille" LIQUID = gp.ind.all_indices[gp.LS[1].geoL.cap[:,:,5] .> (1-1e-16)] MIXED = gp.ind.all_indices[gp.LS[1].geoL.cap[:,:,5] .<= (1-1e-16) .&& gp.LS[1].geoL.cap[:,:,5] .> 1e-16] if study.compute_errors == "Poiseuille" norm_all = relative_errors(phL.v, vPoiseuille, vcat(LIQUID, MIXED), gp.LS[1].geoL.cap[:,:,5], num.Δ) norm_mixed = relative_errors(phL.v, vPoiseuille, MIXED, gp.LS[1].geoL.cap[:,:,5], num.Δ) norm_full = relative_errors(phL.v, vPoiseuille, LIQUID, gp.LS[1].geoL.cap[:,:,5], num.Δ) elseif study.compute_errors == "diffusion_full_cell" # Parameters L = mesh.xmax-mesh.xmin D = num.diffusion_coeff[2] #3.2e-9 # diffusivity analytical_nx = 1000 #accumulating error if n too small so if N big, accumulating too much if nx=100 and N=1000 diffusion_time_scale = L^2 / D println("diffusion time scale", diffusion_time_scale) # ele_current_i = minimum(i_current_mag from gradient) # print("i mag min ", ele_current_i) ele_current_i = -1.59e4 F1 = ele_current_i / (2 * num.Faraday * D) #BC t_max = diffusion_time_scale println("diffusion_time_scale", diffusion_time_scale) analytical_dx = L / analytical_nx # Spatial step size analytical_x = collect(0:analytical_dx:L) # Spatial domain c0 = num.concentration0[2] println("max c without convection", c0 - F1 * L / 2, c0 + F1 * L / 2) max_nb_Fourier_series = 1000 # Special function u1 function full_cell_u1(x, t) return F1 * x # return (F2-F1)/(2*L) + F1 * x + D*(F2-F1)/L*t #if different left right end # Initial condition f(x) function full_cell_f(x) return c0 - full_cell_u1(x, 0) # Example initial condition c0-u1? end # Coefficients A_n function full_cell_A_n(n, L, F1, dx,c0) # integral = sum(full_cell_f(xi) * cos(n * π * xi / L) for xi in x) * dx # return (2 / L) * integral #* cos(n * π * x / L) if n == 0 return -F1*L +2*c0 #-c0*L*2/L #analytical if just F1 * x else return (-2*F1*L)/(n*π)^2*((-1)^n-1) #analytical if just F1 * x end end # Solution u2 function full_cell_u2(x, t, L, D, max_nb_Fourier_series,c0,analytical_dx) result = sum(full_cell_A_n(n, L, F1, analytical_dx,c0) * cos(n * π * x / L) * exp(-D * (n * π / L)^2 * t) for n in 1:max_nb_Fourier_series) result += full_cell_A_n(0, L, F1, analytical_dx,c0) * cos(0 * π * x / L) * exp(-D * (0 * π / L)^2 * t) / 2 #case 0 special int cos(0) cos(0)dx = L return result end # Total solution u function full_cell_u(x, t, L, D, max_nb_Fourier_series,c0,analytical_dx) return full_cell_u1(x, t) + full_cell_u2(x, t, L, D, max_nb_Fourier_series,c0,analytical_dx) end # print("\n param ",gp.x[2,:]," | ",num.time," | ",L," | ",num.diffusion_coeff[2]," | ",max_nb_Fourier_series," | ",c0," | ",analytical_dx) #TODO NX vs n ... need to interpolate # concentration_profile = full_cell_u.(analytical_x,num.time,L,num.diffusion_coeff[2],max_nb_Fourier_series,c0,analytical_dx) concentration_profile = full_cell_u.(gp.x[2,:],num.time,L,num.diffusion_coeff[2],max_nb_Fourier_series,c0,analytical_dx) print("\n concentration profile ",phL.trans_scal[2,:,2]) print("\n concentration profile ",concentration_profile) # norm_all = relative_errors(phL.trans_scal[:,:,2], concentration_profile, vcat(LIQUID, MIXED), gp.LS[1].geoL.cap[:,:,5], num.Δ) # norm_mixed = relative_errors(phL.trans_scal[:,:,2], concentration_profile, MIXED, gp.LS[1].geoL.cap[:,:,5], num.Δ) # norm_full = relative_errors(phL.trans_scal[:,:,2], concentration_profile, LIQUID, gp.LS[1].geoL.cap[:,:,5], num.Δ) l1,l2,linfty = relative_errors(phL.trans_scal[:,:,2], concentration_profile, vcat(LIQUID, MIXED), gp.LS[1].geoL.cap[:,:,5], num.Δ) l1_mixed,l2_mixed,linfty_mixed = relative_errors(phL.trans_scal[:,:,2], concentration_profile, MIXED, gp.LS[1].geoL.cap[:,:,5], num.Δ) l1_full,l2_full,linfty_full = relative_errors(phL.trans_scal[:,:,2], concentration_profile, LIQUID, gp.LS[1].geoL.cap[:,:,5], num.Δ) end # error_list_l1[i] = norm_all[1] # error_list_l2[i] = norm_all[2] # error_list_linfty[i] = norm_all[3] # error_list_l1_mixed[i] = norm_mixed[1] # error_list_l2_mixed[i] = norm_mixed[2] # error_list_linfty_mixed[i] = norm_mixed[3] # error_list_l1_full[i] = norm_full[1] # error_list_l2_full[i] = norm_full[2] # error_list_linfty_full[i] = norm_full[3] error_list_l1[i] = l1 error_list_l2[i] = l2 error_list_linfty[i] = linfty error_list_l1_mixed[i] = l1_mixed error_list_l2_mixed[i] = l2_mixed error_list_linfty_mixed[i] = linfty_mixed error_list_l1_full[i] = l1_full error_list_l2_full[i] = l2_full error_list_linfty_full[i] = linfty_full cell_volume_list[i] = minimum(gp.LS[1].geoL.dcap[:,:,5]) min_cell_volume = minimum(gp.LS[1].geoL.dcap[:,:,5]) if study.compute_errors != "None" #"Poiseuille" local PDI_status = @ccall "libpdi".PDI_multi_expose("convergence_study_iter"::Cstring, "study_nb_grid_points"::Cstring, study_nb_grid_points::Ref{Clong}, PDI_OUT::Cint, "study_timestep"::Cstring, timestep::Ref{Cdouble}, PDI_OUT::Cint, # "cell_volume"::Cstring, cell_volume_list::Ptr{Cdouble}, PDI_OUT::Cint, "study_l1_rel_error"::Cstring, l1::Ref{Cdouble}, PDI_OUT::Cint, "study_l2_rel_error"::Cstring, l2::Ref{Cdouble}, PDI_OUT::Cint, "study_linfty_rel_error"::Cstring, linfty::Ref{Cdouble}, PDI_OUT::Cint, "study_l1_rel_error_full_cells"::Cstring, l1_full::Ref{Cdouble}, PDI_OUT::Cint, "study_l2_rel_error_full_cells"::Cstring, l2_full::Ref{Cdouble}, PDI_OUT::Cint, "study_linfty_rel_error_full_cells"::Cstring, linfty_full::Ref{Cdouble}, PDI_OUT::Cint, "study_l1_rel_error_partial_cells"::Cstring, l1_mixed::Ref{Cdouble}, PDI_OUT::Cint, "study_l2_rel_error_partial_cells"::Cstring, l2_mixed::Ref{Cdouble}, PDI_OUT::Cint, "study_linfty_rel_error_partial_cells"::Cstring, linfty_mixed::Ref{Cdouble}, PDI_OUT::Cint, "domain_length"::Cstring, L0::Ref{Cdouble}, PDI_OUT::Cint, "min_cell_volume"::Cstring, min_cell_volume::Ref{Cdouble}, PDI_OUT::Cint, C_NULL::Ptr{Cvoid})::Cint end end current_directory = pwd() if current_directory == base_directory*"/"*timestep_to_string*"/"*mesh_to_string cd("..") else print("\n Error in mesh subdirectory",current_directory," not ",base_directory*"/"*timestep_to_string*"/"*mesh_to_string) exit() end end #endregion mesh convergence current_directory = pwd() if current_directory == base_directory*"/"*timestep_to_string cd("..") else print("\n Error in timestep subdirectory",current_directory," not ",base_directory*"/"*timestep_to_string) exit() end end #endregion timestep convergence min_cell_volume = minimum(gp.LS[1].geoL.cap[:,:,5]) print("\n min_cell_volume ",min_cell_volume," type ",typeof(min_cell_volume)) if study.compute_errors != "None" #"Poiseuille" local PDI_status = @ccall "libpdi".PDI_multi_expose("convergence_study"::Cstring, "n_tests"::Cstring, n_cases::Ref{Clonglong}, PDI_OUT::Cint, "nx_list"::Cstring, nb_grid_points::Ptr{Clonglong}, PDI_OUT::Cint, "cell_volume_list"::Cstring, cell_volume_list::Ptr{Cdouble}, PDI_OUT::Cint, "l1_rel_error"::Cstring, error_list_l1::Ptr{Cdouble}, PDI_OUT::Cint, "l2_rel_error"::Cstring, error_list_l2::Ptr{Cdouble}, PDI_OUT::Cint, "linfty_rel_error"::Cstring, error_list_linfty::Ptr{Cdouble}, PDI_OUT::Cint, "l1_rel_error_full_cells"::Cstring, error_list_l1_full::Ptr{Cdouble}, PDI_OUT::Cint, "l2_rel_error_full_cells"::Cstring, error_list_l2_full::Ptr{Cdouble}, PDI_OUT::Cint, "linfty_rel_error_full_cells"::Cstring, error_list_linfty_full::Ptr{Cdouble}, PDI_OUT::Cint, "l1_rel_error_partial_cells"::Cstring, error_list_l1_mixed::Ptr{Cdouble}, PDI_OUT::Cint, "l2_rel_error_partial_cells"::Cstring, error_list_l2_mixed::Ptr{Cdouble}, PDI_OUT::Cint, "linfty_rel_error_partial_cells"::Cstring, error_list_linfty_mixed::Ptr{Cdouble}, PDI_OUT::Cint, "domain_length"::Cstring, L0::Ref{Cdouble}, PDI_OUT::Cint, "min_cell_volume"::Cstring, min_cell_volume::Ref{Cdouble}, PDI_OUT::Cint, C_NULL::Ptr{Cvoid})::Cint end if io.pdi>0 try # local PDI_status = @ccall "libpdi".PDI_event("close_pycall"::Cstring)::Cint # PDI_event("finalization"); # local PDI_status = @ccall "libpdi".PDI_event("finalization"::Cstring)::Cint local PDI_status = @ccall "libpdi".PDI_finalize()::Cint # printstyled(color=:red, @sprintf "\n PDI end\n" ) catch error printstyled(color=:red, @sprintf "\n PDI error \n") print(error) end end #if io.pdi>0 printstyled(color=:red, @sprintf "\n After PDI \n") #Tests if haskey(macros,"test_end") eval(Meta.parseall(macros.test_end)) end