diff --git a/scripts/Embedded/Bugs/odes_reassemble_stiffness.jl b/scripts/Embedded/Bugs/odes_reassemble_stiffness.jl new file mode 100644 index 00000000..1b320632 --- /dev/null +++ b/scripts/Embedded/Bugs/odes_reassemble_stiffness.jl @@ -0,0 +1,250 @@ +using Test + +using GridapTopOpt +using Gridap, Gridap.Geometry, Gridap.Adaptivity, Gridap.ODEs +using GridapEmbedded + +function update_reuse!(state,reuse_new,op;zero_tF=false) + U, (tF, stateF, state0, uF, odecache) = state + odeslvrcache, odeopcache = odecache + _, ui_pre, slopes, J, r, sysslvrcaches = odeslvrcache + + data = allocate_odecache(ode_solver,ODEOpFromTFEOp(op),tF,stateF) + odeslvrcache_new = (reuse_new, ui_pre, slopes, J, r, data[1][end]) + odecache_new = odeslvrcache_new, odeopcache + _tF = zero_tF ? 0.0 : tF + return U, (_tF, stateF, state0, uF, odecache_new) +end + +order = 1 +n = 50 +_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) +steps = 10 + +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_φ) +velh = interpolate(x->-1,V_φ) + +## Trians + +Ω = Triangulation(model) +dΩ = Measure(Ω,2*order) +Γg = SkeletonTriangulation(Ω) +dΓg = Measure(Γg,2order) +n_Γg = get_normal_vector(Γg) + +## ODE Solver +ode_ls = LUSolver() +ode_nl = NLSolver(ode_ls, show_trace=false, method=:newton, iterations=10) +ode_solver = GridapTopOpt.MutableRungeKutta(ode_nl, ode_ls, 0.1*h, :DIRK_CrankNicolson_2_2) + +## ODE Op +v_norm = maximum(abs,get_free_dof_values(velh)) +β(vh,∇φ) = vh/(1e-20 + v_norm) * ∇φ/(1e-20 + norm(∇φ)) +stiffness(t,u,v) = ∫(((β ∘ (velh,∇(φh))) ⋅ ∇(u)) * v)dΩ + ∫(0.1*h^2*jump(∇(u) ⋅ n_Γg)*jump(∇(v) ⋅ n_Γg))dΓg +mass(t, ∂ₜu, v) = ∫(∂ₜu * v)dΩ +forcing(t,v) = ∫(0v)dΩ + ∫(0*jump(∇(v) ⋅ n_Γg))dΓg +Ut_φ = TransientTrialFESpace(V_φ) + +# Both forms constant +op = TransientLinearFEOperator((stiffness,mass),forcing,Ut_φ,V_φ; + constant_forms=(true,true)) +ode_sol = ODEs.solve(ode_solver,op,0.0,ode_solver.dt*steps,φh) +data0, state0 = Base.iterate(ode_sol) +data1, state1 = Base.iterate(ode_sol,state0) + +# # Update φh and velh +# φh_new = FEFunction(V_φ,copy(get_free_dof_values(data1[2]))) +# velh_new = interpolate(x->-2,V_φ) +# v_norm_new = maximum(abs,get_free_dof_values(velh_new)) +# β_new(vh,∇φ) = vh/(1e-20 + v_norm_new) * ∇φ/(1e-20 + norm(∇φ)) +# stiffness_new(t,u,v) = ∫(((β_new ∘ (velh,∇(φh))) ⋅ ∇(u)) * v)dΩ + ∫(0.1*h^2*jump(∇(u) ⋅ n_Γg)*jump(∇(v) ⋅ n_Γg))dΓg +# op = TransientLinearFEOperator((stiffness_new,mass),forcing,Ut_φ,V_φ; +# constant_forms=(true,true)) +# ode_sol_new = ODEs.solve(ode_solver,op,0.0,ode_solver.dt*steps,φh_new) +# data2, state2 = Base.iterate(ode_sol) +# data3, state3 = Base.iterate(ode_sol,state2) + +# first form non-constant +nc_op = TransientLinearFEOperator((stiffness,mass),forcing,Ut_φ,V_φ; + constant_forms=(false,true)) +nc_ode_sol = ODEs.solve(ode_solver,nc_op,0.0,ode_solver.dt*steps,φh) +nc_data0, nc_state0 = Base.iterate(nc_ode_sol) +nc_data1, nc_state1 = Base.iterate(nc_ode_sol,nc_state0) + +data1[1] == nc_data1[1] +norm(get_free_dof_values(data1[2]) - get_free_dof_values(nc_data1[2]),Inf) + +# first form non-constant then switch to constant +s_op = TransientLinearFEOperator((stiffness,mass),forcing,Ut_φ,V_φ; + constant_forms=(false,true)) +s_ode_sol = ODEs.solve(ode_solver,s_op,0.0,ode_solver.dt*steps,φh) +s_data0, s_state0 = Base.iterate(s_ode_sol) +s_state0_new = update_reuse!(s_state0,true) +s_data1, s_state1 = Base.iterate(s_ode_sol,s_state0_new) + +data1[1] == s_data1[1] +norm(get_free_dof_values(data1[2]) - get_free_dof_values(s_data1[2]),Inf) + +######################################################## +using Gridap.FESpaces, Gridap.Polynomials +using LinearAlgebra +# TrianientLinearFEOpFromWeakForm2 +struct TransientLinearFEOpFromWeakForm2 <: TransientFEOperator{LinearODE} + forms::Tuple{Vararg{Function}} + res::Function + jacs::Tuple{Vararg{Function}} + constant_forms::BitVector + assembler::Assembler + trial::FESpace + test::FESpace + order::Integer +end + +# Constructor with manual jacobians +function TransientLinearFEOperator2( + forms::Tuple{Vararg{Function}}, res::Function, jacs::Tuple{Vararg{Function}}, + trial, test; + constant_forms::BitVector=falses(length(forms)), + assembler=SparseMatrixAssembler(trial, test) +) + order = length(jacs) - 1 + TransientLinearFEOpFromWeakForm2( + forms, res, jacs, constant_forms, + assembler, trial, test, order + ) +end + +# No constructor with flat arguments: would clash with the constructors +# below with flat forms and automatic jacobians, which are more useful + +# Constructor with automatic jacobians +function TransientLinearFEOperator2( + forms::Tuple{Vararg{Function}}, res::Function, + trial, test; + constant_forms::BitVector=falses(length(forms)), + assembler=SparseMatrixAssembler(trial, test) +) + # When the operator is linear, the jacobians are the forms themselves + order = length(forms) - 1 + jacs = ntuple(k -> ((t, u, duk, v) -> forms[k](t, duk, v)), order + 1) + + TransientLinearFEOperator2( + forms, res, jacs, trial, test; + constant_forms, assembler + ) +end + +# Constructor with flat forms and automatic jacobians (orders 0, 1, 2) +function TransientLinearFEOperator2( + mass::Function, res::Function, + trial, test; + constant_forms::BitVector=falses(1), + assembler=SparseMatrixAssembler(trial, test) +) + TransientLinearFEOperator2( + (mass,), res, trial, test; + constant_forms, assembler + ) +end + +function TransientLinearFEOperator2( + stiffness::Function, mass::Function, res::Function, + trial, test; + constant_forms::BitVector=falses(2), + assembler=SparseMatrixAssembler(trial, test) +) + TransientLinearFEOperator2( + (stiffness, mass), res, trial, test; + constant_forms, assembler + ) +end + +function TransientLinearFEOperator2( + stiffness::Function, damping::Function, mass::Function, res::Function, + trial, test; + constant_forms::BitVector=falses(3), + assembler=SparseMatrixAssembler(trial, test) +) + TransientLinearFEOpFromWeakForm2( + (stiffness, damping, mass), res, trial, test; + constant_forms, assembler + ) +end + +# TransientFEOperator interface +FESpaces.get_test(tfeop::TransientLinearFEOpFromWeakForm2) = tfeop.test + +FESpaces.get_trial(tfeop::TransientLinearFEOpFromWeakForm2) = tfeop.trial + +Polynomials.get_order(tfeop::TransientLinearFEOpFromWeakForm2) = tfeop.order + +ODEs.get_res(tfeop::TransientLinearFEOpFromWeakForm2) = (t, u, v) -> tfeop.res(t, v) + +ODEs.get_jacs(tfeop::TransientLinearFEOpFromWeakForm2) = tfeop.jacs + +ODEs.get_forms(tfeop::TransientLinearFEOpFromWeakForm2) = tfeop.forms + +function ODEs.is_form_constant(tfeop::TransientLinearFEOpFromWeakForm2, k::Integer) + tfeop.constant_forms[k+1] +end + +ODEs.get_assembler(tfeop::TransientLinearFEOpFromWeakForm2) = tfeop.assembler + +####### +ss_op = TransientLinearFEOperator2((stiffness,mass),forcing,Ut_φ,V_φ; + constant_forms=BitVector((false,true))) +ss_ode_sol = ODEs.solve(ode_solver,ss_op,0.0,ode_solver.dt*steps,φh) +ss_data0, ss_state0 = Base.iterate(ss_ode_sol) + +_t = ss_data0[1] +stiff_const_form = copy(ss_state0[2][5][1][4]) # <- cache this guy + +#... + +us = (get_free_dof_values(ss_data0[2]),get_free_dof_values(ss_data0[2])) +order = get_order(ss_op) +Ut = get_trial(ss_op) +# U = allocate_space(Ut) +# Uts = (Ut,) +# Us = (U,) +# for k in 1:order +# Uts = (Uts..., ∂t(Uts[k])) +# Us = (Us..., allocate_space(Uts[k+1])) +# end + +_odeopcache = ss_state0[end][end][end] +Us = () +for k in 0:get_order(ss_op) + Us = (Us..., evaluate!(_odeopcache.Us[k+1], _odeopcache.Uts[k+1], _t)) +end + +uh = ODEs._make_uh_from_us(ss_op, us, Us) +du = get_trial_fe_basis(Ut_φ) +V = get_test(ss_op) +v = get_fe_basis(V) + +jacs = ODEs.get_jacs(ss_op) +jac = jacs[1] +dc = jac(_t, uh, du, v) +matdata = collect_cell_matrix(Ut, V, dc) +LinearAlgebra.fillstored!(stiff_const_form, zero(eltype(stiff_const_form))) +assemble_matrix_add!(stiff_const_form, get_assembler(ss_op), matdata) + +## Update +new_const_forms = (stiff_const_form,ss_state0[end][end][end].const_forms[2]) +ss_state0[end][end][end].const_forms = new_const_forms + +ss_state0_new = update_reuse!(ss_state0,true) +ss_op.constant_forms[1] = true +ss_data1, ss_state1 = Base.iterate(ss_ode_sol,ss_state0_new) + +data1[1] == ss_data1[1] +norm(get_free_dof_values(data1[2]) - get_free_dof_values(ss_data1[2]),Inf) + diff --git a/scripts/Embedded/Examples/FCM_2d_compliance_L.jl b/scripts/Embedded/Examples/FCM_2d_compliance_L.jl index a5e4f105..1e9c9b47 100644 --- a/scripts/Embedded/Examples/FCM_2d_compliance_L.jl +++ b/scripts/Embedded/Examples/FCM_2d_compliance_L.jl @@ -4,15 +4,15 @@ using GridapEmbedded, GridapEmbedded.LevelSetCutters using GridapTopOpt: StateParamIntegrandWithMeasure -path="./results/FCM_elastic_compliance_LShape_ALM_diffusion/" +path="./results/FCM_elastic_compliance_LShape_ALM_n100/" rm(path,force=true,recursive=true) mkpath(path) n = 100 order = 1 -γ = 0.1 -max_steps = floor(Int,order*n/5) +γ = 0.2 +max_steps = floor(Int,order*n/10) vf = 0.4 -α_coeff = 3max_steps*γ +α_coeff = max_steps*γ iter_mod = 1 _model = CartesianDiscreteModel((0,1,0,1),(n,n)) @@ -49,7 +49,7 @@ U_reg = TrialFESpace(V_reg,0) V_φ = TestFESpace(model,reffe_scalar) ## Levet-set function -φh = interpolate(x->-cos(8π*x[1])*cos(8π*x[2])-0.2,V_φ) +φh = interpolate(x->-cos(6π*x[1])*cos(8π*x[2])-0.2,V_φ) Ωs = EmbeddedCollection(model,φh) do cutgeo,_ Ωin = DifferentiableTriangulation(Triangulation(cutgeo,PHYSICAL_IN),V_φ) Ωout = DifferentiableTriangulation(Triangulation(cutgeo,PHYSICAL_OUT),V_φ) @@ -84,8 +84,7 @@ g = VectorValue(0,-0.1) const ϵ = (λ + μ)*1e-3 a(u,v,φ) = ∫(ε(v) ⊙ (σ ∘ ε(u)))Ωs.dΩin + - # ∫(ϵ*(ε(v) ⊙ (σ ∘ ε(u))))Ωs.dΩout # Ersatz - ∫(ϵ * ∇(v) ⊙ ∇(u))Ωs.dΩout # Pure diffusion + ∫(ϵ*(ε(v) ⊙ (σ ∘ ε(u))))Ωs.dΩout l(v,φ) = ∫(v⋅g)dΓ_N ## Optimisation functionals @@ -101,13 +100,14 @@ state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh) pcfs = PDEConstrainedFunctionals(J,[Vol],state_map;analytic_dC=(dVol,)) ## Evolution Method -evo = CutFEMEvolve(V_φ,Ωs,dΩ,h;max_steps,γg=0.5) -reinit = StabilisedReinit(V_φ,Ωs,dΩ,h;stabilisation_method=ArtificialViscosity(3.0)) +evo = CutFEMEvolve(V_φ,Ωs,dΩ,h;max_steps,γg=0.1) +reinit1 = StabilisedReinit(V_φ,Ωs,dΩ,h;stabilisation_method=ArtificialViscosity(3.0)) +reinit2 = StabilisedReinit(V_φ,Ωs,dΩ,h;stabilisation_method=GridapTopOpt.InteriorPenalty(V_φ,γg=2.0)) +reinit = GridapTopOpt.MultiStageStabilisedReinit([reinit1,reinit2]) ls_evo = UnfittedFEEvolution(evo,reinit) -reinit!(ls_evo,φh) ## Hilbertian extension-regularisation problems -α = α_coeff*(h_refine/order)^2 +α = (α_coeff*h_refine/order)^2 a_hilb(p,q) =∫(α*∇(p)⋅∇(q) + p*q)dΩ; vel_ext = VelocityExtension(a_hilb,U_reg,V_reg) diff --git a/scripts/Embedded/Examples/FCM_2d_thermal.jl b/scripts/Embedded/Examples/FCM_2d_thermal.jl index eff48bee..16d8613c 100644 --- a/scripts/Embedded/Examples/FCM_2d_thermal.jl +++ b/scripts/Embedded/Examples/FCM_2d_thermal.jl @@ -9,10 +9,10 @@ rm(path,force=true,recursive=true) mkpath(path) n = 50 order = 1 -γ = 0.1 -max_steps = floor(Int,order*n/5) +γ = 0.2 +max_steps = floor(Int,order*n/10) vf = 0.4 -α_coeff = 4max_steps*γ +α_coeff = max_steps*γ iter_mod = 1 _model = CartesianDiscreteModel((0,1,0,1),(n,n)) @@ -79,10 +79,11 @@ state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh) pcfs = PDEConstrainedFunctionals(J,[Vol],state_map;analytic_dC=(dVol,)) ## Evolution Method -evo = CutFEMEvolve(V_φ,Ωs,dΩ,h;max_steps) -reinit = StabilisedReinit(V_φ,Ωs,dΩ,h;stabilisation_method=ArtificialViscosity(3.0)) +evo = CutFEMEvolve(V_φ,Ωs,dΩ,h;max_steps,γg=0.1) +reinit1 = StabilisedReinit(V_φ,Ωs,dΩ,h;stabilisation_method=ArtificialViscosity(3.0)) +reinit2 = StabilisedReinit(V_φ,Ωs,dΩ,h;stabilisation_method=GridapTopOpt.InteriorPenalty(V_φ,γg=2.0)) +reinit = GridapTopOpt.MultiStageStabilisedReinit([reinit1,reinit2]) ls_evo = UnfittedFEEvolution(evo,reinit) -reinit!(ls_evo,φh) ## Hilbertian extension-regularisation problems α = α_coeff*(h_refine/order)^2 diff --git a/scripts/Embedded/Examples/fsi/gmsh/GMSH_TO-6-Brinkmann_stokes_P2-P1_Ersatz_elast_fsi.jl b/scripts/Embedded/Examples/fsi/gmsh/GMSH_TO-6-Brinkmann_stokes_P2-P1_Ersatz_elast_fsi.jl index f2d4e64e..614f04eb 100644 --- a/scripts/Embedded/Examples/fsi/gmsh/GMSH_TO-6-Brinkmann_stokes_P2-P1_Ersatz_elast_fsi.jl +++ b/scripts/Embedded/Examples/fsi/gmsh/GMSH_TO-6-Brinkmann_stokes_P2-P1_Ersatz_elast_fsi.jl @@ -5,13 +5,13 @@ using GridapTopOpt ############# -path = "./results/GMSH-TO-6-Brinkmann_stokes_P2-P1_Ersatz_elast_fsi_multistage_reinit/results/" +path = "./results/GMSH-TO-6-Brinkmann_stokes_P2-P1_Ersatz_elast_fsi/results/" mkpath(path) γ_evo = 0.1 max_steps = 20 vf = 0.025 -α_coeff = 2 +α_coeff = 1 iter_mod = 1 D = 2 @@ -148,16 +148,14 @@ l_coupled((v,q,s),φ) = ∫(0.0q)dΩ_act J_pres((u,p,d),φ) = ∫(p)dΓf_D - ∫(p)dΓf_N J_comp((u,p,d),φ) = ∫(ε(d) ⊙ (σ ∘ ε(d)))Ω.dΩs Vol((u,p,d),φ) = ∫(vol_D)Ω.dΩs - ∫(vf/vol_D)dΩ_act +dVol(q,(u,p,d),φ) = ∫(-1/vol_D*q/(norm ∘ (∇(φ))))Ω.dΓ ## Setup solver and FE operators state_map = AffineFEStateMap(a_coupled,l_coupled,X,Y,V_φ,U_reg,φh) -pcfs = PDEConstrainedFunctionals(J_comp,[Vol],state_map) +pcfs = PDEConstrainedFunctionals(J_comp,[Vol],state_map;analytic_dC=(dVol,)) ## Evolution Method -# evo = CutFEMEvolve(V_φ,Ω,dΩ_act,hmin;max_steps,γg=0.075) -# reinit = StabilisedReinit(V_φ,Ω,dΩ_act,hmin;stabilisation_method=ArtificialViscosity(3.0)) -# ls_evo = UnfittedFEEvolution(evo,reinit) -evo = CutFEMEvolve(V_φ,Ω,dΩ_act,hₕ;max_steps,γg=0.1) +evo = CutFEMEvolve(V_φ,Ω,dΩ_act,hₕ;max_steps,γg=0.01) reinit1 = StabilisedReinit(V_φ,Ω,dΩ_act,hₕ;stabilisation_method=ArtificialViscosity(2.0)) reinit2 = StabilisedReinit(V_φ,Ω,dΩ_act,hₕ;stabilisation_method=GridapTopOpt.InteriorPenalty(V_φ,γg=1.0)) reinit = GridapTopOpt.MultiStageStabilisedReinit([reinit1,reinit2]) @@ -180,7 +178,7 @@ function has_oscillations(m,os_it) all(@.(abs(history[:C,it]) < 0.05vf)) && GridapTopOpt.default_has_oscillations(m,os_it) end optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh; - γ=γ_evo,verbose=true,constraint_names=[:Vol],converged,has_oscillations,Λ_update_tol=0) + γ=γ_evo,verbose=true,constraint_names=[:Vol],converged,has_oscillations) for (it,(uh,ph,dh),φh) in optimiser GC.gc() if iszero(it % iter_mod) diff --git a/scripts/Embedded/MWEs/GMSH_STAGGERED_TESTING-6-Brinkmann_stokes_P2-P1_Ersatz_elast_fsi.jl b/scripts/Embedded/MWEs/GMSH_STAGGERED_TESTING-6-Brinkmann_stokes_P2-P1_Ersatz_elast_fsi.jl new file mode 100644 index 00000000..7034fe5e --- /dev/null +++ b/scripts/Embedded/MWEs/GMSH_STAGGERED_TESTING-6-Brinkmann_stokes_P2-P1_Ersatz_elast_fsi.jl @@ -0,0 +1,221 @@ +using Gridap, Gridap.Geometry, Gridap.Adaptivity, Gridap.MultiField +using GridapEmbedded, GridapEmbedded.LevelSetCutters +using GridapGmsh +using GridapTopOpt + +using GridapSolvers + +############# + +path = "./results/GMSH_STAGGERED_TESTING-6-Brinkmann_stokes_P2-P1_Ersatz_elast_fsi/results/" +mkpath(path) + +γ_evo = 0.1 +max_steps = 20 +vf = 0.025 +α_coeff = 1 +iter_mod = 1 +D = 2 + +# Load gmsh mesh (Currently need to update mesh.geo and these values concurrently) +L = 2.0; +H = 0.5; +x0 = 0.5; +l = 0.4; +w = 0.025; +a = 0.3; +b = 0.01; + +model = GmshDiscreteModel((@__DIR__)*"/mesh.msh") +writevtk(model,path*"model") + +Ω_act = Triangulation(model) +hₕ = CellField(get_element_diameters(model),Ω_act) +hmin = minimum(get_element_diameters(model)) + +# Cut the background model +reffe_scalar = ReferenceFE(lagrangian,Float64,1) +V_φ = TestFESpace(model,reffe_scalar) +V_reg = TestFESpace(model,reffe_scalar;dirichlet_tags=["Omega_NonDesign","Gamma_s_D"]) +U_reg = TrialFESpace(V_reg) + +_e = 1e-3 +f0((x,y),W,H) = max(2/W*abs(x-x0),1/(H/2+1)*abs(y-H/2+1))-1 +f1((x,y),q,r) = - cos(q*π*x)*cos(q*π*y)/q - r/q +fin(x) = f0(x,l*(1+5_e),a*(1+5_e)) +fsolid(x) = min(f0(x,l*(1+_e),b*(1+_e)),f0(x,w*(1+_e),a*(1+_e))) +fholes((x,y),q,r) = max(f1((x,y),q,r),f1((x-1/q,y),q,r)) +φf(x) = min(max(fin(x),fholes(x,22,0.6)),fsolid(x)) +φh = interpolate(φf,V_φ) +writevtk(get_triangulation(φh),path*"initial_lsf",cellfields=["φ"=>φh,"h"=>hₕ]) + +# Setup integration meshes and measures +order = 2 +degree = 2*order + +# Ω_act = Triangulation(model) +dΩ_act = Measure(Ω_act,degree) +vol_D = L*H + +Γf_D = BoundaryTriangulation(model,tags="Gamma_f_D") +Γf_N = BoundaryTriangulation(model,tags="Gamma_f_N") +dΓf_D = Measure(Γf_D,degree) +dΓf_N = Measure(Γf_N,degree) +Ω = EmbeddedCollection(model,φh) do cutgeo,_ + Ωs = DifferentiableTriangulation(Triangulation(cutgeo,PHYSICAL),V_φ) + Ωf = DifferentiableTriangulation(Triangulation(cutgeo,PHYSICAL_OUT),V_φ) + Γ = DifferentiableTriangulation(EmbeddedBoundary(cutgeo),V_φ) + Γg = GhostSkeleton(cutgeo) + Ωact = Triangulation(cutgeo,ACTIVE) + (; + :Ωs => Ωs, + :dΩs => Measure(Ωs,degree), + :Ωf => Ωf, + :dΩf => Measure(Ωf,degree), + :Γg => Γg, + :dΓg => Measure(Γg,degree), + :n_Γg => get_normal_vector(Γg), + :Γ => Γ, + :dΓ => Measure(Γ,degree), + :n_Γ => get_normal_vector(Γ.trian), + :Ωact => Ωact + ) +end + +# Setup FESpace +uin(x) = VectorValue(16x[2]*(H-x[2]),0.0) + +reffe_u = ReferenceFE(lagrangian,VectorValue{D,Float64},order,space=:P) +reffe_p = ReferenceFE(lagrangian,Float64,order-1,space=:P) +reffe_d = ReferenceFE(lagrangian,VectorValue{D,Float64},order-1) + +V = TestFESpace(Ω_act,reffe_u,conformity=:H1, + dirichlet_tags=["Gamma_f_D","Gamma_NoSlipTop","Gamma_NoSlipBottom","Gamma_s_D"]) +Q = TestFESpace(Ω_act,reffe_p,conformity=:C0) +T = TestFESpace(Ω_act ,reffe_d,conformity=:H1,dirichlet_tags=["Gamma_s_D"]) + +U = TrialFESpace(V,[uin,VectorValue(0.0,0.0),VectorValue(0.0,0.0),VectorValue(0.0,0.0)]) +P = TrialFESpace(Q) +R = TrialFESpace(T) + +# mfs_VQ = BlockMultiFieldStyle(2,(1,1)) #<- allows to specify block structure for first problem +VQ = MultiFieldFESpace([V,Q])#;style=mfs_VQ) +UP = MultiFieldFESpace([U,P])#;style=mfs_VQ) + +mfs = BlockMultiFieldStyle(2,(2,1)) +X = MultiFieldFESpace([U,P,R];style=mfs) + +# Weak form +## Fluid +# Properties +Re = 60 # Reynolds number +ρ = 1.0 # Density +cl = H # Characteristic length +u0_max = maximum(abs,get_dirichlet_dof_values(U)) +μ = ρ*cl*u0_max/Re # Viscosity +ν = μ/ρ # Kinematic viscosity +# Stabilization parameters +γ(h) = 1000/h #1e5*μ/h + +# Terms +σf_n(u,p,n) = μ*∇(u) ⋅ n - p*n +a_Ω(u,v) = μ*(∇(u) ⊙ ∇(v)) +b_Ω(v,p) = - (∇⋅v)*p + +a_fluid((),(u,p),(v,q)) = + ∫( a_Ω(u,v)+b_Ω(u,q)+b_Ω(v,p)) * Ω.dΩf + + ∫( a_Ω(u,v)+b_Ω(u,q)+b_Ω(v,p) + (γ ∘ hₕ)*u⋅v) * Ω.dΩs +l_fluid((),(v,q)) = 0.0 + +## Structure +# Stabilization and material parameters +function lame_parameters(E,ν) + λ = (E*ν)/((1+ν)*(1-2*ν)) + μ = E/(2*(1+ν)) + (λ, μ) +end +λs, μs = lame_parameters(0.1,0.05) +ϵ = (λs + 2μs)*1e-3 +# Terms +σ(ε) = λs*tr(ε)*one(ε) + 2*μs*ε +a_solid(((u,p),),d,s) = ∫(ε(s) ⊙ (σ ∘ ε(d)))Ω.dΩs + + ∫(ϵ*(ε(s) ⊙ (σ ∘ ε(d))))Ω.dΩf # Ersatz + +function l_solid(((u,p),),s) + n_AD = get_normal_vector(Ω.Γ) + return ∫(σf_n(u,p,n_AD) ⋅ s)Ω.dΓ +end + +## FEOperator way +fluid_op = AffineFEOperator((x,y)->a_fluid((),x,y),y->l_fluid((),y),UP,VQ) +uh_1,ph_1 = solve(fluid_op) +solid_op = AffineFEOperator((x,y)->a_solid(((uh_1,ph_1),),x,y),y->l_solid(((uh_1,ph_1),),y),R,T) +dh_1 = solve(solid_op) + +## Staggered +op = GridapSolvers.StaggeredAffineFEOperator([a_fluid,a_solid],[l_fluid,l_solid],[UP,R],[VQ,T]) + + +xh = zero(X) +# lsolver = CGSolver(JacobiLinearSolver();rtol=1.e-12,verbose=true) +solver = GridapSolvers.StaggeredFESolver(fill(LUSolver(),2)) +xh,cache = solve!(xh,solver,op) +uh_2,ph_2,dh_2 = xh + +writevtk(Ω.Ωf,path*"Omega_f", + cellfields=["uh"=>uh,"ph"=>ph,"dh"=>dh]) +writevtk(Ω.Ωs,path*"Omega_s", + cellfields=["uh"=>uh,"ph"=>ph,"dh"=>dh]) + +norm(get_free_dof_values(uh_1)-get_free_dof_values(uh_2),Inf) +norm(get_free_dof_values(ph_1)-get_free_dof_values(ph_2),Inf) +norm(get_free_dof_values(dh_1)-get_free_dof_values(dh_2),Inf) + +## Old way +function old_main() + _X = MultiFieldFESpace([U,P,R]) + _Y = MultiFieldFESpace([V,Q,T]) + + # Terms + _a_fluid((u,p),(v,q)) = + ∫( a_Ω(u,v)+b_Ω(u,q)+b_Ω(v,p)) * Ω.dΩf + + ∫( a_Ω(u,v)+b_Ω(u,q)+b_Ω(v,p) + (γ ∘ hₕ)*u⋅v ) * Ω.dΩs + + ## Structure + _a_solid(d,s) = ∫(ε(s) ⊙ (σ ∘ ε(d)))Ω.dΩs + + ∫(ϵ*(ε(s) ⊙ (σ ∘ ε(d))))Ω.dΩf # Ersatz + + ## Full problem + function a_coupled((u,p,d),(v,q,s)) + n_AD = get_normal_vector(Ω.Γ) + return _a_fluid((u,p),(v,q)) + _a_solid(d,s) + + ∫(-σf_n(u,p,n_AD) ⋅ s)Ω.dΓ + end + l_coupled((v,q,s)) = 0.0 + + return AffineFEOperator(a_coupled,l_coupled,_X,_Y) +end + +op_old = old_main() + +_xh = solve(op_old) +_uh,_ph,_dh = _xh + +norm(get_free_dof_values(uh)-get_free_dof_values(_uh),Inf) +norm(get_free_dof_values(ph)-get_free_dof_values(_ph),Inf) +norm(get_free_dof_values(dh)-get_free_dof_values(_dh),Inf) + + + +######## +_f1((),x,ϕ) = ϕ +_f2((x,),y,ϕ) = ϕ +_f3((x,y,),z,ϕ) = ϕ + +_fvec0 = [_f1,_f2,_f3] + +_fvec = map(x->((xs,u)->x(xs,u,0)),_fvec0) + +_fvec[1]((),1) +_fvec[2]((1,),2) +_fvec[3]((1,1),2) \ No newline at end of file diff --git a/src/LevelSetEvolution/UnfittedEvolution/CutFEMEvolve.jl b/src/LevelSetEvolution/UnfittedEvolution/CutFEMEvolve.jl index 5de013a6..ed083501 100644 --- a/src/LevelSetEvolution/UnfittedEvolution/CutFEMEvolve.jl +++ b/src/LevelSetEvolution/UnfittedEvolution/CutFEMEvolve.jl @@ -106,8 +106,12 @@ function get_transient_operator(φh,velh,s::CutFEMEvolve) # In particular, this happens when we have jumps, where some contributions come from two # cells away. Boundary cells then get contributions from cells which are not in the local domain. Ut_φ = TransientTrialFESpace(V_φ) + + # TODO: This has been disabled due to bug. See below discussion. + # ode_op = TransientLinearFEOperator((stiffness,mass),forcing,Ut_φ,V_φ; + # constant_forms=(false,true),assembler) ode_op = TransientLinearFEOperator((stiffness,mass),forcing,Ut_φ,V_φ; - constant_forms=(false,true),assembler) + constant_forms=(true,true),assembler) return ode_op end @@ -135,7 +139,9 @@ function solve!(s::CutFEMEvolve,φh,velh,γ,cache::Nothing) # March march = Base.iterate(ode_sol) data, state = march - state_new = update_reuse!(state,true) + # state_new = update_reuse!(state,true) # TODO: This has been disabled due to bug. See below discussion. + state_new = state + march_new = data, state_new while march_new !== nothing data, state_new = march_new @@ -145,8 +151,12 @@ function solve!(s::CutFEMEvolve,φh,velh,γ,cache::Nothing) # Update φh and cache _, φhF = data copy!(get_free_dof_values(φh),get_free_dof_values(φhF)) - s.cache = state_new - update_collection!(s.Ωs,φh) # TODO: remove? + # TODO: This has been disabled for the time being. Originally when this code + # was written, we expected that changing reuse to false and iterating once + # would update the stiffness matrix. However, this does not appear to be the case. + # See `scripts/Embedded/Bugs/odes_reassemble_stiffness.jl` for a minimal example. + # s.cache = state_new + update_collection!(s.Ωs,φh) return φh end diff --git a/src/LevelSetEvolution/UnfittedEvolution/StabilisedReinit.jl b/src/LevelSetEvolution/UnfittedEvolution/StabilisedReinit.jl index 4e0c5867..6f4d62d8 100644 --- a/src/LevelSetEvolution/UnfittedEvolution/StabilisedReinit.jl +++ b/src/LevelSetEvolution/UnfittedEvolution/StabilisedReinit.jl @@ -163,8 +163,10 @@ struct MultiStageStabilisedReinit <: Reinitialiser stages::Vector{StabilisedReinit} function MultiStageStabilisedReinit(stages::Vector{<:StabilisedReinit}) h = get_element_diameters(first(stages)) + V_φ = get_space(first(stages)) for stage in stages @check h === get_element_diameters(stage) + @check V_φ === get_space(stage) end new(stages) end @@ -172,6 +174,7 @@ end get_cache(s::MultiStageStabilisedReinit) = get_cache.(s.stages) get_element_diameters(s::MultiStageStabilisedReinit) = get_element_diameters(first(s.stages)) +get_space(s::MultiStageStabilisedReinit) = get_space(first(s.stages)) function solve!(s::MultiStageStabilisedReinit,φh,caches) for (stage,cache) in zip(s.stages,caches)