corr

Signed-off-by: Jean-Michel Batto <jmbatto@eldarsoft.com>

Dockerfile

@@ -130,6 +130,7 @@ ENV GKSwstype=100
 ENV JULIA_PKG_PRECOMPILE_AUTO=0
 ENV JULIA_PKG_USE_CLI_GIT=true
 ENV DISPLAY=host.docker.internal:0.0
+ENV PYTHONPATH=/usr/local/lib/python3/dist-packages
 
 # Switch to non-root user for Package Installation
 # Ensures ~/.julia permissions are correctly set for the user 'coder'.

convergence.jl

@@ -0,0 +1,1043 @@
+# 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

redo-20260123.txt

@@ -132,3 +132,60 @@ include("test_functions.jl")
 
 include("post_processing.jl")
 
+[06:59:48][PDI] *** error: Error while triggering event `write_initialization': ModuleNotFoundError: No module named 'pdi'
+
+ENV PYTHONPATH=/usr/local/lib/python3/dist-packages
+
+ERROR: LoadError: UndefVarError: `BC_uL` not defined in `Main`
+The binding may be too new: running in world age 38879, while current world is 38905.
+Stacktrace:
+ [1] top-level scope
+   @ /usr/local/var/Flower.jl/examples/convergence.jl:766
+ [2] include(mod::Module, _path::String)
+   @ Base ./Base.jl:306
+ [3] exec_options(opts::Base.JLOptions)
+   @ Base ./client.jl:317
+ [4] _start()
+   @ Base ./client.jl:550
+in expression starting at /usr/local/var/Flower.jl/examples/convergence.jl:228
+
+erreur ligne 766 convergence.jl
+dans le run forward
+
+# CORRECTION WORLD AGE :
+        # On utilise getfield(Main, :NomVariable) pour forcer la récupération dynamique
+        # des variables créées par eval() dans la boucle courante.
+        run_forward!(
+            num, gp, gu, gv, op, phS, phL;
+            periodic_x = (sim.periodic_x == 1),
+            periodic_y = (sim.periodic_y == 1),
+            # --- CORRECTION DES LIENS DYNAMIQUES ICI ---
+            BC_uL = getfield(Main, :BC_uL),
+            BC_uS = getfield(Main, :BC_uS),
+            BC_vL = getfield(Main, :BC_vL),
+            BC_vS = getfield(Main, :BC_vS),
+            BC_pL = getfield(Main, :BC_pL),
+            BC_pS = getfield(Main, :BC_pS),
+            BC_u  = getfield(Main, :BC_u),
+            BC_int = getfield(Main, :BC_int),
+            BC_trans_scal = getfield(Main, :BC_trans_scal),
+            BC_phi_ele = getfield(Main, :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,
+            non_dimensionalize=sim.non_dimensionalize,
+            mode_2d = sim.mode_2d,
+            breakup = sim.breakup,    
+        ) breakup = sim.breakup,    
+        )