added bug scripts

scripts/Embedded/Bugs/odes_skeleton_distributed_bug.jl

@@ -0,0 +1,52 @@
+using Test
+
+using GridapTopOpt
+using Gridap, Gridap.Geometry, Gridap.Adaptivity
+using GridapEmbedded
+
+using GridapDistributed, PartitionedArrays
+
+parts = (2,2)
+ranks = DebugArray(LinearIndices((prod(parts),)))
+
+n = 200
+order = 1
+
+_model = CartesianDiscreteModel(ranks,parts,(0,1,0,1),(n,n))
+base_model = UnstructuredDiscreteModel(_model)
+ref_model = refine(base_model, refinement_method = "barycentric")
+model = Adaptivity.get_model(ref_model)
+el_Δ = get_el_Δ(_model)
+h = maximum(el_Δ)
+h_refine = maximum(el_Δ)/2
+
+Ω = Triangulation(model)
+dΩ = Measure(Ω,2*order)
+reffe_scalar = ReferenceFE(lagrangian,Float64,order)
+V_φ = TestFESpace(model,reffe_scalar)
+Γg = SkeletonTriangulation(get_triangulation(V_φ))
+dΓg = Measure(Γg,2order)
+n_Γg = get_normal_vector(Γg)
+
+φh = interpolate(x->-sqrt((x[1]-0.5)^2+(x[2]-0.5)^2)+0.25,V_φ)
+
+velh = interpolate(x->-1,V_φ)
+
+γg = 0.1
+dt = GridapTopOpt._compute_Δt(h,0.1,get_free_dof_values(velh))
+ϵ = 1e-20
+
+ode_ls = LUSolver()
+ode_nl = NLSolver(ode_ls, show_trace=false, method=:newton, iterations=10)
+oode_solver = RungeKutta(ode_nl, ode_ls, 0.1, :DIRK_CrankNicolson_2_2)
+
+β = velh*∇(φh)/(ϵ + norm ∘ ∇(φh))
+stiffness(t,u,v) = ∫((β ⋅ ∇(u)) * v)dΩ + ∫(γg*h^2*jump(∇(u) ⋅ n_Γg)*jump(∇(v) ⋅ n_Γg))dΓg
+mass(t, ∂ₜu, v) = ∫(∂ₜu * v)dΩ
+forcing(t,v) = ∫(0v)dΩ
+Ut_φ = TransientTrialFESpace(V_φ)
+ode_op = TransientLinearFEOperator((stiffness,mass),forcing,Ut_φ,V_φ;
+  constant_forms=(false,true))
+ode_sol = solve(oode_solver,ode_op,0.0,dt*10,φh)
+
+march = Base.iterate(ode_sol)

src/LevelSetEvolution/UnfittedEvolution/StabilisedReinit.jl

@@ -3,24 +3,24 @@ abstract type StabilisationMethod end
 """
     mutable struct StabilisedReinit{V,M} <: Reinitialiser
 
-Stabilised FE method for level-set reinitialisation. Artificial viscosity 
-approach (`ArtificialViscosity`) based on Mallon et al. (2023). DOI: `10.48550/arXiv.2303.13672`. 
+Stabilised FE method for level-set reinitialisation. Artificial viscosity
+approach (`ArtificialViscosity`) based on Mallon et al. (2023). DOI: `10.48550/arXiv.2303.13672`.
 Interior jump penalty approach (`InteriorPenalty`) adapted from that work and replaces
-the artifical viscosity term with an interior jump penalty term.  
+the artifical viscosity term with an interior jump penalty term.
 
 # Parameters
 - `nls::NonlinearSolver`: Nonlinear solver for solving the reinitialisation equation
-- `stabilisation_method::A`: A `StabilisationMethod` method for stabilising the problem 
+- `stabilisation_method::A`: A `StabilisationMethod` method for stabilising the problem
 - `Ωs::B`: `EmbeddedCollection` holding updatable triangulation and measures from GridapEmbedded
 - `dΩ_bg::D`: Background mesh measure
 - `space::C`: FE space for level-set function
-- `assembler::Assembler`: Assembler for LS FE space 
+- `assembler::Assembler`: Assembler for LS FE space
 - `params::E`: Tuple of Nitsche parameter `γd` and mesh size `h`
 - `cache`: Cache for evolver, initially `nothing`.
 
 # Note
-- We expect the EmbeddedCollection `Ωs` to contain `:dΓ`. If this is not 
-  available we add it to the recipe list in `Ωs` and a warning will appear. 
+- We expect the EmbeddedCollection `Ωs` to contain `:dΓ`. If this is not
+  available we add it to the recipe list in `Ωs` and a warning will appear.
 """
 mutable struct StabilisedReinit{A,B,C,D} <: Reinitialiser
   nls::NonlinearSolver
@@ -36,16 +36,16 @@ mutable struct StabilisedReinit{A,B,C,D} <: Reinitialiser
       nls = NewtonSolver(LUSolver();maxiter=50,rtol=1.e-14,verbose=true),
       assembler=SparseMatrixAssembler(V_φ,V_φ),
       stabilisation_method::A = InteriorPenalty(V_φ)) where {A,B,C}
-    
+
     if !(:dΓ ∈ keys(Ωs.objects))
-      @warn "Expected measure ':dΓ' not found in the 
-      EmbeddedCollection. This has been added to the recipe list. 
-      
+      @warn "Expected measure ':dΓ' not found in the
+      EmbeddedCollection. This has been added to the recipe list.
+
       Ensure that you are not using ':dΓ' under a different
       name to avoid additional computation for cutting."
       function dΓ_recipe(cutgeo)
         Γ = DifferentiableTriangulation(EmbeddedBoundary(cutgeo),V_φ)
-        (; 
+        (;
           :dΓ => Measure(Γ,2get_order(V_φ))
         )
       end
@@ -92,7 +92,7 @@ function solve!(s::StabilisedReinit,φh,nls_cache)
   return φh
 end
 
-struct ArtificialViscosity <: StabilisationMethod 
+struct ArtificialViscosity <: StabilisationMethod
   stabilisation_coefficent::Number
 end
 
@@ -116,7 +116,7 @@ function get_residual_and_jacobian(s::StabilisedReinit{ArtificialViscosity},φh)
 end
 
 struct InteriorPenalty <: StabilisationMethod
-  dΛ::Measure
+  dΛ
 end
 function InteriorPenalty(V_φ::FESpace)
   Λ = SkeletonTriangulation(get_triangulation(V_φ))