From 618409684280085c35e0066c1ed4bb4c3e02c9ad Mon Sep 17 00:00:00 2001 From: Z J Wegert <60646897+zjwegert@users.noreply.github.com> Date: Wed, 13 Nov 2024 15:22:09 +1100 Subject: [PATCH] Added UnfittedEvolution to src and tests --- .../{ => deprecated}/evolve.jl | 216 +++++++++--------- .../Unfitted_Evolve&Reinit/evolution_test.jl | 88 ------- .../Unfitted_Evolve&Reinit/reinit_test.jl | 56 ----- src/GridapTopOpt.jl | 7 + src/LevelSetEvolution/LevelSetEvolution.jl | 3 +- .../UnfittedEvolution/CutFEMEvolve.jl | 92 +++++--- .../UnfittedEvolution/StabilisedReinit.jl | 10 +- ...ArtificialViscosityStabilisedReinitTest.jl | 54 +++++ .../CutFEMEvolveTest.jl | 54 +++++ .../InteriorPenaltyStabilisedReinitTest.jl | 55 +++++ test/seq/UnfittedEvolutionTests/runtests.jl | 9 + test/seq/runtests.jl | 1 + 12 files changed, 356 insertions(+), 289 deletions(-) rename scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/{ => deprecated}/evolve.jl (96%) delete mode 100644 scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/evolution_test.jl delete mode 100644 scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/reinit_test.jl create mode 100644 test/seq/UnfittedEvolutionTests/ArtificialViscosityStabilisedReinitTest.jl create mode 100644 test/seq/UnfittedEvolutionTests/CutFEMEvolveTest.jl create mode 100644 test/seq/UnfittedEvolutionTests/InteriorPenaltyStabilisedReinitTest.jl create mode 100644 test/seq/UnfittedEvolutionTests/runtests.jl diff --git a/scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/evolve.jl b/scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/deprecated/evolve.jl similarity index 96% rename from scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/evolve.jl rename to scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/deprecated/evolve.jl index 3ac907fb..528376e7 100644 --- a/scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/evolve.jl +++ b/scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/deprecated/evolve.jl @@ -1,108 +1,108 @@ -# Based on: -# @article{ -# Burman_Elfverson_Hansbo_Larson_Larsson_2018, -# title={Shape optimization using the cut finite element method}, -# volume={328}, -# ISSN={00457825}, -# DOI={10.1016/j.cma.2017.09.005}, -# journal={Computer Methods in Applied Mechanics and Engineering}, -# author={Burman, Erik and Elfverson, Daniel and Hansbo, Peter and Larson, Mats G. and Larsson, Karl}, -# year={2018}, -# month=jan, -# pages={242–261}, -# language={en} -# } - -using GridapTopOpt - -using Gridap - -using GridapEmbedded -using GridapEmbedded.LevelSetCutters -using Gridap.Geometry, Gridap.FESpaces, Gridap.CellData -import Gridap.Geometry: get_node_coordinates, collect1d - -include("../../differentiable_trians.jl") - -order = 1 -n = 101 - -_model = CartesianDiscreteModel((0,1,0,1),(n,n)) -cd = Gridap.Geometry.get_cartesian_descriptor(_model) -h = maximum(cd.sizes) - -model = simplexify(_model) -Ω = Triangulation(model) -dΩ = Measure(Ω,2*order) -reffe_scalar = ReferenceFE(lagrangian,Float64,order) -V_φ = TestFESpace(model,reffe_scalar) - -φh = interpolate(x->-sqrt((x[1]-0.5)^2+(x[2]-0.5)^2)+0.25,V_φ) # <- Already SDF -# φh = interpolate(x->cos(2π*x[1])*cos(2π*x[2])-0.11,V_φ) # This needs to reinitialised, can use `reinit_based_on_connor_paper.jl` - -geo = DiscreteGeometry(φh,model) -cutgeo = cut(model,geo) -Γ = EmbeddedBoundary(cutgeo) -dΓ = Measure(Γ,2*order) - -F = EmbeddedFacetDiscretization(LevelSetCutter(),model,geo) -FΓ = SkeletonTriangulation(F) -dFΓ = Measure(FΓ,2*order) - -dt = 0.01 -Tn = 5*dt -sysslvr_l = LUSolver() -sysslvr_nl = NLSolver(sysslvr_l, show_trace=true, method=:newton, iterations=10) -ode_solver = RungeKutta(sysslvr_nl, sysslvr_l, dt, :DIRK_CrankNicolson_2_2) - -c = 0.1 - -Ut_φ = TransientTrialFESpace(V_φ,t -> (x->-1)) -n = get_normal_vector(FΓ) -velh = interpolate(x->-1,V_φ) -d1(∇u) = 1 / ( ϵ + norm(∇u) ) -_n(∇u) = ∇u/(10^-20+norm(∇u)) -β = velh*∇(φh)/(10^-20+norm ∘ ∇(φh)) -stiffness(t,u,v) = ∫((β ⋅ ∇(u)) * v)dΩ + ∫(c*h^2*jump(∇(u) ⋅ n)*jump(∇(v) ⋅ n))dFΓ -mass(t, ∂ₜu, v) = ∫(∂ₜu * v)dΩ -forcing(t,v) = ∫(0v)dΩ -op = TransientLinearFEOperator((stiffness,mass),forcing,Ut_φ,V_φ,constant_forms=(true,true)) -uht = solve(ode_solver,op,0.0,Tn,φh) -for (t,uh) in uht - if t ≈ Tn - copyto!(get_free_dof_values(φh),get_free_dof_values(uh)) - end -end - -# You can compare the zero-iso curve generated by the following -φh_expected_contour = interpolate(x->-sqrt((x[1]-0.5)^2+(x[2]-0.5)^2)+0.25+Tn,V_φ) - -writevtk( - Ω,"results/test_evolve", - cellfields=["φh"=>φh,"φh_expected_contour"=>φh_expected_contour,"|∇φh|"=>norm ∘ ∇(φh)] -) - - -sol,state = Base.iterate(uht) -_tF = state[2][1] -_stateF = state[2][2] -_state0 = state[2][3] -_uF = state[2][4] -_odecache = state[2][5] - -_odeslvrcache, _odeopcache = _odecache - -_odeslvrcache -_odeopcache.const_forms # contains linear systems - -# For each call of evolve need to: -# 1. set _tF to zero -# 2. copy φh to state0 -# 3. reassemble matrices in _odeopcache.const_forms -# 4. solve -# 5. copy stateF to φh - -# User struct: `UnfittedFEEvolution <: LevelSetEvolution` with methods evolve! and reinit! -# Internal structs: `CutFEMEvolve <: Evolver`, `PicardReinit <: Reinitialiser` -# Types that inherit from `UnfittedFEMethod` iterate to re-use cache etc. - +# Based on: +# @article{ +# Burman_Elfverson_Hansbo_Larson_Larsson_2018, +# title={Shape optimization using the cut finite element method}, +# volume={328}, +# ISSN={00457825}, +# DOI={10.1016/j.cma.2017.09.005}, +# journal={Computer Methods in Applied Mechanics and Engineering}, +# author={Burman, Erik and Elfverson, Daniel and Hansbo, Peter and Larson, Mats G. and Larsson, Karl}, +# year={2018}, +# month=jan, +# pages={242–261}, +# language={en} +# } + +using GridapTopOpt + +using Gridap + +using GridapEmbedded +using GridapEmbedded.LevelSetCutters +using Gridap.Geometry, Gridap.FESpaces, Gridap.CellData +import Gridap.Geometry: get_node_coordinates, collect1d + +include("../../differentiable_trians.jl") + +order = 1 +n = 101 + +_model = CartesianDiscreteModel((0,1,0,1),(n,n)) +cd = Gridap.Geometry.get_cartesian_descriptor(_model) +h = maximum(cd.sizes) + +model = simplexify(_model) +Ω = Triangulation(model) +dΩ = Measure(Ω,2*order) +reffe_scalar = ReferenceFE(lagrangian,Float64,order) +V_φ = TestFESpace(model,reffe_scalar) + +φh = interpolate(x->-sqrt((x[1]-0.5)^2+(x[2]-0.5)^2)+0.25,V_φ) # <- Already SDF +# φh = interpolate(x->cos(2π*x[1])*cos(2π*x[2])-0.11,V_φ) # This needs to reinitialised, can use `reinit_based_on_connor_paper.jl` + +geo = DiscreteGeometry(φh,model) +cutgeo = cut(model,geo) +Γ = EmbeddedBoundary(cutgeo) +dΓ = Measure(Γ,2*order) + +F = EmbeddedFacetDiscretization(LevelSetCutter(),model,geo) +FΓ = SkeletonTriangulation(F) +dFΓ = Measure(FΓ,2*order) + +dt = 0.01 +Tn = 5*dt +sysslvr_l = LUSolver() +sysslvr_nl = NLSolver(sysslvr_l, show_trace=true, method=:newton, iterations=10) +ode_solver = RungeKutta(sysslvr_nl, sysslvr_l, dt, :DIRK_CrankNicolson_2_2) + +c = 0.1 + +Ut_φ = TransientTrialFESpace(V_φ,t -> (x->-1)) +n = get_normal_vector(FΓ) +velh = interpolate(x->-1,V_φ) +d1(∇u) = 1 / ( ϵ + norm(∇u) ) +_n(∇u) = ∇u/(10^-20+norm(∇u)) +β = velh*∇(φh)/(10^-20+norm ∘ ∇(φh)) +stiffness(t,u,v) = ∫((β ⋅ ∇(u)) * v)dΩ + ∫(c*h^2*jump(∇(u) ⋅ n)*jump(∇(v) ⋅ n))dFΓ +mass(t, ∂ₜu, v) = ∫(∂ₜu * v)dΩ +forcing(t,v) = ∫(0v)dΩ +op = TransientLinearFEOperator((stiffness,mass),forcing,Ut_φ,V_φ,constant_forms=(true,true)) +uht = solve(ode_solver,op,0.0,Tn,φh) +for (t,uh) in uht + if t ≈ Tn + copyto!(get_free_dof_values(φh),get_free_dof_values(uh)) + end +end + +# You can compare the zero-iso curve generated by the following +φh_expected_contour = interpolate(x->-sqrt((x[1]-0.5)^2+(x[2]-0.5)^2)+0.25+Tn,V_φ) + +writevtk( + Ω,"results/test_evolve", + cellfields=["φh"=>φh,"φh_expected_contour"=>φh_expected_contour,"|∇φh|"=>norm ∘ ∇(φh)] +) + + +sol,state = Base.iterate(uht) +_tF = state[2][1] +_stateF = state[2][2] +_state0 = state[2][3] +_uF = state[2][4] +_odecache = state[2][5] + +_odeslvrcache, _odeopcache = _odecache + +_odeslvrcache +_odeopcache.const_forms # contains linear systems + +# For each call of evolve need to: +# 1. set _tF to zero +# 2. copy φh to state0 +# 3. reassemble matrices in _odeopcache.const_forms +# 4. solve +# 5. copy stateF to φh + +# User struct: `UnfittedFEEvolution <: LevelSetEvolution` with methods evolve! and reinit! +# Internal structs: `CutFEMEvolve <: Evolver`, `PicardReinit <: Reinitialiser` +# Types that inherit from `UnfittedFEMethod` iterate to re-use cache etc. + diff --git a/scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/evolution_test.jl b/scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/evolution_test.jl deleted file mode 100644 index 7077ea4d..00000000 --- a/scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/evolution_test.jl +++ /dev/null @@ -1,88 +0,0 @@ -using GridapTopOpt -using GridapSolvers.NonlinearSolvers - -using Gridap -using Gridap.Adaptivity - -using GridapEmbedded -using GridapEmbedded.LevelSetCutters -using Gridap.Geometry, Gridap.FESpaces, Gridap.CellData -import Gridap.Geometry: get_node_coordinates, collect1d - -include("../../differentiable_trians.jl") - -using GridapTopOpt: LevelSetEvolution, instantiate_caches -using Gridap.Algebra: NonlinearSolver -using Gridap.ODEs: ODESolver -using GridapDistributed: DistributedFESpace -import Gridap: solve! - -include("../../../../src/LevelSetEvolution/UnfittedEvolution/UnfittedEvolution.jl") - -order = 1 -n = 101 -_model = CartesianDiscreteModel((0,1,0,1),(n,n)) -cd = Gridap.Geometry.get_cartesian_descriptor(_model) -base_model = UnstructuredDiscreteModel(_model) -ref_model = refine(base_model, refinement_method = "barycentric") -model = ref_model.model -# model = simplexify(_model) -h = maximum(cd.sizes) - -Ω = Triangulation(model) -dΩ = Measure(Ω,2*order) -reffe_scalar = ReferenceFE(lagrangian,Float64,order) -V_φ = TestFESpace(model,reffe_scalar) - -α = 2.0*h -a_hilb(p,q) = ∫(α^2*∇(p)⋅∇(q) + p*q)dΩ; -vel_ext = VelocityExtension(a_hilb,V_φ,V_φ) - -φh = interpolate(x->-sqrt((x[1]-0.5)^2+(x[2]-0.5)^2)+0.25,V_φ) # <- Already SDF -# φh = interpolate(x->-(x[1]-0.5)^2-(x[2]-0.5)^2+0.25^2,V_φ) -# φh = interpolate(x->cos(4π*x[1])*cos(4π*x[2])-0.11,V_φ) - -function compute_test_velh(φh) - geo = DiscreteGeometry(φh,model) - cutgeo = cut(model,geo) - Ωin = Triangulation(cutgeo,PHYSICAL_IN).a - diff_Ωin = DifferentiableTriangulation(Ωin) - dΩin = Measure(diff_Ωin,2*order) - j_in(φ) = ∫(1)dΩin - dj_in = gradient(j_in,φh) - dj_vec_in = assemble_vector(dj_in,V_φ) - project!(vel_ext,dj_vec_in) - vel = dj_vec_in/sqrt(dj_vec_in'*vel_ext.K*dj_vec_in) - return FEFunction(V_φ,vel) -end - -ls_evo = CutFEMEvolve(model,V_φ,dΩ,h) -ls_reinit = StabilisedReinit(model,V_φ,dΩ,h) -evo = UnfittedFEEvolution(ls_evo,ls_reinit) -# reinit!(evo,φh); - -φ0 = copy(get_free_dof_values(φh)) -φh0 = FEFunction(V_φ,φ0) - -# velh = compute_test_velh(φh) -velh = interpolate(x->-1,V_φ) -evolve!(evo,φh,velh,0.1) - -ls_evo = CutFEMEvolve(model,V_φ,dΩ,h) -ls_reinit = StabilisedReinit(model,V_φ,dΩ,h) -evo = UnfittedFEEvolution(ls_evo,ls_reinit) - -velh = interpolate(x->1,V_φ) -evolve!(evo,φh,velh,0.1) - -Δt = _compute_Δt(h,0.1,get_free_dof_values(velh)) -φh_expected_contour = interpolate(x->-sqrt((x[1]-0.5)^2+(x[2]-0.5)^2)+0.25+evo.evolver.params.max_steps*Δt,V_φ) - -writevtk( - Ω,"results/test_evolve", - cellfields=[ - "φh0"=>φh0,"|∇φh0|"=>norm ∘ ∇(φh0), - "φh_advect"=>φh,"|∇φh_advect|"=>norm ∘ ∇(φh), - # "φh_expected_contour"=>φh_expected_contour - ] -) \ No newline at end of file diff --git a/scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/reinit_test.jl b/scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/reinit_test.jl deleted file mode 100644 index 94a077ec..00000000 --- a/scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/reinit_test.jl +++ /dev/null @@ -1,56 +0,0 @@ -using GridapTopOpt -using GridapSolvers.NonlinearSolvers - -using Gridap -using Gridap.Adaptivity - -using GridapEmbedded -using GridapEmbedded.LevelSetCutters -using Gridap.Geometry, Gridap.FESpaces, Gridap.CellData -import Gridap.Geometry: get_node_coordinates, collect1d - -include("../../differentiable_trians.jl") - -using GridapTopOpt: LevelSetEvolution, instantiate_caches -using Gridap.Algebra: NonlinearSolver -using Gridap.ODEs: ODESolver -using GridapDistributed: DistributedFESpace -import Gridap: solve! - -include("../../../../src/LevelSetEvolution/UnfittedEvolution/UnfittedEvolution.jl") - -order = 1 -n = 101 -_model = CartesianDiscreteModel((0,1,0,1),(n,n)) -cd = Gridap.Geometry.get_cartesian_descriptor(_model) -base_model = UnstructuredDiscreteModel(_model) -ref_model = refine(base_model, refinement_method = "barycentric") -model = ref_model.model -h = maximum(cd.sizes) - -Ω = Triangulation(model) -dΩ = Measure(Ω,2*order) -reffe_scalar = ReferenceFE(lagrangian,Float64,order) -V_φ = TestFESpace(model,reffe_scalar) - -φh = interpolate(x->-(x[1]-0.5)^2-(x[2]-0.5)^2+0.25^2,V_φ) - -ls_evo = CutFEMEvolve(model,V_φ,dΩ,h) -ls_reinit = StabilisedReinit(model,V_φ,dΩ,h) - -φ0 = copy(get_free_dof_values(φh)) -φh0 = FEFunction(V_φ,φ0) - -evo = UnfittedFEEvolution(ls_evo,ls_reinit) -reinit!(evo,φh); - -L2error(u) = ∫(u*u)dΩ -@assert sum(L2error((norm ∘ ∇(φh_reinit)))) - 1 < 10^-4 - -writevtk( - Ω,"results/test_evolve", - cellfields=[ - "φh0"=>φh0,"|∇φh0|"=>norm ∘ ∇(φh0), - "φh_reinit"=>φh_reinit,"|∇φh_reinit|"=>norm ∘ ∇(φh_reinit) - ] -) \ No newline at end of file diff --git a/src/GridapTopOpt.jl b/src/GridapTopOpt.jl index aa25bcac..20d72590 100644 --- a/src/GridapTopOpt.jl +++ b/src/GridapTopOpt.jl @@ -16,6 +16,7 @@ using Gridap.ReferenceFEs, Gridap.FESpaces, Gridap.MultiField using Gridap.Geometry: get_faces, num_nodes using Gridap.FESpaces: get_assembly_strategy +using Gridap.ODEs: ODESolver using Gridap: writevtk using GridapDistributed @@ -38,6 +39,7 @@ using GridapEmbedded.Interfaces: SubFacetData, SubCellTriangulation, SubFacetTri using JLD2: save_object, load_object, jldsave import Base: + +import Gridap: solve! include("GridapExtensions.jl") @@ -65,6 +67,11 @@ export isotropic_elast_tensor include("LevelSetEvolution/LevelSetEvolution.jl") export HamiltonJacobiEvolution export FirstOrderStencil +export UnfittedFEEvolution +export CutFEMEvolve +export StabilisedReinit +export ArtificialViscosity +export InteriorPenalty export evolve! export reinit! diff --git a/src/LevelSetEvolution/LevelSetEvolution.jl b/src/LevelSetEvolution/LevelSetEvolution.jl index 6f9c4f32..c1081a08 100644 --- a/src/LevelSetEvolution/LevelSetEvolution.jl +++ b/src/LevelSetEvolution/LevelSetEvolution.jl @@ -36,4 +36,5 @@ function get_dof_Δ(::LevelSetEvolution) end include("Stencil.jl") -include("HamiltonJacobiEvolution.jl") \ No newline at end of file +include("HamiltonJacobiEvolution.jl") +include("UnfittedEvolution/UnfittedEvolution.jl") \ No newline at end of file diff --git a/src/LevelSetEvolution/UnfittedEvolution/CutFEMEvolve.jl b/src/LevelSetEvolution/UnfittedEvolution/CutFEMEvolve.jl index 626329a9..c50c3729 100644 --- a/src/LevelSetEvolution/UnfittedEvolution/CutFEMEvolve.jl +++ b/src/LevelSetEvolution/UnfittedEvolution/CutFEMEvolve.jl @@ -46,7 +46,10 @@ get_measure(s::CutFEMEvolve) = s.dΩ get_params(s::CutFEMEvolve) = s.params get_cache(s::CutFEMEvolve) = s.cache -function get_stiffness_matrix_map(φh,velh,dΩ,model,order,Γg,h) +function get_transient_operator(φh,velh,s::CutFEMEvolve) + model, V_φ, dΩ, assembler, params = s.model, s.space, s.dΩ, s.assembler, s.params + Γg, h, order = params.Γg, params.h, params.order + geo = DiscreteGeometry(φh,model) cutgeo = cut(model,geo) Ω_act = Triangulation(cutgeo,ACTIVE) @@ -56,22 +59,33 @@ function get_stiffness_matrix_map(φh,velh,dΩ,model,order,Γg,h) ϵ = 1e-20 β = velh*∇(φh)/(ϵ + norm ∘ ∇(φh)) - return (t,u,v) -> ∫((β ⋅ ∇(u)) * v)dΩ + ∫(Γg*h^2*jump(∇(u) ⋅ n)*jump(∇(v) ⋅ n))dF_act + stiffness(t,u,v) = ∫((β ⋅ ∇(u)) * v)dΩ + ∫(Γg*h^2*jump(∇(u) ⋅ n)*jump(∇(v) ⋅ n))dF_act + 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),assembler) + return ode_op end -function solve!(s::CutFEMEvolve,φh,velh,γ,cache::Nothing) - ode_solver, model, V_φ, dΩ, assembler = s.ode_solver, s.model, s.space, s.dΩ, s.assembler - Γg, h, max_steps, order = s.params +function update_reuse!(state,reuse_new;zero_tF=false) + U, (tF, stateF, state0, uF, odecache) = state + odeslvrcache, odeopcache = odecache + _, ui_pre, slopes, J, r, sysslvrcaches = odeslvrcache - # Weak form - stiffness = get_stiffness_matrix_map(φh,velh,dΩ,model,order,Γg,h) - mass(t, ∂ₜu, v) = ∫(∂ₜu * v)dΩ - forcing(t,v) = ∫(0v)dΩ + odeslvrcache_new = (reuse_new, ui_pre, slopes, J, r, sysslvrcaches) + odecache_new = odeslvrcache_new, odeopcache + _tF = zero_tF ? 0.0 : tF + return U, (_tF, stateF, state0, uF, odecache_new) +end + +function solve!(s::CutFEMEvolve,φh,velh,γ,cache::Nothing) + ode_solver = get_ode_solver(s) + params = get_params(s) + h, max_steps = params.h, params.max_steps # Setup FE operator and solver - Ut_φ = TransientTrialFESpace(V_φ) - ode_op = TransientLinearFEOperator((stiffness,mass),forcing,Ut_φ,V_φ; - constant_forms=(true,true),assembler) + ode_op = get_transient_operator(φh,velh,s) dt = _compute_Δt(h,γ,get_free_dof_values(velh)) ode_solver.dt = dt ode_sol = solve(ode_solver,ode_op,0.0,dt*max_steps,φh) @@ -79,48 +93,56 @@ function solve!(s::CutFEMEvolve,φh,velh,γ,cache::Nothing) # March march = Base.iterate(ode_sol) data, state = march - while march !== nothing - data, state = march - march = Base.iterate(ode_sol,state) + state_new = update_reuse!(state,true) + march_new = data, state_new + while march_new !== nothing + data, state_new = march_new + march_new = Base.iterate(ode_sol,state_new) end # Update φh and cache _, φhF = data copy!(get_free_dof_values(φh),get_free_dof_values(φhF)) - s.cache = (ode_op, state) + s.cache = state_new return φh end function solve!(s::CutFEMEvolve,φh,velh,γ,cache) - ode_solver, model, dΩ, V_φ, assembler, = s.ode_solver, s.model, s.dΩ, s.space, s.assembler - Γg, h, max_steps, order = s.params - ode_op, state0 = cache - - U, (tF, stateF, state0, uF, odecache) = state0 - odeslvrcache, odeopcache = odecache - K, M = odeopcache.const_forms - - # Update odeopcache and solver - stiffness = get_stiffness_matrix_map(φh,velh,dΩ,model,order,Γg,h) - assemble_matrix!((u,v)->stiffness(0.0,u,v),K,assembler,V_φ,V_φ) - state = U, (0.0, stateF, state0, uF, odecache) + ode_solver = get_ode_solver(s) + params = get_params(s) + h, max_steps = params.h, params.max_steps + + ## Update state + # `get_transient_operator` re-creates the entire TransientLinearFEOperator wrapper. + # We do this so that the first iterate of ODESolution always recomputes the + # stiffness matrix and associated the Jacboian, numerical setups, etc via + # `constant_forms = (false,true)`. + ode_op = get_transient_operator(φh,velh,s) + # Between the first iterate and subsequent iterates we use the function + # `update_reuse!` to update the iterator state so that we re-use + # the stiffness matrix, etc. The Optional argument `zero_tF` indicates + # whether we are solving a new ODE with the same functional form but + # updated coefficients in the weak form. If so, we want to re-use the cache. + state_inter = update_reuse!(cache,false;zero_tF=true) dt = _compute_Δt(h,γ,get_free_dof_values(velh)) ode_solver.dt = dt - # March + ## March ode_sol = solve(ode_solver,ode_op,0.0,dt*max_steps,φh) - march = Base.iterate(ode_sol,state) + march = Base.iterate(ode_sol,state_inter) # First step includes stiffness matrix update data, state = march - while march !== nothing - data, state = march - march = Base.iterate(ode_sol,state) + state_updated = update_reuse!(state,true) # Fix the stiffness matrix for remaining march + march_updated = data, state_updated + while march_updated !== nothing + data, state_updated = march_updated + march_updated = Base.iterate(ode_sol,state_updated) end - # Update φh and cache + ## Update φh and cache _, φhF = data copy!(get_free_dof_values(φh),get_free_dof_values(φhF)) - s.cache = (ode_op, state) + s.cache = state_updated return φh end diff --git a/src/LevelSetEvolution/UnfittedEvolution/StabilisedReinit.jl b/src/LevelSetEvolution/UnfittedEvolution/StabilisedReinit.jl index 2c89de8e..fb1e4aba 100644 --- a/src/LevelSetEvolution/UnfittedEvolution/StabilisedReinit.jl +++ b/src/LevelSetEvolution/UnfittedEvolution/StabilisedReinit.jl @@ -9,7 +9,15 @@ Interior jump penalty approach (`InteriorPenalty`) adapted from that work and re 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 +- `model::B`: FE model +- `space::C`: FE space for level-set function +- `dΩ::D`: Background mesh measure +- `assembler::Assembler`: Assembler for LS FE space +- `params::E`: Tuple of Nitsche parameter `γd`, mesh size `h`, + and FE space `order` +- `cache`: Cache for evolver, initially `nothing`. """ mutable struct StabilisedReinit{A,B,C,D,E} <: Reinitialiser diff --git a/test/seq/UnfittedEvolutionTests/ArtificialViscosityStabilisedReinitTest.jl b/test/seq/UnfittedEvolutionTests/ArtificialViscosityStabilisedReinitTest.jl new file mode 100644 index 00000000..34bcb62f --- /dev/null +++ b/test/seq/UnfittedEvolutionTests/ArtificialViscosityStabilisedReinitTest.jl @@ -0,0 +1,54 @@ +module ArtificialViscosityStabilisedReinitTest +using Test + +using GridapTopOpt +using Gridap, Gridap.Geometry, Gridap.Adaptivity +using GridapEmbedded, GridapEmbedded.LevelSetCutters + +order = 1 +n = 201 +_model = CartesianDiscreteModel((0,1,0,1),(n,n)) +cd = Gridap.Geometry.get_cartesian_descriptor(_model) +base_model = UnstructuredDiscreteModel(_model) +ref_model = refine(base_model, refinement_method = "barycentric") +model = ref_model.model +h = maximum(cd.sizes) + +Ω = Triangulation(model) +dΩ = Measure(Ω,2*order) +reffe_scalar = ReferenceFE(lagrangian,Float64,order) +V_φ = TestFESpace(model,reffe_scalar) + +φh = interpolate(x->-(x[1]-0.5)^2-(x[2]-0.5)^2+0.25^2,V_φ) +φh0 = interpolate(x->-sqrt((x[1]-0.5)^2+(x[2]-0.5)^2)+0.25,V_φ) + +ls_evo = CutFEMEvolve(model,V_φ,dΩ,h) +ls_reinit = StabilisedReinit(model,V_φ,dΩ,h;stabilisation_method=ArtificialViscosity(1.5h)) +evo = UnfittedFEEvolution(ls_evo,ls_reinit) +reinit!(evo,φh); + +L2error(u) = sqrt(sum(∫(u ⋅ u)dΩ)) +# Check |∇(φh)| +@test abs(L2error(norm ∘ ∇(φh))-1) < 1e-4 + +# Check φh error +@test L2error(φh-φh0) < 1e-4 + +# Check facet coords +geo = DiscreteGeometry(φh,model) +geo0 = DiscreteGeometry(φh0,model) +cutgeo = cut(model,geo) +cutgeo0 = cut(model,geo0) +Γ = EmbeddedBoundary(cutgeo) +Γ0 = EmbeddedBoundary(cutgeo0) +@test norm(Γ.subfacets.point_to_coords - Γ0.subfacets.point_to_coords,Inf) < 1e-4 + +# writevtk( +# Ω,"results/test_evolve", +# cellfields=[ +# "φh0"=>φh0,"|∇φh0|"=>norm ∘ ∇(φh0), +# "φh_reinit"=>φh_reinit,"|∇φh_reinit|"=>norm ∘ ∇(φh_reinit) +# ] +# ) + +end \ No newline at end of file diff --git a/test/seq/UnfittedEvolutionTests/CutFEMEvolveTest.jl b/test/seq/UnfittedEvolutionTests/CutFEMEvolveTest.jl new file mode 100644 index 00000000..a73ee3d6 --- /dev/null +++ b/test/seq/UnfittedEvolutionTests/CutFEMEvolveTest.jl @@ -0,0 +1,54 @@ +module CutFEMEvolveTest +using Test + +using GridapTopOpt +using Gridap, Gridap.Geometry, Gridap.Adaptivity + +order = 1 +n = 201 +_model = CartesianDiscreteModel((0,1,0,1),(n,n)) +cd = Gridap.Geometry.get_cartesian_descriptor(_model) +base_model = UnstructuredDiscreteModel(_model) +ref_model = refine(base_model, refinement_method = "barycentric") +model = ref_model.model +h = maximum(cd.sizes) + +Ω = Triangulation(model) +dΩ = Measure(Ω,2*order) +reffe_scalar = ReferenceFE(lagrangian,Float64,order) +V_φ = TestFESpace(model,reffe_scalar) + +φh = interpolate(x->-sqrt((x[1]-0.5)^2+(x[2]-0.5)^2)+0.25,V_φ) + +ls_evo = CutFEMEvolve(model,V_φ,dΩ,h) +ls_reinit = StabilisedReinit(model,V_φ,dΩ,h) +evo = UnfittedFEEvolution(ls_evo,ls_reinit) + +φ0 = copy(get_free_dof_values(φh)) +φh0 = FEFunction(V_φ,φ0) + +velh = interpolate(x->-1,V_φ) +evolve!(evo,φh,velh,0.1) + +Δt = GridapTopOpt._compute_Δt(h,0.1,get_free_dof_values(velh)) +φh_expected_lsf = interpolate(x->-sqrt((x[1]-0.5)^2+(x[2]-0.5)^2)+0.25+evo.evolver.params.max_steps*Δt,V_φ) + +# Test advected LSF mataches expected LSF +L2error(u) = sqrt(sum(∫(u ⋅ u)dΩ)) +@test L2error(φh_expected_lsf-φh) < 1e-4 + +# # Test advected LSF mataches original LSF when going backwards +velh = interpolate(x->1,V_φ) +evolve!(evo,φh,velh,0.1) +@test L2error(φh0-φh) < 1e-5 + +# writevtk( +# Ω,"results/test_evolve", +# cellfields=[ +# "φh0"=>φh0,"|∇φh0|"=>norm ∘ ∇(φh0), +# "φh_advect"=>φh,"|∇φh_advect|"=>norm ∘ ∇(φh), +# "φh_expected_lsf"=>φh_expected_lsf +# ] +# ) + +end \ No newline at end of file diff --git a/test/seq/UnfittedEvolutionTests/InteriorPenaltyStabilisedReinitTest.jl b/test/seq/UnfittedEvolutionTests/InteriorPenaltyStabilisedReinitTest.jl new file mode 100644 index 00000000..c914fc3e --- /dev/null +++ b/test/seq/UnfittedEvolutionTests/InteriorPenaltyStabilisedReinitTest.jl @@ -0,0 +1,55 @@ +module InteriorPenaltyStabilisedReinitTest +using Test + +using GridapTopOpt +using Gridap, Gridap.Geometry, Gridap.Adaptivity +using GridapEmbedded, GridapEmbedded.LevelSetCutters + +order = 1 +n = 201 +_model = CartesianDiscreteModel((0,1,0,1),(n,n)) +cd = Gridap.Geometry.get_cartesian_descriptor(_model) +base_model = UnstructuredDiscreteModel(_model) +ref_model = refine(base_model, refinement_method = "barycentric") +model = ref_model.model +h = maximum(cd.sizes) + +Ω = Triangulation(model) +dΩ = Measure(Ω,2*order) +reffe_scalar = ReferenceFE(lagrangian,Float64,order) +V_φ = TestFESpace(model,reffe_scalar) + +φh = interpolate(x->-(x[1]-0.5)^2-(x[2]-0.5)^2+0.25^2,V_φ) +φh0 = interpolate(x->-sqrt((x[1]-0.5)^2+(x[2]-0.5)^2)+0.25,V_φ) + +ls_evo = CutFEMEvolve(model,V_φ,dΩ,h) +ls_reinit = StabilisedReinit(model,V_φ,dΩ,h;stabilisation_method=InteriorPenalty(V_φ)) +evo = UnfittedFEEvolution(ls_evo,ls_reinit) +reinit!(evo,φh); + +L2error(u) = sqrt(sum(∫(u ⋅ u)dΩ)) +# Check |∇(φh)| +abs(L2error(norm ∘ ∇(φh))-1) < 1e-4 + +# Check φh error +@test L2error(φh-φh0) < 1e-4 + +# Check facet coords +geo = DiscreteGeometry(φh,model) +geo0 = DiscreteGeometry(φh0,model) +cutgeo = cut(model,geo) +cutgeo0 = cut(model,geo0) +Γ = EmbeddedBoundary(cutgeo) +Γ0 = EmbeddedBoundary(cutgeo0) +@test norm(Γ.subfacets.point_to_coords - Γ0.subfacets.point_to_coords,Inf) < 1e-4 + + +# writevtk( +# Ω,"results/test_evolve", +# cellfields=[ +# "φh0"=>φh0,"|∇φh0|"=>norm ∘ ∇(φh0), +# "φh_reinit"=>φh_reinit,"|∇φh_reinit|"=>norm ∘ ∇(φh_reinit) +# ] +# ) + +end \ No newline at end of file diff --git a/test/seq/UnfittedEvolutionTests/runtests.jl b/test/seq/UnfittedEvolutionTests/runtests.jl new file mode 100644 index 00000000..6e264cb8 --- /dev/null +++ b/test/seq/UnfittedEvolutionTests/runtests.jl @@ -0,0 +1,9 @@ +module UnfittedEvolutionTests + +using Test + +@testset "Reinitialisation (artificial viscosity)" begin include("ArtificialViscosityStabilisedReinitTest.jl") end +@testset "Reinitialisation (interior penalty)" begin include("InteriorPenaltyStabilisedReinitTest.jl") end +@testset "CutFEM evolution" begin include("CutFEMEvolveTest.jl") end + +end # module \ No newline at end of file diff --git a/test/seq/runtests.jl b/test/seq/runtests.jl index cea7a334..ce314c1e 100644 --- a/test/seq/runtests.jl +++ b/test/seq/runtests.jl @@ -14,5 +14,6 @@ using Test include("EmbeddedDifferentiationTests.jl") include("EmbeddedCollectionsTests.jl") end +@time @testset "UnfittedEvolution" begin include("UnfittedEvolutionTests/runtests.jl") end end # module \ No newline at end of file