diff --git a/compile/warmup.jl b/compile/warmup.jl index c79ecbd1..c2b6e5e3 100644 --- a/compile/warmup.jl +++ b/compile/warmup.jl @@ -70,7 +70,7 @@ function main(mesh_partition,distribute) dJ = (q,u,φ,dΩ,dΓ_N) -> ∫((ξ-D*∇(u)⋅∇(u))*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ; ## Finite difference solver and level set function - stencil = AdvectionStencil(FirstOrderStencil(2,Float64),model,V_φ,Δ./order,max_steps,tol) + ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,Δ./order,max_steps,tol) reinit!(stencil,φ,γ_reinit) ## Setup solver and FE operators @@ -84,7 +84,7 @@ function main(mesh_partition,distribute) ## Optimiser _conv_cond = t->LevelSetTopOpt.conv_cond(t;coef=1/50); - optimiser = AugmentedLagrangian(φ,pcfs,stencil,vel_ext,interp,el_size,γ,γ_reinit,conv_criterion=_conv_cond); + optimiser = AugmentedLagrangian(φ,pcfs,ls_evo,vel_ext,interp,el_size,γ,γ_reinit,conv_criterion=_conv_cond); for history in optimiser it,Ji,_,_ = last(history) print_history(it,["J"=>Ji];ranks=ranks) diff --git a/docs/make.jl b/docs/make.jl index 67bba876..6edb2324 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -30,7 +30,7 @@ makedocs( "Reference" => [ "reference/optimisers.md", "reference/chainrules.md", - "reference/advection.md", + "reference/levelsetevolution.md", "reference/velext.md", "reference/io.md", "reference/utilities.md", diff --git a/docs/src/reference/advection.md b/docs/src/reference/advection.md deleted file mode 100644 index 46166d58..00000000 --- a/docs/src/reference/advection.md +++ /dev/null @@ -1,25 +0,0 @@ -# Advection - -## `AdvectionStencil` -```@docs -LevelSetTopOpt.AdvectionStencil -LevelSetTopOpt.AdvectionStencil(stencil::LevelSetTopOpt.Stencil,model,space,tol=1.e-3,max_steps=100,max_steps_reinit=2000) -LevelSetTopOpt.advect! -LevelSetTopOpt.reinit! -``` - -## Stencils - -```@docs -LevelSetTopOpt.FirstOrderStencil -``` - -## Custom `Stencil` - -```@docs -LevelSetTopOpt.Stencil -LevelSetTopOpt.advect!(::LevelSetTopOpt.Stencil,φ,vel,Δt,Δx,isperiodic,caches) -LevelSetTopOpt.reinit!(::LevelSetTopOpt.Stencil,φ_new,φ,vel,Δt,Δx,isperiodic,caches) -LevelSetTopOpt.compute_Δt(::LevelSetTopOpt.Stencil,φ,vel) -LevelSetTopOpt.allocate_caches(::LevelSetTopOpt.Stencil,φ,vel) -``` \ No newline at end of file diff --git a/docs/src/reference/levelsetevolution.md b/docs/src/reference/levelsetevolution.md new file mode 100644 index 00000000..3f41a926 --- /dev/null +++ b/docs/src/reference/levelsetevolution.md @@ -0,0 +1,39 @@ +# LevelSetEvolution +In LevelSetTopOpt, the level set is evolved and reinitialised using a `LevelSetEvolution` method. The most standard of these is the Hamilton-Jacobi evolution equation solved using a first order upwind finite difference scheme. A forward Euler in time method is provided below via `HamiltonJacobiEvolution <: LevelSetEvolution` along with an upwind finite difference stencil for the spatial discretisation via `FirstOrderStencil`. + +This can be extended in several ways. For example, higher order spatial stencils can be implemented by extending the `Stencil` interface below. In addition, more advanced ODE solvers could be implemented (e.g., Runge–Kutta methods) or entirely different level set evolution methods by extending the `LevelSetEvolution` interface below. + +## `HamiltonJacobiEvolution` +```@docs +LevelSetTopOpt.HamiltonJacobiEvolution +LevelSetTopOpt.HamiltonJacobiEvolution(stencil::LevelSetTopOpt.Stencil,model,space,tol=1.e-3,max_steps=100,max_steps_reinit=2000) +LevelSetTopOpt.evolve! +LevelSetTopOpt.reinit! +LevelSetTopOpt.get_dof_Δ(m::HamiltonJacobiEvolution) +``` + +## Spatial stencils for `HamiltonJacobiEvolution` + +```@docs +LevelSetTopOpt.FirstOrderStencil +``` + +## Custom `Stencil` + +```@docs +LevelSetTopOpt.Stencil +LevelSetTopOpt.evolve!(::LevelSetTopOpt.Stencil,φ,vel,Δt,Δx,isperiodic,caches) +LevelSetTopOpt.reinit!(::LevelSetTopOpt.Stencil,φ_new,φ,vel,Δt,Δx,isperiodic,caches) +LevelSetTopOpt.allocate_caches(::LevelSetTopOpt.Stencil,φ,vel) +LevelSetTopOpt.check_order +``` + +## Custom `LevelSetEvolution` +To implement a custom level set evolution method, we can extend the methods below. For example, one could consider Reaction-Diffusion-based evolution of the level set function. This can be solved with a finite element method and so we can implement a new type that inherits from `LevelSetEvolution` independently of the `Stencil` types. + +```@docs +LevelSetTopOpt.LevelSetEvolution +LevelSetTopOpt.evolve!(::LevelSetTopOpt.LevelSetEvolution,φ,args...) +LevelSetTopOpt.reinit!(::LevelSetTopOpt.LevelSetEvolution,φ,args...) +LevelSetTopOpt.get_dof_Δ(::LevelSetTopOpt.LevelSetEvolution) +``` \ No newline at end of file diff --git a/docs/src/reference/optimisers.md b/docs/src/reference/optimisers.md index f886ccf2..502759d3 100644 --- a/docs/src/reference/optimisers.md +++ b/docs/src/reference/optimisers.md @@ -23,7 +23,7 @@ LevelSetTopOpt.OptimiserHistory LevelSetTopOpt.OptimiserHistorySlice ``` -## Custom optimiser +## Custom `Optimiser` ```@docs LevelSetTopOpt.Optimiser LevelSetTopOpt.iterate(::LevelSetTopOpt.Optimiser) diff --git a/docs/src/tutorials/minimum_thermal_compliance.md b/docs/src/tutorials/minimum_thermal_compliance.md index ddbb4337..3d52db3c 100644 --- a/docs/src/tutorials/minimum_thermal_compliance.md +++ b/docs/src/tutorials/minimum_thermal_compliance.md @@ -346,10 +346,10 @@ Both of these equations can be solved numerically on a Cartesian mesh using a fi ```julia # Finite difference scheme scheme = FirstOrderStencil(length(el_size),Float64) -stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps) +ls_evo = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps) ``` -In the above we first build an object [`FirstOrderStencil`](@ref) that represents a finite difference stencil for a single step of the Hamilton-Jacobi evolution equation and reinitialisation equation. We use `length(el_size)` to indicate the dimension of the problem. We then create an [`AdvectionStencil`](@ref) which enables finite differencing on order `O` finite elements in serial or parallel. The [`AdvectionStencil`](@ref) object provides two important methods [`advect!`](@ref) and [`reinit!`](@ref) that correspond to solving the Hamilton-Jacobi evolution equation and reinitialisation equation, respectively. +In the above we first build an object [`FirstOrderStencil`](@ref) that represents a finite difference stencil for a single step of the Hamilton-Jacobi evolution equation and reinitialisation equation. We use `length(el_size)` to indicate the dimension of the problem. We then create an [`HamiltonJacobiEvolution`](@ref) which enables finite differencing on order `O` finite elements in serial or parallel. The [`HamiltonJacobiEvolution`](@ref) object provides two important methods [`evolve!`](@ref) and [`reinit!`](@ref) that correspond to solving the Hamilton-Jacobi evolution equation and reinitialisation equation, respectively. ### Optimiser, visualisation and IO @@ -357,7 +357,7 @@ We may now create the optimiser object. This structure holds all information reg ```julia # Optimiser -optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit,verbose=true,constraint_names=[:Vol]) +optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;γ,γ_reinit,verbose=true,constraint_names=[:Vol]) ``` As optimisers inheriting from [`LevelSetTopOpt.Optimiser`](@ref) implement Julia's iterator functionality, we can solve the optimisation problem to convergence by iterating over the optimiser: @@ -461,9 +461,9 @@ a_hilb(p,q) =∫(α^2*∇(p)⋅∇(q) + p*q)dΩ vel_ext = VelocityExtension(a_hilb,U_reg,V_reg) # Finite difference scheme scheme = FirstOrderStencil(length(el_size),Float64) -stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps) +ls_evo = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps) # Optimiser -optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit,verbose=true,constraint_names=[:Vol]) +optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;γ,γ_reinit,verbose=true,constraint_names=[:Vol]) # Solve for (it,uh,φh) in optimiser data = ["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh] @@ -604,9 +604,9 @@ a_hilb(p,q) =∫(α^2*∇(p)⋅∇(q) + p*q)dΩ vel_ext = VelocityExtension(a_hilb,U_reg,V_reg) # Finite difference scheme scheme = FirstOrderStencil(length(el_size),Float64) -stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps) +ls_evo = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps) # Optimiser -optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit,verbose=true,constraint_names=[:Vol]) +optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;γ,γ_reinit,verbose=true,constraint_names=[:Vol]) # Solve for (it,uh,φh) in optimiser data = ["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh] @@ -754,9 +754,9 @@ function main() ) # Finite difference scheme scheme = FirstOrderStencil(length(el_size),Float64) - stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps) + ls_evo = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps) # Optimiser - optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit,verbose=true,constraint_names=[:Vol]) + optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;γ,γ_reinit,verbose=true,constraint_names=[:Vol]) # Solve for (it,uh,φh) in optimiser data = ["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh] @@ -827,7 +827,7 @@ We then adjust lines 28-31 as follows: The function `i_am_main` returns true only on the first processor. This function is useful for ensuring certain operations only happen once instead of several times across each executable. In addition, we now create a partitioned Cartesian model using `CartesianDiscreteModel(ranks,mesh_partition,dom,el_size)`. Finally, we adjust line 82 to ensure that verbosity only happens on the first processors: ```julia -optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh; +optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh; γ,γ_reinit,verbose=i_am_main(ranks),constraint_names=[:Vol]) ``` @@ -918,9 +918,9 @@ function main(mesh_partition,distribute) ) # Finite difference scheme scheme = FirstOrderStencil(length(el_size),Float64) - stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps) + ls_evo = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps) # Optimiser - optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh; + optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh; γ,γ_reinit,verbose=i_am_main(ranks),constraint_names=[:Vol]) # Solve for (it,uh,φh) in optimiser @@ -1070,9 +1070,9 @@ a_hilb(p,q) =∫(α^2*∇(p)⋅∇(q) + p*q)dΩ vel_ext = VelocityExtension(a_hilb,U_reg,V_reg) # Finite difference scheme scheme = FirstOrderStencil(length(el_size),Float64) -stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps) +ls_evo = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps) # Optimiser -optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit,verbose=true,constraint_names=[:Vol]) +optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;γ,γ_reinit,verbose=true,constraint_names=[:Vol]) # Solve for (it,uh,φh) in optimiser data = ["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh] diff --git a/scripts/Benchmarks/benchmark.jl b/scripts/Benchmarks/benchmark.jl index 3c0fe455..862262d7 100644 --- a/scripts/Benchmarks/benchmark.jl +++ b/scripts/Benchmarks/benchmark.jl @@ -102,9 +102,9 @@ function nl_elast(mesh_partition,ranks,el_size,order,verbose) ) # Finite difference scheme scheme = FirstOrderStencil(length(el_size),Float64) - stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps) + ls_evo = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps) # Optimiser - return AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit,verbose) + return AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;γ,γ_reinit,verbose) end function therm(mesh_partition,ranks,el_size,order,verbose) @@ -180,9 +180,9 @@ function therm(mesh_partition,ranks,el_size,order,verbose) ) # Finite difference scheme scheme = FirstOrderStencil(length(el_size),Float64) - stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps) + ls_evo = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps) # Optimiser - return AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit,verbose) + return AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;γ,γ_reinit,verbose) end function elast(mesh_partition,ranks,el_size,order,verbose) @@ -256,9 +256,9 @@ function elast(mesh_partition,ranks,el_size,order,verbose) ) # Finite difference scheme scheme = FirstOrderStencil(length(el_size),Float64) - stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps) + ls_evo = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps) # Optimiser - return AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit,verbose) + return AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;γ,γ_reinit,verbose) end function inverter_HPM(mesh_partition,ranks,el_size,order,verbose) @@ -329,7 +329,7 @@ function inverter_HPM(mesh_partition,ranks,el_size,order,verbose) UΓ_out(u,φ,dΩ,dΓ_in,dΓ_out) = ∫((u⋅-e₁-δₓ)/vol_Γ_out)dΓ_out ## Finite difference solver - stencil = AdvectionStencil(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps) + ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps) ## Setup solver and FE operators Tm = SparseMatrixCSR{0,PetscScalar,PetscInt} @@ -355,7 +355,7 @@ function inverter_HPM(mesh_partition,ranks,el_size,order,verbose) ) ## Optimiser - return HilbertianProjection(pcfs,stencil,vel_ext,φh;γ,γ_reinit,verbose=verbose) + return HilbertianProjection(pcfs,ls_evo,vel_ext,φh;γ,γ_reinit,verbose=verbose) end with_mpi() do distribute diff --git a/scripts/MPI/3d_elastic_compliance_ALM.jl b/scripts/MPI/3d_elastic_compliance_ALM.jl index 58836523..558d05d8 100644 --- a/scripts/MPI/3d_elastic_compliance_ALM.jl +++ b/scripts/MPI/3d_elastic_compliance_ALM.jl @@ -81,7 +81,7 @@ function main(mesh_partition,distribute,el_size) dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ ## Finite difference solver and level set function - stencil = AdvectionStencil(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps) + ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps) ## Setup solver and FE operators Tm = SparseMatrixCSR{0,PetscScalar,PetscInt} @@ -107,7 +107,7 @@ function main(mesh_partition,distribute,el_size) ) ## Optimiser - optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh; + optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh; γ,γ_reinit,verbose=i_am_main(ranks),constraint_names=[:Vol]) for (it, uh, φh) in optimiser data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh] diff --git a/scripts/MPI/3d_hyperelastic_compliance_ALM.jl b/scripts/MPI/3d_hyperelastic_compliance_ALM.jl index 6064ae61..6cd67d04 100644 --- a/scripts/MPI/3d_hyperelastic_compliance_ALM.jl +++ b/scripts/MPI/3d_hyperelastic_compliance_ALM.jl @@ -83,7 +83,7 @@ function main(mesh_partition,distribute,el_size) dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ ## Finite difference solver and level set function - stencil = AdvectionStencil(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps) + ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps) ## Setup solver and FE operators Tm = SparseMatrixCSR{0,PetscScalar,PetscInt} @@ -110,7 +110,7 @@ function main(mesh_partition,distribute,el_size) ) ## Optimiser - optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh; + optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh; γ,γ_reinit,verbose=i_am_main(ranks),constraint_names=[:Vol]) for (it, uh, φh) in optimiser data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh] diff --git a/scripts/MPI/3d_hyperelastic_compliance_neohook_ALM.jl b/scripts/MPI/3d_hyperelastic_compliance_neohook_ALM.jl index 7b02dd0c..8ad44dc4 100644 --- a/scripts/MPI/3d_hyperelastic_compliance_neohook_ALM.jl +++ b/scripts/MPI/3d_hyperelastic_compliance_neohook_ALM.jl @@ -97,7 +97,7 @@ function main(mesh_partition,distribute,el_size) dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ ## Finite difference solver and level set function - stencil = AdvectionStencil(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps) + ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps) ## Setup solver and FE operators Tm = SparseMatrixCSR{0,PetscScalar,PetscInt} @@ -124,7 +124,7 @@ function main(mesh_partition,distribute,el_size) ) ## Optimiser - optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh; + optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh; γ,γ_reinit,verbose=i_am_main(ranks),constraint_names=[:Vol]) for (it, uh, φh) in optimiser data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh] diff --git a/scripts/MPI/3d_inverse_homenisation_ALM.jl b/scripts/MPI/3d_inverse_homenisation_ALM.jl index ec6f7872..1746a942 100644 --- a/scripts/MPI/3d_inverse_homenisation_ALM.jl +++ b/scripts/MPI/3d_inverse_homenisation_ALM.jl @@ -85,7 +85,7 @@ function main(mesh_partition,distribute,el_size) dVol(q,u,φ,dΩ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ ## Finite difference solver and level set function - stencil = AdvectionStencil(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps) + ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps) ## Setup solver and FE operators Tm = SparseMatrixCSR{0,PetscScalar,PetscInt} @@ -111,7 +111,7 @@ function main(mesh_partition,distribute,el_size) ) ## Optimiser - optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh; + optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh; γ,γ_reinit,verbose=i_am_main(ranks),constraint_names=[:Vol]) for (it, uh, φh) in optimiser data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh))] diff --git a/scripts/MPI/3d_inverter_ALM.jl b/scripts/MPI/3d_inverter_ALM.jl index 0b600003..4fa6ee26 100644 --- a/scripts/MPI/3d_inverter_ALM.jl +++ b/scripts/MPI/3d_inverter_ALM.jl @@ -99,7 +99,7 @@ function main(mesh_partition,distribute,el_size) UΓ_out(u,φ,dΩ,dΓ_in,dΓ_out) = ∫((u⋅-e₁-δₓ)/vol_Γ_out)dΓ_out ## Finite difference solver and level set function - stencil = AdvectionStencil(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps) + ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps) ## Setup solver and FE operators Tm = SparseMatrixCSR{0,PetscScalar,PetscInt} @@ -125,7 +125,7 @@ function main(mesh_partition,distribute,el_size) ) ## Optimiser - optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh; + optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh; γ,γ_reinit,verbose=i_am_main(ranks),constraint_names=[:Vol,:UΓ_out]) for (it, uh, φh) in optimiser data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh] diff --git a/scripts/MPI/3d_inverter_HPM.jl b/scripts/MPI/3d_inverter_HPM.jl index bd6d16e9..3aacbfe1 100644 --- a/scripts/MPI/3d_inverter_HPM.jl +++ b/scripts/MPI/3d_inverter_HPM.jl @@ -99,7 +99,7 @@ function main(mesh_partition,distribute,el_size) UΓ_out(u,φ,dΩ,dΓ_in,dΓ_out) = ∫((u⋅-e₁-δₓ)/vol_Γ_out)dΓ_out ## Finite difference solver and level set function - stencil = AdvectionStencil(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps) + ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps) ## Setup solver and FE operators Tm = SparseMatrixCSR{0,PetscScalar,PetscInt} @@ -126,7 +126,7 @@ function main(mesh_partition,distribute,el_size) ## Optimiser ls_enabled=false # Setting to true will use a line search instead of oscillation detection - optimiser = HilbertianProjection(pcfs,stencil,vel_ext,φh;γ,γ_reinit,ls_enabled,α_min=0.7, + optimiser = HilbertianProjection(pcfs,ls_evo,vel_ext,φh;γ,γ_reinit,ls_enabled,α_min=0.7, ls_γ_max=0.05,verbose=i_am_main(ranks),constraint_names=[:Vol,:UΓ_out]) for (it, uh, φh) in optimiser data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh] diff --git a/scripts/MPI/3d_nonlinear_thermal_compliance_ALM.jl b/scripts/MPI/3d_nonlinear_thermal_compliance_ALM.jl index f4f110d6..e3bbfde2 100644 --- a/scripts/MPI/3d_nonlinear_thermal_compliance_ALM.jl +++ b/scripts/MPI/3d_nonlinear_thermal_compliance_ALM.jl @@ -85,7 +85,7 @@ function main(mesh_partition,distribute,el_size) dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ ## Finite difference solver and level set function - stencil = AdvectionStencil(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps) + ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps) ## Setup solver and FE operators Tm = SparseMatrixCSR{0,PetscScalar,PetscInt} @@ -113,7 +113,7 @@ function main(mesh_partition,distribute,el_size) ) ## Optimiser - optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh; + optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh; γ,γ_reinit,verbose=i_am_main(ranks),constraint_names=[:Vol]) for (it, uh, φh) in optimiser data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh] diff --git a/scripts/MPI/3d_thermal_compliance_ALM.jl b/scripts/MPI/3d_thermal_compliance_ALM.jl index 0fa01deb..2c63d3e0 100644 --- a/scripts/MPI/3d_thermal_compliance_ALM.jl +++ b/scripts/MPI/3d_thermal_compliance_ALM.jl @@ -81,7 +81,7 @@ function main(mesh_partition,distribute,el_size) dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ ## Finite difference solver and level set function - stencil = AdvectionStencil(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps) + ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps) ## Setup solver and FE operators Tm = SparseMatrixCSR{0,PetscScalar,PetscInt} @@ -107,7 +107,7 @@ function main(mesh_partition,distribute,el_size) ) ## Optimiser - optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh; + optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh; γ,γ_reinit,verbose=i_am_main(ranks),constraint_names=[:Vol]) for (it, uh, φh) in optimiser data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh] diff --git a/scripts/MPI/elastic_compliance_ALM.jl b/scripts/MPI/elastic_compliance_ALM.jl index fcfb451f..76a48e2b 100644 --- a/scripts/MPI/elastic_compliance_ALM.jl +++ b/scripts/MPI/elastic_compliance_ALM.jl @@ -74,7 +74,7 @@ function main(mesh_partition,distribute,el_size) dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ ## Finite difference solver and level set function - stencil = AdvectionStencil(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) + ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) ## Setup solver and FE operators Tm = SparseMatrixCSR{0,PetscScalar,PetscInt} @@ -98,7 +98,7 @@ function main(mesh_partition,distribute,el_size) ls=PETScLinearSolver()) ## Optimiser - optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh; + optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh; γ,γ_reinit,verbose=i_am_main(ranks),constraint_names=[:Vol]) for (it, uh, φh) in optimiser data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh] diff --git a/scripts/MPI/inverse_homenisation_ALM.jl b/scripts/MPI/inverse_homenisation_ALM.jl index 7a068021..55fe3af0 100644 --- a/scripts/MPI/inverse_homenisation_ALM.jl +++ b/scripts/MPI/inverse_homenisation_ALM.jl @@ -73,7 +73,7 @@ function main(mesh_partition,distribute,el_size) dVol(q,u,φ,dΩ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ ## Finite difference solver and level set function - stencil = AdvectionStencil(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) + ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) ## Setup solver and FE operators Tm = SparseMatrixCSR{0,PetscScalar,PetscInt} @@ -99,7 +99,7 @@ function main(mesh_partition,distribute,el_size) ) ## Optimiser - optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh; + optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh; γ,γ_reinit,verbose=i_am_main(ranks),constraint_names=[:Vol]) for (it, uh, φh) in optimiser data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh))] diff --git a/scripts/MPI/thermal_compliance_ALM.jl b/scripts/MPI/thermal_compliance_ALM.jl index 8bb1e9aa..db881e3e 100644 --- a/scripts/MPI/thermal_compliance_ALM.jl +++ b/scripts/MPI/thermal_compliance_ALM.jl @@ -72,7 +72,7 @@ function main(mesh_partition,distribute) dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ ## Finite difference solver and level set function - stencil = AdvectionStencil(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) + ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) ## Setup solver and FE operators state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N) @@ -84,7 +84,7 @@ function main(mesh_partition,distribute) vel_ext = VelocityExtension(a_hilb,U_reg,V_reg) ## Optimiser - optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit, + optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;γ,γ_reinit, verbose=i_am_main(ranks),constraint_names=[:Vol]) for (it, uh, φh) in optimiser data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh] diff --git a/scripts/MPI/thermal_compliance_ALM_PETSc.jl b/scripts/MPI/thermal_compliance_ALM_PETSc.jl index 39bd0271..366f13ef 100644 --- a/scripts/MPI/thermal_compliance_ALM_PETSc.jl +++ b/scripts/MPI/thermal_compliance_ALM_PETSc.jl @@ -72,7 +72,7 @@ function main(mesh_partition,distribute,el_size) dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ ## Finite difference solver and level set function - stencil = AdvectionStencil(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) + ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) ## Setup solver and FE operators Tm = SparseMatrixCSR{0,PetscScalar,PetscInt} @@ -98,7 +98,7 @@ function main(mesh_partition,distribute,el_size) ) ## Optimiser - optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit, + optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;γ,γ_reinit, verbose=i_am_main(ranks),constraint_names=[:Vol]) for (it, uh, φh) in optimiser data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh] diff --git a/scripts/Serial/elastic_compliance_ALM.jl b/scripts/Serial/elastic_compliance_ALM.jl index c313ef68..6b8e4d17 100644 --- a/scripts/Serial/elastic_compliance_ALM.jl +++ b/scripts/Serial/elastic_compliance_ALM.jl @@ -71,7 +71,7 @@ function main() dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ ## Finite difference solver and level set function - stencil = AdvectionStencil(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) + ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) ## Setup solver and FE operators state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N) @@ -83,7 +83,7 @@ function main() vel_ext = VelocityExtension(a_hilb,U_reg,V_reg) ## Optimiser - optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh; + optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh; γ,γ_reinit,verbose=true,constraint_names=[:Vol]) for (it,uh,φh) in optimiser data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh] diff --git a/scripts/Serial/elastic_compliance_HPM.jl b/scripts/Serial/elastic_compliance_HPM.jl index e7fe84c8..d13be2c1 100644 --- a/scripts/Serial/elastic_compliance_HPM.jl +++ b/scripts/Serial/elastic_compliance_HPM.jl @@ -72,7 +72,7 @@ function main() dVol = (q,u,φ,dΩ,dΓ_N) -> ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ ## Finite difference solver and level set function - stencil = AdvectionStencil(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) + ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) ## Setup solver and FE operators state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N) @@ -84,7 +84,7 @@ function main() vel_ext = VelocityExtension(a_hilb,U_reg,V_reg) ## Optimiser - optimiser = HilbertianProjection(pcfs,stencil,vel_ext,φh; + optimiser = HilbertianProjection(pcfs,ls_evo,vel_ext,φh; γ,γ_reinit,verbose=true,constraint_names=[:Vol]) for (it,uh,φh) in optimiser data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh] diff --git a/scripts/Serial/elastic_compliance_LM.jl b/scripts/Serial/elastic_compliance_LM.jl index 004e2119..3accdfac 100644 --- a/scripts/Serial/elastic_compliance_LM.jl +++ b/scripts/Serial/elastic_compliance_LM.jl @@ -68,7 +68,7 @@ function main() dJ = (q,u,φ,dΩ,dΓ_N) -> ∫((- ξ + C ⊙ ε(u) ⊙ ε(u))*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ; ## Finite difference solver and level set function - stencil = AdvectionStencil(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) + ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) ## Setup solver and FE operators state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N) @@ -80,7 +80,7 @@ function main() vel_ext = VelocityExtension(a_hilb,U_reg,V_reg) ## Optimiser - optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit,verbose=true) + optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;γ,γ_reinit,verbose=true) for (it, uh, φh) in optimiser data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh] iszero(it % iter_mod) && writevtk(Ω,path*"out$it",cellfields=data) diff --git a/scripts/Serial/hyperelastic_compliance_ALM.jl b/scripts/Serial/hyperelastic_compliance_ALM.jl index fe45ee8e..63b9e1bf 100644 --- a/scripts/Serial/hyperelastic_compliance_ALM.jl +++ b/scripts/Serial/hyperelastic_compliance_ALM.jl @@ -71,7 +71,7 @@ function main() dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ ## Finite difference solver and level set function - stencil = AdvectionStencil(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) + ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) ## Setup solver and FE operators state_map = NonlinearFEStateMap(res,U,V,V_φ,U_reg,φh,dΩ,dΓ_N) @@ -83,7 +83,7 @@ function main() vel_ext = VelocityExtension(a_hilb,U_reg,V_reg) ## Optimiser - optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh; + optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh; γ,γ_reinit,verbose=true,constraint_names=[:Vol]) for (it,uh,φh) in optimiser data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh] diff --git a/scripts/Serial/hyperelastic_compliance_neohook_ALM.jl b/scripts/Serial/hyperelastic_compliance_neohook_ALM.jl index c212f0ef..72738c4e 100644 --- a/scripts/Serial/hyperelastic_compliance_neohook_ALM.jl +++ b/scripts/Serial/hyperelastic_compliance_neohook_ALM.jl @@ -85,7 +85,7 @@ function main() dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ ## Finite difference solver and level set function - stencil = AdvectionStencil(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) + ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) ## Setup solver and FE operators state_map = NonlinearFEStateMap(res,U,V,V_φ,U_reg,φh,dΩ,dΓ_N) @@ -97,7 +97,7 @@ function main() vel_ext = VelocityExtension(a_hilb,U_reg,V_reg) ## Optimiser - optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh; + optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh; γ,γ_reinit,verbose=true,constraint_names=[:Vol]) for (it,uh,φh) in optimiser data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh] diff --git a/scripts/Serial/inverse_homenisation_ALM.jl b/scripts/Serial/inverse_homenisation_ALM.jl index 7779b851..0ffdd411 100644 --- a/scripts/Serial/inverse_homenisation_ALM.jl +++ b/scripts/Serial/inverse_homenisation_ALM.jl @@ -71,7 +71,7 @@ function main() dVol(q,u,φ,dΩ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ ## Finite difference solver and level set function - stencil = AdvectionStencil(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) + ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) ## Setup solver and FE operators state_map = RepeatingAffineFEStateMap(3,a,l,U,V,V_φ,U_reg,φh,dΩ) @@ -83,7 +83,7 @@ function main() vel_ext = VelocityExtension(a_hilb,U_reg,V_reg) ## Optimiser - optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh; + optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh; γ,γ_reinit,verbose=true,constraint_names=[:Vol]) for (it,uh,φh) in optimiser data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh))] diff --git a/scripts/Serial/inverter_ALM.jl b/scripts/Serial/inverter_ALM.jl index d1144e2f..132a32e3 100644 --- a/scripts/Serial/inverter_ALM.jl +++ b/scripts/Serial/inverter_ALM.jl @@ -86,7 +86,7 @@ function main() UΓ_out(u,φ,dΩ,dΓ_in,dΓ_out) = ∫((u⋅-e₁-δₓ)/vol_Γ_out)dΓ_out ## Finite difference solver and level set function - stencil = AdvectionStencil(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) + ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) ## Setup solver and FE operators state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_in,dΓ_out) @@ -98,7 +98,7 @@ function main() vel_ext = VelocityExtension(a_hilb,U_reg,V_reg) ## Optimiser - optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh; + optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh; γ,γ_reinit,verbose=true,constraint_names=[:Vol,:UΓ_out]) for (it,uh,φh) in optimiser data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh] diff --git a/scripts/Serial/inverter_HPM.jl b/scripts/Serial/inverter_HPM.jl index ae21a683..39ed22a0 100644 --- a/scripts/Serial/inverter_HPM.jl +++ b/scripts/Serial/inverter_HPM.jl @@ -87,7 +87,7 @@ function main() UΓ_out(u,φ,dΩ,dΓ_in,dΓ_out) = ∫((u⋅-e₁-δₓ)/vol_Γ_out)dΓ_out ## Finite difference solver and level set function - stencil = AdvectionStencil(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) + ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) ## Setup solver and FE operators state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_in,dΓ_out) @@ -100,7 +100,7 @@ function main() ## Optimiser ls_enabled=false # Setting to true will use a line search instead of oscillation detection - optimiser = HilbertianProjection(pcfs,stencil,vel_ext,φh; + optimiser = HilbertianProjection(pcfs,ls_evo,vel_ext,φh; γ,γ_reinit,α_min=0.4,ls_enabled,verbose=true,constraint_names=[:Vol,:UΓ_out]) for (it,uh,φh) in optimiser data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh] diff --git a/scripts/Serial/nonlinear_thermal_compliance_ALM.jl b/scripts/Serial/nonlinear_thermal_compliance_ALM.jl index edbc0b1b..621060a2 100644 --- a/scripts/Serial/nonlinear_thermal_compliance_ALM.jl +++ b/scripts/Serial/nonlinear_thermal_compliance_ALM.jl @@ -74,7 +74,7 @@ function main() dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ ## Finite difference solver and level set function - stencil = AdvectionStencil(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) + ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) ## Setup solver and FE operators state_map = NonlinearFEStateMap(res,U,V,V_φ,U_reg,φh,dΩ,dΓ_N) @@ -86,7 +86,7 @@ function main() vel_ext = VelocityExtension(a_hilb,U_reg,V_reg) ## Optimiser - optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh; + optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh; γ,γ_reinit,verbose=true,constraint_names=[:Vol]) for (it, uh, φh) in optimiser data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh] diff --git a/scripts/Serial/thermal_compliance_ALM.jl b/scripts/Serial/thermal_compliance_ALM.jl index 44903c27..9cd0883e 100644 --- a/scripts/Serial/thermal_compliance_ALM.jl +++ b/scripts/Serial/thermal_compliance_ALM.jl @@ -72,7 +72,7 @@ function main() dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ ## Finite difference solver and level set function - stencil = AdvectionStencil(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) + ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) ## Setup solver and FE operators state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N) @@ -84,7 +84,7 @@ function main() vel_ext = VelocityExtension(a_hilb,U_reg,V_reg) ## Optimiser - optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh; + optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh; γ,γ_reinit,verbose=true,constraint_names=[:Vol]) for (it,uh,φh) in optimiser data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh] diff --git a/scripts/Serial/thermal_compliance_ALM_higherorderlsf.jl b/scripts/Serial/thermal_compliance_ALM_higherorderlsf.jl index 9123a832..1cbc6344 100644 --- a/scripts/Serial/thermal_compliance_ALM_higherorderlsf.jl +++ b/scripts/Serial/thermal_compliance_ALM_higherorderlsf.jl @@ -78,7 +78,7 @@ function main() dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ ## Finite difference solver and level set function - stencil = AdvectionStencil(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) + ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) ## Setup solver and FE operators state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N) @@ -90,7 +90,7 @@ function main() vel_ext = VelocityExtension(a_hilb,U_reg,V_reg) ## Optimiser - optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh; + optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh; γ,γ_reinit,verbose=true,constraint_names=[:Vol]) for (it,uh,φh) in optimiser data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh] diff --git a/scripts/_dev/3d_hyperelastic_compliance_ALM.jl b/scripts/_dev/3d_hyperelastic_compliance_ALM.jl index 2565536b..ddc16243 100644 --- a/scripts/_dev/3d_hyperelastic_compliance_ALM.jl +++ b/scripts/_dev/3d_hyperelastic_compliance_ALM.jl @@ -153,7 +153,7 @@ function main_alt(mesh_partition,distribute,el_size,gz) res(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*((T ∘ ∇(u)) ⊙ ∇(v)))dΩ - ∫(g⋅v)dΓ_N ## Finite difference solver and level set function - # stencil = AdvectionStencil(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps) + # ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps) # reinit!(stencil,φh,γ_reinit) ## Setup solver and FE operators diff --git a/scripts/_dev/Paper_scripts/elast_comp_MPI_3D.jl b/scripts/_dev/Paper_scripts/elast_comp_MPI_3D.jl index 52871792..d177d1da 100644 --- a/scripts/_dev/Paper_scripts/elast_comp_MPI_3D.jl +++ b/scripts/_dev/Paper_scripts/elast_comp_MPI_3D.jl @@ -78,9 +78,9 @@ function main(mesh_partition,distribute) ) # Finite difference scheme scheme = FirstOrderStencil(length(el_size),Float64) - stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps) + ls_evo = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps) # Optimiser - optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit, + optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;γ,γ_reinit, verbose=i_am_main(ranks)) # Solve for (it,uh,φh) in optimiser diff --git a/scripts/_dev/Paper_scripts/hyperelast_comp_MPI_3D.jl b/scripts/_dev/Paper_scripts/hyperelast_comp_MPI_3D.jl index 5b4730c5..bef8a90d 100644 --- a/scripts/_dev/Paper_scripts/hyperelast_comp_MPI_3D.jl +++ b/scripts/_dev/Paper_scripts/hyperelast_comp_MPI_3D.jl @@ -92,9 +92,9 @@ function main(mesh_partition,distribute) ) # Finite difference scheme scheme = FirstOrderStencil(length(el_size),Float64) - stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps) + ls_evo = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps) # Optimiser - optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit, + optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;γ,γ_reinit, verbose=i_am_main(ranks)) # Solve for (it,uh,φh) in optimiser diff --git a/scripts/_dev/Paper_scripts/inverse_hom_MPI_3D.jl b/scripts/_dev/Paper_scripts/inverse_hom_MPI_3D.jl index d15c424c..af7c2735 100644 --- a/scripts/_dev/Paper_scripts/inverse_hom_MPI_3D.jl +++ b/scripts/_dev/Paper_scripts/inverse_hom_MPI_3D.jl @@ -82,9 +82,9 @@ function main(mesh_partition,distribute) ) # Finite difference scheme scheme = FirstOrderStencil(length(el_size),Float64) - stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps) + ls_evo = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps) # Optimiser - optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit, + optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;γ,γ_reinit, verbose=i_am_main(ranks)) # Solve for (it,uh,φh) in optimiser diff --git a/scripts/_dev/Paper_scripts/inverter_MPI_3D.jl b/scripts/_dev/Paper_scripts/inverter_MPI_3D.jl index b9f3269c..d2e25582 100644 --- a/scripts/_dev/Paper_scripts/inverter_MPI_3D.jl +++ b/scripts/_dev/Paper_scripts/inverter_MPI_3D.jl @@ -88,9 +88,9 @@ function main(mesh_partition,distribute) ) # Finite difference scheme scheme = FirstOrderStencil(length(el_size),Float64) - stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps) + ls_evo = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps) # Optimiser - optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit, + optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;γ,γ_reinit, verbose=i_am_main(ranks)) # Solve for (it,uh,φh) in optimiser diff --git a/scripts/_dev/Paper_scripts/therm_comp_MPI.jl b/scripts/_dev/Paper_scripts/therm_comp_MPI.jl index 9b76c455..902742bf 100644 --- a/scripts/_dev/Paper_scripts/therm_comp_MPI.jl +++ b/scripts/_dev/Paper_scripts/therm_comp_MPI.jl @@ -78,9 +78,9 @@ function main(mesh_partition,distribute) ) # Finite difference scheme scheme = FirstOrderStencil(length(el_size),Float64) - stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps) + ls_evo = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps) # Optimiser - optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit,verbose=i_am_main(ranks)) + optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;γ,γ_reinit,verbose=i_am_main(ranks)) # Solve for (it,uh,φh) in optimiser data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh] diff --git a/scripts/_dev/Paper_scripts/therm_comp_MPI_3D.jl b/scripts/_dev/Paper_scripts/therm_comp_MPI_3D.jl index 99be6ec8..5104e0d6 100644 --- a/scripts/_dev/Paper_scripts/therm_comp_MPI_3D.jl +++ b/scripts/_dev/Paper_scripts/therm_comp_MPI_3D.jl @@ -80,9 +80,9 @@ function main(mesh_partition,distribute) ) # Finite difference scheme scheme = FirstOrderStencil(length(el_size),Float64) - stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps) + ls_evo = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps) # Optimiser - optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit, + optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;γ,γ_reinit, verbose=i_am_main(ranks)) # Solve for (it,uh,φh) in optimiser diff --git a/scripts/_dev/Paper_scripts/therm_comp_serial.jl b/scripts/_dev/Paper_scripts/therm_comp_serial.jl index 78cd651c..0b53164f 100644 --- a/scripts/_dev/Paper_scripts/therm_comp_serial.jl +++ b/scripts/_dev/Paper_scripts/therm_comp_serial.jl @@ -63,9 +63,9 @@ function main() vel_ext = VelocityExtension(a_hilb,U_reg,V_reg) # Finite difference scheme scheme = FirstOrderStencil(length(el_size),Float64) - stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps) + ls_evo = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps) # Optimiser - optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit, + optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;γ,γ_reinit, verbose=true) # Solve for (it,uh,φh) in optimiser diff --git a/scripts/_dev/benchmark_tests/MPI_benchmark_test_2D_themal.jl b/scripts/_dev/benchmark_tests/MPI_benchmark_test_2D_themal.jl index e4a309ad..9a0c3e3f 100644 --- a/scripts/_dev/benchmark_tests/MPI_benchmark_test_2D_themal.jl +++ b/scripts/_dev/benchmark_tests/MPI_benchmark_test_2D_themal.jl @@ -62,7 +62,7 @@ function main(mesh_partition,el_size,ranks) dJ(q,u,φ,dΩ,dΓ_N) = ∫((-ξ+D*∇(u)⋅∇(u))*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ; ## Finite difference solver and level set function - stencil = AdvectionStencil(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) + ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) reinit!(stencil,φh,γ_reinit) ## Setup solver and FE operators @@ -75,7 +75,7 @@ function main(mesh_partition,el_size,ranks) vel_ext = VelocityExtension(a_hilb,U_reg,V_reg) ## Optimiser - return AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit) + return AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;γ,γ_reinit) end with_debug() do distribute diff --git a/scripts/_dev/benchmark_tests/MPI_benchmark_test_3D_themal.jl b/scripts/_dev/benchmark_tests/MPI_benchmark_test_3D_themal.jl index 37b149c8..b04bbac8 100644 --- a/scripts/_dev/benchmark_tests/MPI_benchmark_test_3D_themal.jl +++ b/scripts/_dev/benchmark_tests/MPI_benchmark_test_3D_themal.jl @@ -68,7 +68,7 @@ function main(mesh_partition,distribute,el_size) dJ(q,u,φ,dΩ,dΓ_N) = ∫((ξ-D*∇(u)⋅∇(u))*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ ## Finite difference solver and level set function - stencil = AdvectionStencil(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps) + ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps) reinit!(stencil,φh,γ_reinit) ## Setup solver and FE operators @@ -96,7 +96,7 @@ function main(mesh_partition,distribute,el_size) ## Optimiser make_dir(path;ranks=ranks) - return AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit) + return AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;γ,γ_reinit) end with_mpi() do distribute diff --git a/scripts/_dev/benchmark_tests/serial_benchmark_test.jl b/scripts/_dev/benchmark_tests/serial_benchmark_test.jl index 8de90168..17eb71cf 100644 --- a/scripts/_dev/benchmark_tests/serial_benchmark_test.jl +++ b/scripts/_dev/benchmark_tests/serial_benchmark_test.jl @@ -56,7 +56,7 @@ function main() dJ(q,u,φ,dΩ,dΓ_N) = ∫((ξ-D*∇(u)⋅∇(u))*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ ## Finite difference solver and level set function - stencil = AdvectionStencil(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) + ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) reinit!(stencil,φh,γ_reinit) ## Setup solver and FE operators @@ -69,7 +69,7 @@ function main() vel_ext = VelocityExtension(a_hilb,U_reg,V_reg) ## Optimiser - return AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit) + return AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;γ,γ_reinit) end function main_benchmark() diff --git a/scripts/_dev/bug_issue_46/thermal_MPI.jl b/scripts/_dev/bug_issue_46/thermal_MPI.jl index 2e4596ee..f318893d 100644 --- a/scripts/_dev/bug_issue_46/thermal_MPI.jl +++ b/scripts/_dev/bug_issue_46/thermal_MPI.jl @@ -68,9 +68,9 @@ function main(mesh_partition,distribute,use_l::Bool) vel_ext = VelocityExtension(a_hilb, U_reg, V_reg) # Finite difference scheme scheme = FirstOrderStencil(length(el_size),Float64) - stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps) + ls_evo = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps) # Optimiser - optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh; + optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh; γ,γ_reinit,verbose=i_am_main(ranks),constraint_names=[:Vol]) # Solve for (it,uh,φh) in optimiser diff --git a/scripts/_dev/bug_issue_46/thermal_serial.jl b/scripts/_dev/bug_issue_46/thermal_serial.jl index e9c29a6c..2d0e2017 100644 --- a/scripts/_dev/bug_issue_46/thermal_serial.jl +++ b/scripts/_dev/bug_issue_46/thermal_serial.jl @@ -67,9 +67,9 @@ function main(use_l::Bool) vel_ext = VelocityExtension(a_hilb,U_reg,V_reg) # Finite difference scheme scheme = FirstOrderStencil(length(el_size),Float64) - stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps) + ls_evo = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps) # Optimiser - optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh; + optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh; γ,γ_reinit,verbose=true,constraint_names=[:Vol]) # Solve for (it,uh,φh) in optimiser diff --git a/scripts/_dev/full_piezo_script.jl b/scripts/_dev/full_piezo_script.jl index 7cb7c149..7ed4c1ab 100644 --- a/scripts/_dev/full_piezo_script.jl +++ b/scripts/_dev/full_piezo_script.jl @@ -92,7 +92,7 @@ function main(mesh_partition,distribute,el_size) dVol(q,u,φ,dΩ) = ∫(1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ ## Finite difference solver and level set function - stencil = AdvectionStencil(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps) + ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps) reinit!(stencil,φh,γ_reinit) ## Setup solver and FE operators @@ -126,7 +126,7 @@ function main(mesh_partition,distribute,el_size) # ## Optimiser # make_dir(path;ranks=ranks) - # optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit,verbose=i_am_main(ranks)) + # optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;γ,γ_reinit,verbose=i_am_main(ranks)) # for (it, uh, φh) in optimiser # write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh]) # write_history(path*"/history.txt",optimiser.history;ranks=ranks) diff --git a/scripts/_dev/inv_hom_block_assem_testing.jl b/scripts/_dev/inv_hom_block_assem_testing.jl index acd4c897..14d2b24a 100644 --- a/scripts/_dev/inv_hom_block_assem_testing.jl +++ b/scripts/_dev/inv_hom_block_assem_testing.jl @@ -74,7 +74,7 @@ Vol = (u,φ,dΩ) -> ∫(((ρ ∘ φ) - 0.5)/vol_D)dΩ; dVol = (q,u,φ,dΩ) -> ∫(1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ ## Finite difference solver and level set function -stencil = AdvectionStencil(FirstOrderStencil(2,Float64),model,V_φ,el_size./order,max_steps,tol); +ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,el_size./order,max_steps,tol); reinit!(stencil,φ,γ_reinit) ## Initialise op diff --git a/scripts/_dev/inv_hom_block_assem_testing_3D_MPI.jl b/scripts/_dev/inv_hom_block_assem_testing_3D_MPI.jl index 4ab969a6..3c7e57a6 100644 --- a/scripts/_dev/inv_hom_block_assem_testing_3D_MPI.jl +++ b/scripts/_dev/inv_hom_block_assem_testing_3D_MPI.jl @@ -86,7 +86,7 @@ function main(mesh_partition,distribute,el_size,diag_assem::Bool) dVol = (q,u,φ,dΩ) -> ∫(1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ ## Finite difference solver and level set function - stencil = AdvectionStencil(FirstOrderStencil(3,Float64),model,V_φ,el_size./order,max_steps,tol) + ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(3,Float64),model,V_φ,el_size./order,max_steps,tol) reinit!(stencil,φ,γ_reinit) ## Setup solver and FE operators @@ -130,7 +130,7 @@ function main(mesh_partition,distribute,el_size,diag_assem::Bool) ## Optimiser make_dir(path;ranks=ranks) - optimiser = AugmentedLagrangian(φ,pcfs,stencil,vel_ext,interp,el_size,γ,γ_reinit); + optimiser = AugmentedLagrangian(φ,pcfs,ls_evo,vel_ext,interp,el_size,γ,γ_reinit); for history in optimiser it,Ji,Ci,Li = last(history) λi = optimiser.λ; Λi = optimiser.Λ diff --git a/scripts/_dev/inv_hom_block_assem_testing_full_script.jl b/scripts/_dev/inv_hom_block_assem_testing_full_script.jl index f3c31fc2..46d9baf3 100644 --- a/scripts/_dev/inv_hom_block_assem_testing_full_script.jl +++ b/scripts/_dev/inv_hom_block_assem_testing_full_script.jl @@ -77,7 +77,7 @@ function main(diag_block) dVol = (q,u,φ,dΩ) -> ∫(1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ ## Finite difference solver and level set function - stencil = AdvectionStencil(FirstOrderStencil(2,Float64),model,V_φ,el_size./order,max_steps,tol) + ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,el_size./order,max_steps,tol) reinit!(stencil,φ,γ_reinit) ## Setup solver and FE operators @@ -109,7 +109,7 @@ function main(diag_block) ## Optimiser make_dir(path) - optimiser = AugmentedLagrangian(φ,pcfs,stencil,vel_ext,interp,el_size,γ,γ_reinit); + optimiser = AugmentedLagrangian(φ,pcfs,ls_evo,vel_ext,interp,el_size,γ,γ_reinit); for history in optimiser it,Ji,Ci,Li = last(history) λi = optimiser.λ; Λi = optimiser.Λ diff --git a/scripts/_dev/mem_leak.jl b/scripts/_dev/mem_leak.jl index af9a2b28..278ec667 100644 --- a/scripts/_dev/mem_leak.jl +++ b/scripts/_dev/mem_leak.jl @@ -62,9 +62,9 @@ a_hilb(p,q) =∫(α^2*∇(p)⋅∇(q) + p*q)dΩ vel_ext = VelocityExtension(a_hilb,U_reg,V_reg) # Finite difference scheme scheme = FirstOrderStencil(length(el_size),Float64) -stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps) +ls_evo = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps) # Optimiser -optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit,verbose=true,constraint_names=[:Vol]) +optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;γ,γ_reinit,verbose=true,constraint_names=[:Vol]) # Solve for i = 1:10 LevelSetTopOpt.rrule(state_map,φh) diff --git a/scripts/_dev/nonlinear_adjoint_MWE.jl b/scripts/_dev/nonlinear_adjoint_MWE.jl index a84705a2..df5fe18a 100644 --- a/scripts/_dev/nonlinear_adjoint_MWE.jl +++ b/scripts/_dev/nonlinear_adjoint_MWE.jl @@ -70,7 +70,7 @@ function main() J = (u,φ,dΩ,dΓ_N) -> ∫((I ∘ φ)*(D ∘ u)*∇(u)⋅∇(u) + ξ*(ρ ∘ φ))dΩ ## Finite difference solver and level set function - stencil = AdvectionStencil(FirstOrderStencil(2,Float64),model,V_φ,el_size./order,max_steps,tol) + ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,el_size./order,max_steps,tol) reinit!(stencil,φ,γ_reinit) ## Setup solver and FE operators diff --git a/scripts/_dev/tutorial_test/minimum_thermal_compliance_implementation.jl b/scripts/_dev/tutorial_test/minimum_thermal_compliance_implementation.jl index 6867340e..91dbabaa 100644 --- a/scripts/_dev/tutorial_test/minimum_thermal_compliance_implementation.jl +++ b/scripts/_dev/tutorial_test/minimum_thermal_compliance_implementation.jl @@ -61,9 +61,9 @@ a_hilb(p,q) =∫(α^2*∇(p)⋅∇(q) + p*q)dΩ vel_ext = VelocityExtension(a_hilb,U_reg,V_reg) # Finite difference scheme scheme = FirstOrderStencil(length(el_size),Float64) -stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps) +ls_evo = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps) # Optimiser -optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit,verbose=true,constraint_names=[:Vol]) +optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;γ,γ_reinit,verbose=true,constraint_names=[:Vol]) # Solve for (it,uh,φh) in optimiser data = ["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh] diff --git a/scripts/_dev/tutorial_test/minimum_thermal_compliance_implementation_3d.jl b/scripts/_dev/tutorial_test/minimum_thermal_compliance_implementation_3d.jl index 0f167337..8218dabf 100644 --- a/scripts/_dev/tutorial_test/minimum_thermal_compliance_implementation_3d.jl +++ b/scripts/_dev/tutorial_test/minimum_thermal_compliance_implementation_3d.jl @@ -63,9 +63,9 @@ a_hilb(p,q) =∫(α^2*∇(p)⋅∇(q) + p*q)dΩ vel_ext = VelocityExtension(a_hilb,U_reg,V_reg) # Finite difference scheme scheme = FirstOrderStencil(length(el_size),Float64) -stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps) +ls_evo = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps) # Optimiser -optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit,verbose=true,constraint_names=[:Vol]) +optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;γ,γ_reinit,verbose=true,constraint_names=[:Vol]) # Solve for (it,uh,φh) in optimiser data = ["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh] diff --git a/scripts/_dev/tutorial_test/minimum_thermal_compliance_implementation_3d_petsc.jl b/scripts/_dev/tutorial_test/minimum_thermal_compliance_implementation_3d_petsc.jl index a8db1de4..10c1dacf 100644 --- a/scripts/_dev/tutorial_test/minimum_thermal_compliance_implementation_3d_petsc.jl +++ b/scripts/_dev/tutorial_test/minimum_thermal_compliance_implementation_3d_petsc.jl @@ -77,9 +77,9 @@ function main() ) # Finite difference scheme scheme = FirstOrderStencil(length(el_size),Float64) - stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps) + ls_evo = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps) # Optimiser - optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit,verbose=true,constraint_names=[:Vol]) + optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;γ,γ_reinit,verbose=true,constraint_names=[:Vol]) # Solve for (it,uh,φh) in optimiser data = ["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh] diff --git a/scripts/_dev/tutorial_test/minimum_thermal_compliance_implementation_3d_petsc_mpi.jl b/scripts/_dev/tutorial_test/minimum_thermal_compliance_implementation_3d_petsc_mpi.jl index d9fde797..1b88fd91 100644 --- a/scripts/_dev/tutorial_test/minimum_thermal_compliance_implementation_3d_petsc_mpi.jl +++ b/scripts/_dev/tutorial_test/minimum_thermal_compliance_implementation_3d_petsc_mpi.jl @@ -78,9 +78,9 @@ function main(mesh_partition,distribute) ) # Finite difference scheme scheme = FirstOrderStencil(length(el_size),Float64) - stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps) + ls_evo = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps) # Optimiser - optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh; + optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh; γ,γ_reinit,verbose=i_am_main(ranks),constraint_names=[:Vol]) # Solve for (it,uh,φh) in optimiser diff --git a/scripts/_dev/tutorial_test/minimum_thermal_compliance_implementation_nonlinear.jl b/scripts/_dev/tutorial_test/minimum_thermal_compliance_implementation_nonlinear.jl index 058f1f4f..3642fb95 100644 --- a/scripts/_dev/tutorial_test/minimum_thermal_compliance_implementation_nonlinear.jl +++ b/scripts/_dev/tutorial_test/minimum_thermal_compliance_implementation_nonlinear.jl @@ -59,9 +59,9 @@ a_hilb(p,q) =∫(α^2*∇(p)⋅∇(q) + p*q)dΩ vel_ext = VelocityExtension(a_hilb,U_reg,V_reg) # Finite difference scheme scheme = FirstOrderStencil(length(el_size),Float64) -stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps) +ls_evo = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps) # Optimiser -optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit,verbose=true,constraint_names=[:Vol]) +optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;γ,γ_reinit,verbose=true,constraint_names=[:Vol]) # Solve for (it,uh,φh) in optimiser data = ["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh] diff --git a/src/Benchmarks.jl b/src/Benchmarks.jl index 87cc4a6e..0c052a44 100644 --- a/src/Benchmarks.jl +++ b/src/Benchmarks.jl @@ -116,15 +116,15 @@ function benchmark_forward_problem(m::AbstractFEStateMap, φh, ranks; nreps = 10 end """ - benchmark_advection(stencil::AdvectionStencil, φ0, v0, γ, ranks; nreps) + benchmark_advection(stencil::LevelSetEvolution, φ0, v0, γ, ranks; nreps) Benchmark solving the Hamilton-Jacobi evolution equation given a `stencil`, level-set function `φ0`, velocity function `v0`, and time step coefficient `γ`. -See [`advect!`](@ref) for input types. +See [`evolve!`](@ref) for input types. """ -function benchmark_advection(stencil::AdvectionStencil, φ0, v0, γ, ranks; nreps = 10) +function benchmark_advection(stencil::LevelSetEvolution, φ0, v0, γ, ranks; nreps = 10) function f(stencil,φ,v,γ) - advect!(stencil,φ,v,γ) + evolve!(stencil,φ,v,γ) end function reset!(stencil,φ,v,γ) copy!(φ,φ0) @@ -136,12 +136,12 @@ function benchmark_advection(stencil::AdvectionStencil, φ0, v0, γ, ranks; nrep end """ - benchmark_reinitialisation(stencil::AdvectionStencil, φ0, γ_reinit, ranks; nreps) + benchmark_reinitialisation(stencil::LevelSetEvolution, φ0, γ_reinit, ranks; nreps) Benchmark solving the reinitialisation equation given a `stencil`, level-set function `φ0`, and time step coefficient `γ`. See [`reinit!`](@ref) for input types. """ -function benchmark_reinitialisation(stencil::AdvectionStencil, φ0, γ_reinit, ranks; nreps = 10) +function benchmark_reinitialisation(stencil::LevelSetEvolution, φ0, γ_reinit, ranks; nreps = 10) function f(stencil,φ,γ_reinit) reinit!(stencil,φ,γ_reinit) end diff --git a/src/Advection.jl b/src/LevelSetEvolution/HamiltonJacobiEvolution.jl similarity index 51% rename from src/Advection.jl rename to src/LevelSetEvolution/HamiltonJacobiEvolution.jl index 241af032..06dba88b 100644 --- a/src/Advection.jl +++ b/src/LevelSetEvolution/HamiltonJacobiEvolution.jl @@ -1,189 +1,15 @@ """ - abstract type Stencil + struct HamiltonJacobiEvolution{O} -Finite difference stencil for a single step of the Hamilton-Jacobi -evolution equation and reinitialisation equation. - -Your own stencil can be implemented by extending the methods below. -""" -abstract type Stencil end - -""" - allocate_caches(::Stencil,φ,vel) - -Allocate caches for a given `Stencil`. -""" -function allocate_caches(::Stencil,φ,vel) - nothing # By default, no caches are required. -end - -""" - reinit!(::Stencil,φ_new,φ,vel,Δt,Δx,isperiodic,caches) -> φ - -Single finite difference step of the reinitialisation equation for a given `Stencil`. -""" -function reinit!(::Stencil,φ_new,φ,vel,Δt,Δx,isperiodic,caches) - @abstractmethod -end - -""" - advect!(::Stencil,φ,vel,Δt,Δx,isperiodic,caches) -> φ - -Single finite difference step of the Hamilton-Jacobi evoluation equation for a given -`Stencil`. -""" -function advect!(::Stencil,φ,vel,Δt,Δx,isperiodic,caches) - @abstractmethod -end - -""" - compute_Δt(::Stencil,φ,vel) - -Compute the time step for the `Stencil`. -""" -function compute_Δt(::Stencil,φ,vel) - @abstractmethod -end - -""" - struct FirstOrderStencil{D,T} <: Stencil end - -A first order Godunov upwind difference scheme based on Osher and Fedkiw -([link](https://doi.org/10.1007/b98879)). -""" -struct FirstOrderStencil{D,T} <: Stencil - function FirstOrderStencil(D::Integer,::Type{T}) where T<:Real - new{D,T}() - end -end - -function allocate_caches(::FirstOrderStencil{2},φ,vel) - D⁺ʸ = similar(φ); D⁺ˣ = similar(φ) - D⁻ʸ = similar(φ); D⁻ˣ = similar(φ) - ∇⁺ = similar(φ); ∇⁻ = similar(φ) - return D⁺ʸ, D⁺ˣ, D⁻ʸ, D⁻ˣ, ∇⁺, ∇⁻ -end - -function reinit!(::FirstOrderStencil{2,T},φ_new,φ,vel,Δt,Δx,isperiodic,caches) where T - D⁺ʸ, D⁺ˣ, D⁻ʸ, D⁻ˣ, ∇⁺, ∇⁻ = caches - Δx, Δy = Δx - xperiodic,yperiodic = isperiodic - # Prepare shifted lsf - circshift!(D⁺ʸ,φ,(0,-1)); circshift!(D⁻ʸ,φ,(0,1)) - circshift!(D⁺ˣ,φ,(-1,0)); circshift!(D⁻ˣ,φ,(1,0)) - # Sign approximation - ∇⁺ .= @. (D⁺ʸ - D⁻ʸ)/(2Δy); ~yperiodic ? ∇⁺[:,[1,end]] .= zero(T) : 0; - ∇⁻ .= @. (D⁺ˣ - D⁻ˣ)/(2Δx); ~xperiodic ? ∇⁻[[1,end],:] .= zero(T) : 0; - ϵₛ = minimum((Δx,Δy)) - vel .= @. φ/sqrt(φ^2 + ϵₛ^2*(∇⁺^2+∇⁻^2)) - # Forward (+) & Backward (-) - D⁺ʸ .= @. (D⁺ʸ - φ)/Δy; ~yperiodic ? D⁺ʸ[:,end] .= zero(T) : 0; - D⁺ˣ .= @. (D⁺ˣ - φ)/Δx; ~xperiodic ? D⁺ˣ[end,:] .= zero(T) : 0; - D⁻ʸ .= @. (φ - D⁻ʸ)/Δy; ~yperiodic ? D⁻ʸ[:,1] .= zero(T) : 0; - D⁻ˣ .= @. (φ - D⁻ˣ)/Δx; ~xperiodic ? D⁻ˣ[1,:] .= zero(T) : 0; - # Operators - ∇⁺ .= @. sqrt(max(D⁻ʸ,0)^2 + min(D⁺ʸ,0)^2 + max(D⁻ˣ,0)^2 + min(D⁺ˣ,0)^2); - ∇⁻ .= @. sqrt(max(D⁺ʸ,0)^2 + min(D⁻ʸ,0)^2 + max(D⁺ˣ,0)^2 + min(D⁻ˣ,0)^2); - # Update - φ_new .= @. φ - Δt*(max(vel,0)*∇⁺ + min(vel,0)*∇⁻ - vel) - return φ_new -end - -function advect!(::FirstOrderStencil{2,T},φ,vel,Δt,Δx,isperiodic,caches) where T - D⁺ʸ, D⁺ˣ, D⁻ʸ, D⁻ˣ, ∇⁺, ∇⁻ = caches - Δx, Δy = Δx - xperiodic,yperiodic = isperiodic - # Prepare shifted lsf - circshift!(D⁺ʸ,φ,(0,-1)); circshift!(D⁻ʸ,φ,(0,1)) - circshift!(D⁺ˣ,φ,(-1,0)); circshift!(D⁻ˣ,φ,(1,0)) - # Forward (+) & Backward (-) - D⁺ʸ .= @. (D⁺ʸ - φ)/Δy; ~yperiodic ? D⁺ʸ[:,end] .= zero(T) : 0; - D⁺ˣ .= @. (D⁺ˣ - φ)/Δx; ~xperiodic ? D⁺ˣ[end,:] .= zero(T) : 0; - D⁻ʸ .= @. (φ - D⁻ʸ)/Δy; ~yperiodic ? D⁻ʸ[:,1] .= zero(T) : 0; - D⁻ˣ .= @. (φ - D⁻ˣ)/Δx; ~xperiodic ? D⁻ˣ[1,:] .= zero(T) : 0; - # Operators - ∇⁺ .= @. sqrt(max(D⁻ʸ,0)^2 + min(D⁺ʸ,0)^2 + max(D⁻ˣ,0)^2 + min(D⁺ˣ,0)^2) - ∇⁻ .= @. sqrt(max(D⁺ʸ,0)^2 + min(D⁻ʸ,0)^2 + max(D⁺ˣ,0)^2 + min(D⁻ˣ,0)^2) - # Update - φ .= @. φ - Δt*(max(vel,0)*∇⁺ + min(vel,0)*∇⁻) - return φ -end - -function allocate_caches(::FirstOrderStencil{3},φ,vel) - D⁺ᶻ = similar(φ); D⁺ʸ = similar(φ); D⁺ˣ = similar(φ) - D⁻ᶻ = similar(φ); D⁻ʸ = similar(φ); D⁻ˣ = similar(φ) - ∇⁺ = similar(φ); ∇⁻ = similar(φ) - return D⁺ᶻ, D⁺ʸ, D⁺ˣ, D⁻ᶻ, D⁻ʸ, D⁻ˣ, ∇⁺, ∇⁻ -end - -function reinit!(::FirstOrderStencil{3,T},φ_new,φ,vel,Δt,Δx,isperiodic,caches) where T - D⁺ᶻ, D⁺ʸ, D⁺ˣ, D⁻ᶻ, D⁻ʸ, D⁻ˣ, ∇⁺, ∇⁻=caches - Δx, Δy, Δz = Δx - xperiodic,yperiodic,zperiodic = isperiodic - # Prepare shifted lsf - circshift!(D⁺ʸ,φ,(0,-1,0)); circshift!(D⁻ʸ,φ,(0,1,0)) - circshift!(D⁺ˣ,φ,(-1,0,0)); circshift!(D⁻ˣ,φ,(1,0,0)) - circshift!(D⁺ᶻ,φ,(0,0,-1)); circshift!(D⁻ᶻ,φ,(0,0,1)) - # Sign approximation - ∇⁺ .= @. (D⁺ʸ - D⁻ʸ)/(2Δy); ~yperiodic ? ∇⁺[:,[1,end],:] .= zero(T) : 0; - ∇⁻ .= @. (D⁺ˣ - D⁻ˣ)/(2Δx); ~xperiodic ? ∇⁻[[1,end],:,:] .= zero(T) : 0; - ∇⁺ .= @. ∇⁺^2+∇⁻^2 # |∇φ|² (partially computed) - ∇⁻ .= @. (D⁺ᶻ - D⁻ᶻ)/(2Δz); ~xperiodic ? ∇⁻[:,:,[1,end]] .= zero(T) : 0; - ϵₛ = minimum((Δx,Δy,Δz)) - vel .= @. φ/sqrt(φ^2 + ϵₛ^2*(∇⁺+∇⁻^2)) - # Forward (+) & Backward (-) - D⁺ʸ .= (D⁺ʸ - φ)/Δy; ~yperiodic ? D⁺ʸ[:,end,:] .= zero(T) : 0; - D⁺ˣ .= (D⁺ˣ - φ)/Δx; ~xperiodic ? D⁺ˣ[end,:,:] .= zero(T) : 0; - D⁺ᶻ .= (D⁺ᶻ - φ)/Δz; ~zperiodic ? D⁺ᶻ[:,:,end] .= zero(T) : 0; - D⁻ʸ .= (φ - D⁻ʸ)/Δy; ~yperiodic ? D⁻ʸ[:,1,:] .= zero(T) : 0; - D⁻ˣ .= (φ - D⁻ˣ)/Δx; ~xperiodic ? D⁻ˣ[1,:,:] .= zero(T) : 0; - D⁻ᶻ .= (φ - D⁻ᶻ)/Δz; ~zperiodic ? D⁻ᶻ[:,:,1] .= zero(T) : 0; - # Operators - ∇⁺ .= @. sqrt(max(D⁻ʸ,0)^2 + min(D⁺ʸ,0)^2 + max(D⁻ˣ,0)^2 + min(D⁺ˣ,0)^2 + max(D⁻ᶻ,0)^2 + min(D⁺ᶻ,0)^2) - ∇⁻ .= @. sqrt(max(D⁺ʸ,0)^2 + min(D⁻ʸ,0)^2 + max(D⁺ˣ,0)^2 + min(D⁻ˣ,0)^2 + max(D⁺ᶻ,0)^2 + min(D⁻ᶻ,0)^2) - # Update - φ_new .= @. φ - Δt*(max(vel,0)*∇⁺ + min(vel,0)*∇⁻ - vel) - return φ_new -end - -function advect!(::FirstOrderStencil{3,T},φ,vel,Δt,Δx,isperiodic,caches) where T - D⁺ᶻ, D⁺ʸ, D⁺ˣ, D⁻ᶻ, D⁻ʸ, D⁻ˣ, ∇⁺, ∇⁻=caches - Δx, Δy, Δz = Δx - xperiodic,yperiodic,zperiodic = isperiodic - # Prepare shifted lsf - circshift!(D⁺ʸ,φ,(0,-1,0)); circshift!(D⁻ʸ,φ,(0,1,0)) - circshift!(D⁺ˣ,φ,(-1,0,0)); circshift!(D⁻ˣ,φ,(1,0,0)) - circshift!(D⁺ᶻ,φ,(0,0,-1)); circshift!(D⁻ᶻ,φ,(0,0,1)) - # Forward (+) & Backward (-) - D⁺ʸ .= (D⁺ʸ - φ)/Δy; ~yperiodic ? D⁺ʸ[:,end,:] .= zero(T) : 0; - D⁺ˣ .= (D⁺ˣ - φ)/Δx; ~xperiodic ? D⁺ˣ[end,:,:] .= zero(T) : 0; - D⁺ᶻ .= (D⁺ᶻ - φ)/Δz; ~zperiodic ? D⁺ᶻ[:,:,end] .= zero(T) : 0; - D⁻ʸ .= (φ - D⁻ʸ)/Δy; ~yperiodic ? D⁻ʸ[:,1,:] .= zero(T) : 0; - D⁻ˣ .= (φ - D⁻ˣ)/Δx; ~xperiodic ? D⁻ˣ[1,:,:] .= zero(T) : 0; - D⁻ᶻ .= (φ - D⁻ᶻ)/Δz; ~zperiodic ? D⁻ᶻ[:,:,1] .= zero(T) : 0; - # Operators - ∇⁺ .= @. sqrt(max(D⁻ʸ,0)^2 + min(D⁺ʸ,0)^2 + max(D⁻ˣ,0)^2 + min(D⁺ˣ,0)^2 + max(D⁻ᶻ,0)^2 + min(D⁺ᶻ,0)^2) - ∇⁻ .= @. sqrt(max(D⁺ʸ,0)^2 + min(D⁻ʸ,0)^2 + max(D⁺ˣ,0)^2 + min(D⁻ˣ,0)^2 + max(D⁺ᶻ,0)^2 + min(D⁻ᶻ,0)^2) - # Update - φ .= @. φ - Δt*(max(vel,0)*∇⁺ + min(vel,0)*∇⁻) - return φ -end - -function compute_Δt(::FirstOrderStencil{D,T},Δ,γ,φ,vel) where {D,T} - v_norm = maximum(abs,vel) - return γ * min(Δ...) / (eps(T)^2 + v_norm) -end - -""" - struct AdvectionStencil{O} - -Wrapper to enable finite differencing to solve the +A standard forward Euler in time finite difference method for solving the Hamilton-Jacobi evolution equation and reinitialisation equation -on order `O` finite elements in serial or parallel. +on order `O` finite elements in serial or parallel. + +Based on the scheme by Osher and Fedkiw ([link](https://doi.org/10.1007/b98879)). # Parameters -- `stencil::Stencil`: Finite difference stencil for a single step HJ +- `stencil::Stencil`: Spatial finite difference stencil for a single step HJ equation and reinitialisation equation. - `model`: A `CartesianDiscreteModel`. - `space`: FE space for level-set function @@ -191,7 +17,7 @@ on order `O` finite elements in serial or parallel. - `params`: Tuple of additional params - `cache`: Stencil cache """ -struct AdvectionStencil{O} +struct HamiltonJacobiEvolution{O} <: LevelSetEvolution stencil :: Stencil model space @@ -201,13 +27,13 @@ struct AdvectionStencil{O} end """ - AdvectionStencil(stencil::Stencil,model,space,tol,max_steps,max_steps_reinit) + HamiltonJacobiEvolution(stencil::Stencil,model,space,tol,max_steps,max_steps_reinit) -Create an instance of `AdvectionStencil` given a stencil, model, FE space, and +Create an instance of `HamiltonJacobiEvolution` given a stencil, model, FE space, and additional optional arguments. This automatically creates the DoF permutation to handle high-order finite elements. """ -function AdvectionStencil( +function HamiltonJacobiEvolution( stencil::Stencil, model, space, @@ -219,6 +45,9 @@ function AdvectionStencil( order, isperiodic, Δ, ndof = get_stencil_params(model,space) params = (;isperiodic,Δ,ndof,max_steps,max_steps_reinit,tol) + # Check that we have sufficient order + check_order(stencil,order) + # Dof permutation perm = create_dof_permutation(model,space,order) @@ -226,139 +55,29 @@ function AdvectionStencil( φ, vel = zero_free_values(space), zero_free_values(space) cache = allocate_caches(stencil,φ,vel,perm,order,ndof) - return AdvectionStencil{order}(stencil,model,space,perm,params,cache) -end - -function get_stencil_params(model::CartesianDiscreteModel,space::FESpace) - order = get_order(first(Gridap.CellData.get_data(get_fe_basis(space)))) - desc = get_cartesian_descriptor(model) - isperiodic = desc.isperiodic - ndof = order .* desc.partition .+ 1 .- isperiodic - Δ = desc.sizes ./ order - return order, isperiodic, Δ, ndof -end - -function get_stencil_params(model::DistributedDiscreteModel,space::DistributedFESpace) - order, isperiodic, Δ, ndof = map(local_views(model),local_views(space)) do model, space - get_stencil_params(model,space) - end |> PartitionedArrays.tuple_of_arrays - - isperiodic = getany(isperiodic) - order = getany(order) - Δ = getany(Δ) - return order, isperiodic, Δ, ndof -end - -Gridap.ReferenceFEs.get_order(f::Gridap.Fields.LinearCombinationFieldVector) = get_order(f.fields) - -# Create dof permutation vector to enable finite differences on -# higher order Lagrangian finite elements on a Cartesian mesh. -function create_dof_permutation(model::CartesianDiscreteModel{Dc}, - space::UnconstrainedFESpace, - order::Integer) where Dc - function get_terms(poly::Polytope, orders) - _nodes, facenodes = Gridap.ReferenceFEs._compute_nodes(poly, orders) - terms = Gridap.ReferenceFEs._coords_to_terms(_nodes, orders) - return terms - end - desc = get_cartesian_descriptor(model) - - periodic = desc.isperiodic - ncells = desc.partition - ndofs = order .* ncells .+ 1 .- periodic - @check prod(ndofs) == num_free_dofs(space) - - new_dof_ids = CircularArray(LinearIndices(ndofs)) - n2o_dof_map = fill(-1,num_free_dofs(space)) - - terms = get_terms(first(get_polytopes(model)), fill(order,Dc)) - cell_dof_ids = get_cell_dof_ids(space) - cache_cell_dof_ids = array_cache(cell_dof_ids) - for (iC,cell) in enumerate(CartesianIndices(ncells)) - first_new_dof = order .* (Tuple(cell) .- 1) .+ 1 - new_dofs_range = map(i -> i:i+order,first_new_dof) - new_dofs = view(new_dof_ids,new_dofs_range...) - - cell_dofs = getindex!(cache_cell_dof_ids,cell_dof_ids,iC) - for (iDof, dof) in enumerate(cell_dofs) - t = terms[iDof] - #o2n_dof_map[dof] = new_dofs[t] - n2o_dof_map[new_dofs[t]] = dof - end - end - - return n2o_dof_map -end - -function create_dof_permutation(model::GridapDistributed.DistributedDiscreteModel, - space::GridapDistributed.DistributedFESpace, - order::Integer) - local_perms = map(local_views(model),local_views(space)) do model, space - create_dof_permutation(model,space,order) - end - return local_perms -end - -function PartitionedArrays.permute_indices(indices::LocalIndices,perm) - id = part_id(indices) - n_glob = global_length(indices) - l2g = view(local_to_global(indices),perm) - l2o = view(local_to_owner(indices),perm) - return LocalIndices(n_glob,id,l2g,l2o) -end - -function allocate_caches(s::Stencil,φ::Vector,vel::Vector,perm,order,ndofs) - stencil_caches = allocate_caches(s,reshape(φ,ndofs),reshape(vel,ndofs)) - φ_tmp = similar(φ) - vel_tmp = similar(vel) - perm_caches = (order >= 2) ? (similar(φ), similar(vel)) : nothing - return φ_tmp, vel_tmp, perm_caches, stencil_caches -end - -function allocate_caches(s::Stencil,φ::PVector,vel::PVector,perm,order,local_ndofs) - local_stencil_caches = map(local_views(φ),local_views(vel),local_views(local_ndofs)) do φ,vel,ndofs - allocate_caches(s,reshape(φ,ndofs),reshape(vel,ndofs)) - end - - perm_indices = map(permute_indices,partition(axes(φ,1)),perm) - perm_caches = (order >= 2) ? (pfill(0.0,perm_indices),pfill(0.0,perm_indices)) : nothing - - φ_tmp = (order >= 2) ? pfill(0.0,perm_indices) : similar(φ) - vel_tmp = (order >= 2) ? pfill(0.0,perm_indices) : similar(vel) - return φ_tmp, vel_tmp, perm_caches, local_stencil_caches + return HamiltonJacobiEvolution{order}(stencil,model,space,perm,params,cache) end -function permute!(x_out,x_in,perm) - for (i_new,i_old) in enumerate(perm) - x_out[i_new] = x_in[i_old] - end - return x_out -end - -function permute!(x_out::PVector,x_in::PVector,perm) - map(permute!,partition(x_out),partition(x_in),perm) - return x_out -end +""" + get_dof_Δ(m::HamiltonJacobiEvolution) -function permute_inv!(x_out,x_in,perm) - for (i_new,i_old) in enumerate(perm) - x_out[i_old] = x_in[i_new] - end - return x_out -end -function permute_inv!(x_out::PVector,x_in::PVector,perm) - map(permute_inv!,partition(x_out),partition(x_in),perm) - return x_out +Return the distance betweem degree of freedom +""" +function get_dof_Δ(m::HamiltonJacobiEvolution) + return m.params.Δ end -function advect!(s::AdvectionStencil,φh,args...) - advect!(s,get_free_dof_values(φh),args...) +# Compute the time step for the `HamiltonJacobiEvolution`. +function compute_Δt(::HamiltonJacobiEvolution,Δ,γ,φ,vel) + T = eltype(γ) + v_norm = maximum(abs,vel) + return γ * min(Δ...) / (eps(T)^2 + v_norm) end """ - advect!(s::AdvectionStencil{O},φ,vel,γ) where O + evolve!(s::HamiltonJacobiEvolution{O},φ,vel,γ) where O -Solve the Hamilton-Jacobi evolution equation using the `AdvectionStencil`. +Solve the Hamilton-Jacobi evolution equation using the `HamiltonJacobiEvolution`. # Hamilton-Jacobi evolution equation ``\\frac{\\partial\\phi}{\\partial t} + V(\\boldsymbol{x})\\lVert\\boldsymbol{\\nabla}\\phi\\rVert = 0,`` @@ -367,12 +86,12 @@ with ``\\phi(0,\\boldsymbol{x})=\\phi_0(\\boldsymbol{x})`` and ``\\boldsymbol{x} # Arguments -- `s::AdvectionStencil{O}`: Stencil for computation +- `s::HamiltonJacobiEvolution{O}`: Method - `φ`: level set function as a vector of degrees of freedom - `vel`: the normal velocity as a vector of degrees of freedom - `γ`: coeffient on the time step size. """ -function advect!(s::AdvectionStencil{O},φ::PVector,vel::PVector,γ) where O +function evolve!(s::HamiltonJacobiEvolution{O},φ::PVector,vel::PVector,γ) where O _, _, perm_caches, stencil_cache = s.cache Δ, isperiodic, = s.params.Δ, s.params.isperiodic ndof, max_steps = s.params.ndof, s.params.max_steps @@ -381,13 +100,13 @@ function advect!(s::AdvectionStencil{O},φ::PVector,vel::PVector,γ) where O _vel = (O >= 2) ? permute!(perm_caches[2],vel,s.perm) : vel ## CFL Condition (requires γ≤1.0) - Δt = compute_Δt(s.stencil,Δ,γ,φ,vel) + Δt = compute_Δt(s,Δ,γ,φ,vel) for _ in 1:max_steps # Apply operations across partitions map(local_views(_φ),local_views(_vel),stencil_cache,ndof) do _φ,_vel,stencil_cache,S φ_mat = reshape(_φ,S) vel_mat = reshape(_vel,S) - advect!(s.stencil,φ_mat,vel_mat,Δt,Δ,isperiodic,stencil_cache) + evolve!(s.stencil,φ_mat,vel_mat,Δt,Δ,isperiodic,stencil_cache) end # Update ghost nodes consistent!(_φ) |> fetch @@ -396,7 +115,7 @@ function advect!(s::AdvectionStencil{O},φ::PVector,vel::PVector,γ) where O return φ end -function advect!(s::AdvectionStencil{O},φ::Vector,vel::Vector,γ) where O +function evolve!(s::HamiltonJacobiEvolution{O},φ::Vector,vel::Vector,γ) where O _, _, perm_caches, stencil_cache = s.cache Δ, isperiodic, = s.params.Δ, s.params.isperiodic ndof, max_steps = s.params.ndof, s.params.max_steps @@ -405,24 +124,24 @@ function advect!(s::AdvectionStencil{O},φ::Vector,vel::Vector,γ) where O _vel = (O >= 2) ? permute!(perm_caches[2],vel,s.perm) : vel ## CFL Condition (requires γ≤1.0) - Δt = compute_Δt(s.stencil,Δ,γ,φ,vel) + Δt = compute_Δt(s,Δ,γ,φ,vel) for _ in 1:max_steps φ_mat = reshape(_φ,ndof) vel_mat = reshape(_vel,ndof) - advect!(s.stencil,φ_mat,vel_mat,Δt,Δ,isperiodic,stencil_cache) + evolve!(s.stencil,φ_mat,vel_mat,Δt,Δ,isperiodic,stencil_cache) end φ = (O >= 2) ? permute_inv!(φ,_φ,s.perm) : _φ return φ end -function reinit!(s::AdvectionStencil,φh,args...) - reinit!(s,get_free_dof_values(φh),args...) +function evolve!(s::HamiltonJacobiEvolution,φh,args...) + evolve!(s,get_free_dof_values(φh),args...) end """ - reinit!(s::AdvectionStencil{O},φ,γ) where O + reinit!(s::HamiltonJacobiEvolution{O},φ,γ) where O -Solve the reinitialisation equation using the `AdvectionStencil`. +Solve the reinitialisation equation using the `HamiltonJacobiEvolution`. # Reinitialisation equation ``\\frac{\\partial\\phi}{\\partial t} + \\mathrm{sign}(\\phi_0)(\\lVert\\boldsymbol{\\nabla}\\phi\\rVert-1) = 0,`` @@ -431,19 +150,19 @@ with ``\\phi(0,\\boldsymbol{x})=\\phi_0(\\boldsymbol{x})`` and ``\\boldsymbol{x} # Arguments -- `s::AdvectionStencil{O}`: Stencil for computation +- `s::HamiltonJacobiEvolution{O}`: Method - `φ`: level set function as a vector of degrees of freedom - `γ`: coeffient on the time step size. """ -function reinit!(s::AdvectionStencil{O},φ::PVector,γ) where O +function reinit!(s::HamiltonJacobiEvolution{O},φ::PVector,γ) where O φ_tmp, vel_tmp, perm_caches, stencil_cache = s.cache Δ, isperiodic, ndof = s.params.Δ, s.params.isperiodic, s.params.ndof tol, max_steps = s.params.tol, s.params.max_steps_reinit _φ = (O >= 2) ? permute!(perm_caches[1],φ,s.perm) : φ - ## CFL Condition (requires γ≤0.5). Note inform(vel_tmp) = 1.0 - Δt = compute_Δt(s.stencil,Δ,γ,_φ,1.0) + ## CFL Condition (requires γ≤0.5). Note infnorm(S) = 1.0 + Δt = compute_Δt(s,Δ,γ,_φ,1.0) # Apply operations across partitions step = 1; err = maximum(abs,φ); fill!(φ_tmp,0.0) @@ -469,7 +188,7 @@ function reinit!(s::AdvectionStencil{O},φ::PVector,γ) where O return φ end -function reinit!(s::AdvectionStencil{O},φ::Vector,γ) where O +function reinit!(s::HamiltonJacobiEvolution{O},φ::Vector,γ) where O φ_tmp, vel_tmp, perm_caches, stencil_cache = s.cache Δ, isperiodic, ndof = s.params.Δ, s.params.isperiodic, s.params.ndof tol, max_steps = s.params.tol, s.params.max_steps_reinit @@ -477,7 +196,7 @@ function reinit!(s::AdvectionStencil{O},φ::Vector,γ) where O _φ = (O >= 2) ? permute!(perm_caches[1],φ,s.perm) : φ ## CFL Condition (requires γ≤0.5) - Δt = compute_Δt(s.stencil,Δ,γ,_φ,1.0) + Δt = compute_Δt(s,Δ,γ,_φ,1.0) # Apply operations across partitions step = 1; err = maximum(abs,φ); fill!(φ_tmp,0.0) @@ -499,3 +218,130 @@ function reinit!(s::AdvectionStencil{O},φ::Vector,γ) where O φ = (O >= 2) ? permute_inv!(φ,_φ,s.perm) : _φ return φ end + +function reinit!(s::HamiltonJacobiEvolution,φh,args...) + reinit!(s,get_free_dof_values(φh),args...) +end + +## Utilities +function get_stencil_params(model::CartesianDiscreteModel,space::FESpace) + order = get_order(first(Gridap.CellData.get_data(get_fe_basis(space)))) + desc = get_cartesian_descriptor(model) + isperiodic = desc.isperiodic + ndof = order .* desc.partition .+ 1 .- isperiodic + Δ = desc.sizes ./ order + return order, isperiodic, Δ, ndof +end + +function get_stencil_params(model::DistributedDiscreteModel,space::DistributedFESpace) + order, isperiodic, Δ, ndof = map(local_views(model),local_views(space)) do model, space + get_stencil_params(model,space) + end |> PartitionedArrays.tuple_of_arrays + + isperiodic = getany(isperiodic) + order = getany(order) + Δ = getany(Δ) + return order, isperiodic, Δ, ndof +end + +Gridap.ReferenceFEs.get_order(f::Gridap.Fields.LinearCombinationFieldVector) = get_order(f.fields) + +# Create dof permutation vector to enable finite differences on +# higher order Lagrangian finite elements on a Cartesian mesh. +function create_dof_permutation(model::CartesianDiscreteModel{Dc}, + space::UnconstrainedFESpace, + order::Integer) where Dc + function get_terms(poly::Polytope, orders) + _nodes, facenodes = Gridap.ReferenceFEs._compute_nodes(poly, orders) + terms = Gridap.ReferenceFEs._coords_to_terms(_nodes, orders) + return terms + end + desc = get_cartesian_descriptor(model) + + periodic = desc.isperiodic + ncells = desc.partition + ndofs = order .* ncells .+ 1 .- periodic + @check prod(ndofs) == num_free_dofs(space) + + new_dof_ids = CircularArray(LinearIndices(ndofs)) + n2o_dof_map = fill(-1,num_free_dofs(space)) + + terms = get_terms(first(get_polytopes(model)), fill(order,Dc)) + cell_dof_ids = get_cell_dof_ids(space) + cache_cell_dof_ids = array_cache(cell_dof_ids) + for (iC,cell) in enumerate(CartesianIndices(ncells)) + first_new_dof = order .* (Tuple(cell) .- 1) .+ 1 + new_dofs_range = map(i -> i:i+order,first_new_dof) + new_dofs = view(new_dof_ids,new_dofs_range...) + + cell_dofs = getindex!(cache_cell_dof_ids,cell_dof_ids,iC) + for (iDof, dof) in enumerate(cell_dofs) + t = terms[iDof] + #o2n_dof_map[dof] = new_dofs[t] + n2o_dof_map[new_dofs[t]] = dof + end + end + + return n2o_dof_map +end + +function create_dof_permutation(model::GridapDistributed.DistributedDiscreteModel, + space::GridapDistributed.DistributedFESpace, + order::Integer) + local_perms = map(local_views(model),local_views(space)) do model, space + create_dof_permutation(model,space,order) + end + return local_perms +end + +function PartitionedArrays.permute_indices(indices::LocalIndices,perm) + id = part_id(indices) + n_glob = global_length(indices) + l2g = view(local_to_global(indices),perm) + l2o = view(local_to_owner(indices),perm) + return LocalIndices(n_glob,id,l2g,l2o) +end + +function allocate_caches(s::Stencil,φ::Vector,vel::Vector,perm,order,ndofs) + stencil_caches = allocate_caches(s,reshape(φ,ndofs),reshape(vel,ndofs)) + φ_tmp = similar(φ) + vel_tmp = similar(vel) + perm_caches = (order >= 2) ? (similar(φ), similar(vel)) : nothing + return φ_tmp, vel_tmp, perm_caches, stencil_caches +end + +function allocate_caches(s::Stencil,φ::PVector,vel::PVector,perm,order,local_ndofs) + local_stencil_caches = map(local_views(φ),local_views(vel),local_views(local_ndofs)) do φ,vel,ndofs + allocate_caches(s,reshape(φ,ndofs),reshape(vel,ndofs)) + end + + perm_indices = map(permute_indices,partition(axes(φ,1)),perm) + perm_caches = (order >= 2) ? (pfill(0.0,perm_indices),pfill(0.0,perm_indices)) : nothing + + φ_tmp = (order >= 2) ? pfill(0.0,perm_indices) : similar(φ) + vel_tmp = (order >= 2) ? pfill(0.0,perm_indices) : similar(vel) + return φ_tmp, vel_tmp, perm_caches, local_stencil_caches +end + +function permute!(x_out,x_in,perm) + for (i_new,i_old) in enumerate(perm) + x_out[i_new] = x_in[i_old] + end + return x_out +end + +function permute!(x_out::PVector,x_in::PVector,perm) + map(permute!,partition(x_out),partition(x_in),perm) + return x_out +end + +function permute_inv!(x_out,x_in,perm) + for (i_new,i_old) in enumerate(perm) + x_out[i_old] = x_in[i_new] + end + return x_out +end +function permute_inv!(x_out::PVector,x_in::PVector,perm) + map(permute_inv!,partition(x_out),partition(x_in),perm) + return x_out +end \ No newline at end of file diff --git a/src/LevelSetEvolution/LevelSetEvolution.jl b/src/LevelSetEvolution/LevelSetEvolution.jl new file mode 100644 index 00000000..6f9c4f32 --- /dev/null +++ b/src/LevelSetEvolution/LevelSetEvolution.jl @@ -0,0 +1,39 @@ +""" + abstract type LevelSetEvolution + +Your own evolution method can be implemented by implementing +concrete functionality of the below. +""" +abstract type LevelSetEvolution end + +""" + evolve!(::LevelSetEvolution,φ,args...) + +Evolve the level set function φ according to an evolution +method LevelSetEvolution. +""" +function evolve!(::LevelSetEvolution,φ,args...) + @abstractmethod +end + +""" + reinit!(::LevelSetEvolution,φ,args...) + +Reinitialise the level set function φ according to an +evolution method LevelSetEvolution. +""" +function reinit!(::LevelSetEvolution,φ,args...) + @abstractmethod +end + +""" + get_dof_Δ(::LevelSetEvolution) + +Return the distance betweem degree of freedom +""" +function get_dof_Δ(::LevelSetEvolution) + @abstractmethod +end + +include("Stencil.jl") +include("HamiltonJacobiEvolution.jl") \ No newline at end of file diff --git a/src/LevelSetEvolution/Stencil.jl b/src/LevelSetEvolution/Stencil.jl new file mode 100644 index 00000000..4b15e82e --- /dev/null +++ b/src/LevelSetEvolution/Stencil.jl @@ -0,0 +1,175 @@ +""" + abstract type Stencil + +Finite difference stencil for a single step of the Hamilton-Jacobi +evolution equation and reinitialisation equation. + +Your own spatial stencil can be implemented by extending the methods below. +""" +abstract type Stencil end + +""" + allocate_caches(::Stencil,φ,vel) + +Allocate caches for a given `Stencil`. +""" +function allocate_caches(::Stencil,φ,vel) + nothing # By default, no caches are required. +end + +""" + check_order(::Stencil,order) + +Throw error if insufficient reference element order +to implement stencil in parallel. +""" +function check_order(::Stencil,order) + @abstractmethod +end + +""" + reinit!(::Stencil,φ_new,φ,vel,Δt,Δx,isperiodic,caches) -> φ + +Single finite difference step of the reinitialisation equation for a given `Stencil`. +""" +function reinit!(::Stencil,φ_new,φ,vel,Δt,Δx,isperiodic,caches) + @abstractmethod +end + +""" + evolve!(::Stencil,φ,vel,Δt,Δx,isperiodic,caches) -> φ + +Single finite difference step of the Hamilton-Jacobi evoluation equation for a given +`Stencil`. +""" +function evolve!(::Stencil,φ,vel,Δt,Δx,isperiodic,caches) + @abstractmethod +end + +""" + struct FirstOrderStencil{D,T} <: Stencil end + +A first order upwind difference scheme based on Osher and Fedkiw +([link](https://doi.org/10.1007/b98879)). +""" +struct FirstOrderStencil{D,T} <: Stencil + function FirstOrderStencil(D::Integer,::Type{T}) where T<:Real + new{D,T}() + end +end + +function check_order(::FirstOrderStencil,order) + @check order >= 1 "FirstOrderStencil requires reference element to have order >= 1" +end + +function allocate_caches(::FirstOrderStencil{2},φ,vel) + D⁺ʸ = similar(φ); D⁺ˣ = similar(φ) + D⁻ʸ = similar(φ); D⁻ˣ = similar(φ) + ∇⁺ = similar(φ); ∇⁻ = similar(φ) + return D⁺ʸ, D⁺ˣ, D⁻ʸ, D⁻ˣ, ∇⁺, ∇⁻ +end + +function reinit!(::FirstOrderStencil{2,T},φ_new,φ,vel,Δt,Δx,isperiodic,caches) where T + D⁺ʸ, D⁺ˣ, D⁻ʸ, D⁻ˣ, ∇⁺, ∇⁻ = caches + Δx, Δy = Δx + xperiodic,yperiodic = isperiodic + # Prepare shifted lsf + circshift!(D⁺ʸ,φ,(0,-1)); circshift!(D⁻ʸ,φ,(0,1)) + circshift!(D⁺ˣ,φ,(-1,0)); circshift!(D⁻ˣ,φ,(1,0)) + # Sign approximation + ∇⁺ .= @. (D⁺ʸ - D⁻ʸ)/(2Δy); ~yperiodic ? ∇⁺[:,[1,end]] .= zero(T) : 0; + ∇⁻ .= @. (D⁺ˣ - D⁻ˣ)/(2Δx); ~xperiodic ? ∇⁻[[1,end],:] .= zero(T) : 0; + ϵₛ = minimum((Δx,Δy)) + vel .= @. φ/sqrt(φ^2 + ϵₛ^2*(∇⁺^2+∇⁻^2)) + # Forward (+) & Backward (-) + D⁺ʸ .= @. (D⁺ʸ - φ)/Δy; ~yperiodic ? D⁺ʸ[:,end] .= zero(T) : 0; + D⁺ˣ .= @. (D⁺ˣ - φ)/Δx; ~xperiodic ? D⁺ˣ[end,:] .= zero(T) : 0; + D⁻ʸ .= @. (φ - D⁻ʸ)/Δy; ~yperiodic ? D⁻ʸ[:,1] .= zero(T) : 0; + D⁻ˣ .= @. (φ - D⁻ˣ)/Δx; ~xperiodic ? D⁻ˣ[1,:] .= zero(T) : 0; + # Operators + ∇⁺ .= @. sqrt(max(D⁻ʸ,0)^2 + min(D⁺ʸ,0)^2 + max(D⁻ˣ,0)^2 + min(D⁺ˣ,0)^2); + ∇⁻ .= @. sqrt(max(D⁺ʸ,0)^2 + min(D⁻ʸ,0)^2 + max(D⁺ˣ,0)^2 + min(D⁻ˣ,0)^2); + # Update + φ_new .= @. φ - Δt*(max(vel,0)*∇⁺ + min(vel,0)*∇⁻ - vel) + return φ_new +end + +function evolve!(::FirstOrderStencil{2,T},φ,vel,Δt,Δx,isperiodic,caches) where T + D⁺ʸ, D⁺ˣ, D⁻ʸ, D⁻ˣ, ∇⁺, ∇⁻ = caches + Δx, Δy = Δx + xperiodic,yperiodic = isperiodic + # Prepare shifted lsf + circshift!(D⁺ʸ,φ,(0,-1)); circshift!(D⁻ʸ,φ,(0,1)) + circshift!(D⁺ˣ,φ,(-1,0)); circshift!(D⁻ˣ,φ,(1,0)) + # Forward (+) & Backward (-) + D⁺ʸ .= @. (D⁺ʸ - φ)/Δy; ~yperiodic ? D⁺ʸ[:,end] .= zero(T) : 0; + D⁺ˣ .= @. (D⁺ˣ - φ)/Δx; ~xperiodic ? D⁺ˣ[end,:] .= zero(T) : 0; + D⁻ʸ .= @. (φ - D⁻ʸ)/Δy; ~yperiodic ? D⁻ʸ[:,1] .= zero(T) : 0; + D⁻ˣ .= @. (φ - D⁻ˣ)/Δx; ~xperiodic ? D⁻ˣ[1,:] .= zero(T) : 0; + # Operators + ∇⁺ .= @. sqrt(max(D⁻ʸ,0)^2 + min(D⁺ʸ,0)^2 + max(D⁻ˣ,0)^2 + min(D⁺ˣ,0)^2) + ∇⁻ .= @. sqrt(max(D⁺ʸ,0)^2 + min(D⁻ʸ,0)^2 + max(D⁺ˣ,0)^2 + min(D⁻ˣ,0)^2) + # Update + φ .= @. φ - Δt*(max(vel,0)*∇⁺ + min(vel,0)*∇⁻) + return φ +end + +function allocate_caches(::FirstOrderStencil{3},φ,vel) + D⁺ᶻ = similar(φ); D⁺ʸ = similar(φ); D⁺ˣ = similar(φ) + D⁻ᶻ = similar(φ); D⁻ʸ = similar(φ); D⁻ˣ = similar(φ) + ∇⁺ = similar(φ); ∇⁻ = similar(φ) + return D⁺ᶻ, D⁺ʸ, D⁺ˣ, D⁻ᶻ, D⁻ʸ, D⁻ˣ, ∇⁺, ∇⁻ +end + +function reinit!(::FirstOrderStencil{3,T},φ_new,φ,vel,Δt,Δx,isperiodic,caches) where T + D⁺ᶻ, D⁺ʸ, D⁺ˣ, D⁻ᶻ, D⁻ʸ, D⁻ˣ, ∇⁺, ∇⁻=caches + Δx, Δy, Δz = Δx + xperiodic,yperiodic,zperiodic = isperiodic + # Prepare shifted lsf + circshift!(D⁺ʸ,φ,(0,-1,0)); circshift!(D⁻ʸ,φ,(0,1,0)) + circshift!(D⁺ˣ,φ,(-1,0,0)); circshift!(D⁻ˣ,φ,(1,0,0)) + circshift!(D⁺ᶻ,φ,(0,0,-1)); circshift!(D⁻ᶻ,φ,(0,0,1)) + # Sign approximation + ∇⁺ .= @. (D⁺ʸ - D⁻ʸ)/(2Δy); ~yperiodic ? ∇⁺[:,[1,end],:] .= zero(T) : 0; + ∇⁻ .= @. (D⁺ˣ - D⁻ˣ)/(2Δx); ~xperiodic ? ∇⁻[[1,end],:,:] .= zero(T) : 0; + ∇⁺ .= @. ∇⁺^2+∇⁻^2 # |∇φ|² (partially computed) + ∇⁻ .= @. (D⁺ᶻ - D⁻ᶻ)/(2Δz); ~xperiodic ? ∇⁻[:,:,[1,end]] .= zero(T) : 0; + ϵₛ = minimum((Δx,Δy,Δz)) + vel .= @. φ/sqrt(φ^2 + ϵₛ^2*(∇⁺+∇⁻^2)) + # Forward (+) & Backward (-) + D⁺ʸ .= (D⁺ʸ - φ)/Δy; ~yperiodic ? D⁺ʸ[:,end,:] .= zero(T) : 0; + D⁺ˣ .= (D⁺ˣ - φ)/Δx; ~xperiodic ? D⁺ˣ[end,:,:] .= zero(T) : 0; + D⁺ᶻ .= (D⁺ᶻ - φ)/Δz; ~zperiodic ? D⁺ᶻ[:,:,end] .= zero(T) : 0; + D⁻ʸ .= (φ - D⁻ʸ)/Δy; ~yperiodic ? D⁻ʸ[:,1,:] .= zero(T) : 0; + D⁻ˣ .= (φ - D⁻ˣ)/Δx; ~xperiodic ? D⁻ˣ[1,:,:] .= zero(T) : 0; + D⁻ᶻ .= (φ - D⁻ᶻ)/Δz; ~zperiodic ? D⁻ᶻ[:,:,1] .= zero(T) : 0; + # Operators + ∇⁺ .= @. sqrt(max(D⁻ʸ,0)^2 + min(D⁺ʸ,0)^2 + max(D⁻ˣ,0)^2 + min(D⁺ˣ,0)^2 + max(D⁻ᶻ,0)^2 + min(D⁺ᶻ,0)^2) + ∇⁻ .= @. sqrt(max(D⁺ʸ,0)^2 + min(D⁻ʸ,0)^2 + max(D⁺ˣ,0)^2 + min(D⁻ˣ,0)^2 + max(D⁺ᶻ,0)^2 + min(D⁻ᶻ,0)^2) + # Update + φ_new .= @. φ - Δt*(max(vel,0)*∇⁺ + min(vel,0)*∇⁻ - vel) + return φ_new +end + +function evolve!(::FirstOrderStencil{3,T},φ,vel,Δt,Δx,isperiodic,caches) where T + D⁺ᶻ, D⁺ʸ, D⁺ˣ, D⁻ᶻ, D⁻ʸ, D⁻ˣ, ∇⁺, ∇⁻=caches + Δx, Δy, Δz = Δx + xperiodic,yperiodic,zperiodic = isperiodic + # Prepare shifted lsf + circshift!(D⁺ʸ,φ,(0,-1,0)); circshift!(D⁻ʸ,φ,(0,1,0)) + circshift!(D⁺ˣ,φ,(-1,0,0)); circshift!(D⁻ˣ,φ,(1,0,0)) + circshift!(D⁺ᶻ,φ,(0,0,-1)); circshift!(D⁻ᶻ,φ,(0,0,1)) + # Forward (+) & Backward (-) + D⁺ʸ .= (D⁺ʸ - φ)/Δy; ~yperiodic ? D⁺ʸ[:,end,:] .= zero(T) : 0; + D⁺ˣ .= (D⁺ˣ - φ)/Δx; ~xperiodic ? D⁺ˣ[end,:,:] .= zero(T) : 0; + D⁺ᶻ .= (D⁺ᶻ - φ)/Δz; ~zperiodic ? D⁺ᶻ[:,:,end] .= zero(T) : 0; + D⁻ʸ .= (φ - D⁻ʸ)/Δy; ~yperiodic ? D⁻ʸ[:,1,:] .= zero(T) : 0; + D⁻ˣ .= (φ - D⁻ˣ)/Δx; ~xperiodic ? D⁻ˣ[1,:,:] .= zero(T) : 0; + D⁻ᶻ .= (φ - D⁻ᶻ)/Δz; ~zperiodic ? D⁻ᶻ[:,:,1] .= zero(T) : 0; + # Operators + ∇⁺ .= @. sqrt(max(D⁻ʸ,0)^2 + min(D⁺ʸ,0)^2 + max(D⁻ˣ,0)^2 + min(D⁺ˣ,0)^2 + max(D⁻ᶻ,0)^2 + min(D⁺ᶻ,0)^2) + ∇⁻ .= @. sqrt(max(D⁺ʸ,0)^2 + min(D⁻ʸ,0)^2 + max(D⁺ˣ,0)^2 + min(D⁻ˣ,0)^2 + max(D⁺ᶻ,0)^2 + min(D⁻ᶻ,0)^2) + # Update + φ .= @. φ - Δt*(max(vel,0)*∇⁺ + min(vel,0)*∇⁻) + return φ +end \ No newline at end of file diff --git a/src/LevelSetTopOpt.jl b/src/LevelSetTopOpt.jl index 7b6ebc4d..5de30155 100644 --- a/src/LevelSetTopOpt.jl +++ b/src/LevelSetTopOpt.jl @@ -52,10 +52,10 @@ export initial_lsf export get_el_Δ export isotropic_elast_tensor -include("Advection.jl") -export AdvectionStencil +include("LevelSetEvolution/LevelSetEvolution.jl") +export HamiltonJacobiEvolution export FirstOrderStencil -export advect! +export evolve! export reinit! include("Solvers.jl") diff --git a/src/Optimisers/AugmentedLagrangian.jl b/src/Optimisers/AugmentedLagrangian.jl index 20bcfd19..6c2d48b5 100644 --- a/src/Optimisers/AugmentedLagrangian.jl +++ b/src/Optimisers/AugmentedLagrangian.jl @@ -1,5 +1,5 @@ """ - struct AugmentedLagrangian{N,O} <: Optimiser + struct AugmentedLagrangian <: Optimiser An augmented Lagrangian method based on Nocedal and Wright, 2006 ([link](https://doi.org/10.1007/978-0-387-40065-5)). Note that @@ -8,9 +8,8 @@ are defined in `problem::PDEConstrainedFunctionals`. # Parameters -- `problem::PDEConstrainedFunctionals{N}`: The objective and constraint setup. -- `stencil::AdvectionStencil{O}`: The finite difference stencil for - solving the evolution and reinitisation equations. +- `problem::PDEConstrainedFunctionals`: The objective and constraint setup. +- `ls_evolver::LevelSetEvolution`: Solver for the evolution and reinitisation equations. - `vel_ext::VelocityExtension`: The velocity-extension method for extending shape sensitivities onto the computational domain. - `history::OptimiserHistory{Float64}`: Historical information for optimisation problem. @@ -22,9 +21,9 @@ The `has_oscillations` function has been added to avoid oscillations in the iteration history. By default this uses a mean zero crossing algorithm as implemented in ChaosTools. Oscillations checking can be disabled by taking `has_oscillations = (args...) -> false`. """ -struct AugmentedLagrangian{N,O} <: Optimiser - problem :: PDEConstrainedFunctionals{N} - stencil :: AdvectionStencil{O} +struct AugmentedLagrangian <: Optimiser + problem :: PDEConstrainedFunctionals + ls_evolver :: LevelSetEvolution vel_ext :: VelocityExtension history :: OptimiserHistory{Float64} converged :: Function @@ -34,9 +33,9 @@ struct AugmentedLagrangian{N,O} <: Optimiser @doc """ AugmentedLagrangian( - problem :: PDEConstrainedFunctionals{N}, - stencil :: AdvectionStencil{O}, - vel_ext :: VelocityExtension, + problem :: PDEConstrainedFunctionals{N}, + ls_evolver :: LevelSetEvolution, + vel_ext :: VelocityExtension, φ0; Λ_max = 10^10, ζ = 1.1, update_mod = 5, γ = 0.1, γ_reinit = 0.5, os_γ_mult = 0.75, maxiter = 1000, verbose=false, constraint_names = map(i -> Symbol("C_\$i"),1:N), @@ -48,9 +47,8 @@ struct AugmentedLagrangian{N,O} <: Optimiser # Required - - `problem::PDEConstrainedFunctionals{N}`: The objective and constraint setup. - - `stencil::AdvectionStencil{O}`: The finite difference stencil for - solving the evolution and reinitisation equations. + - `problem::PDEConstrainedFunctionals`: The objective and constraint setup. + - `ls_evolver::LevelSetEvolution`: Solver for the evolution and reinitisation equations. - `vel_ext::VelocityExtension`: The velocity-extension method for extending shape sensitivities onto the computational domain. - `φ0`: An initial level-set function defined as a FEFunction or GridapDistributed equivilent. @@ -75,16 +73,16 @@ struct AugmentedLagrangian{N,O} <: Optimiser - `debug = false`: Debug flag. """ function AugmentedLagrangian( - problem :: PDEConstrainedFunctionals{N}, - stencil :: AdvectionStencil{O}, - vel_ext :: VelocityExtension, + problem :: PDEConstrainedFunctionals{N}, + ls_evolver :: LevelSetEvolution, + vel_ext :: VelocityExtension, φ0; Λ_max = 10^10, ζ = 1.1, update_mod = 5, reinit_mod = 1, γ = 0.1, γ_reinit = 0.5, os_γ_mult = 0.75, maxiter = 1000, verbose=false, constraint_names = map(i -> Symbol("C_$i"),1:N), converged::Function = default_al_converged, debug = false, has_oscillations::Function = default_has_oscillations, initial_parameters::Function = default_al_init_params - ) where {N,O} + ) where N constraint_names = map(Symbol,constraint_names) λ_names = map(i -> Symbol("λ$i"),1:N) @@ -94,7 +92,7 @@ struct AugmentedLagrangian{N,O} <: Optimiser history = OptimiserHistory(Float64,al_keys,al_bundles,maxiter,verbose) params = (;Λ_max,ζ,update_mod,reinit_mod,γ,γ_reinit,os_γ_mult,debug,initial_parameters) - new{N,O}(problem,stencil,vel_ext,history,converged,has_oscillations,params,φ0) + new(problem,ls_evolver,vel_ext,history,converged,has_oscillations,params,φ0) end end @@ -125,7 +123,7 @@ end function default_al_converged( m::AugmentedLagrangian; - L_tol = 0.01*maximum(m.stencil.params.Δ)/(length(m.stencil.params.Δ)-1), + L_tol = 0.01*maximum(get_dof_Δ(m.ls_evolver))/(length(get_dof_Δ(m.ls_evolver))-1), C_tol = 0.01 ) h = m.history @@ -145,7 +143,7 @@ function Base.iterate(m::AugmentedLagrangian) φh, history, params = m.φ0, m.history, m.params ## Reinitialise as SDF - reinit!(m.stencil,φh,params.γ_reinit) + reinit!(m.ls_evolver,φh,params.γ_reinit) ## Compute FE problem and shape derivatives J, C, dJ, dC = evaluate!(m.problem,φh) @@ -176,8 +174,7 @@ end function Base.iterate(m::AugmentedLagrangian,state) it, L, J, C, dL, dJ, dC, uh, φh, vel, λ, Λ, γ, os_it = state params, history = m.params, m.history - update_mod, reinit_mod, ζ, Λ_max, γ_reinit, os_γ_mult = params.update_mod, - params.reinit_mod, params.ζ, params.Λ_max, params.γ_reinit, params.os_γ_mult + Λ_max,ζ,update_mod,reinit_mod,_,γ_reinit,os_γ_mult,_,_ = params ## Periodicially call GC iszero(it % 50) && GC.gc(); @@ -197,8 +194,8 @@ function Base.iterate(m::AugmentedLagrangian,state) U_reg = get_deriv_space(m.problem.state_map) V_φ = get_aux_space(m.problem.state_map) interpolate!(FEFunction(U_reg,dL),vel,V_φ) - advect!(m.stencil,φh,vel,γ) - iszero(it % reinit_mod) && reinit!(m.stencil,φh,γ_reinit) + evolve!(m.ls_evolver,φh,vel,γ) + iszero(it % reinit_mod) && reinit!(m.ls_evolver,φh,γ_reinit) ## Calculate objective, constraints, and shape derivatives J, C, dJ, dC = evaluate!(m.problem,φh) diff --git a/src/Optimisers/HilbertianProjection.jl b/src/Optimisers/HilbertianProjection.jl index 9cb75952..b4f18628 100644 --- a/src/Optimisers/HilbertianProjection.jl +++ b/src/Optimisers/HilbertianProjection.jl @@ -124,7 +124,7 @@ function compute_α(C,dC_orthog,dC,normsq,K,P,λ,α_min,α_max) end """ - struct HilbertianProjection{T,N} <: Optimiser + struct HilbertianProjection <: Optimiser A Hilbertian projection method as described by Wegert et al., 2023 ([link](https://doi.org/10.1007/s00158-023-03663-0)). @@ -132,8 +132,7 @@ A Hilbertian projection method as described by Wegert et al., 2023 # Parameters - `problem::PDEConstrainedFunctionals{N}`: The objective and constraint setup. -- `stencil::AdvectionStencil{O}`: The finite difference stencil for - solving the evolution and reinitisation equations. +- `ls_evolver::LevelSetEvolution`: Solver for the evolution and reinitisation equations. - `vel_ext::VelocityExtension`: The velocity-extension method for extending shape sensitivities onto the computational domain. - `projector::HilbertianProjectionMap`: Sensitivity information projector @@ -142,9 +141,9 @@ A Hilbertian projection method as described by Wegert et al., 2023 - `has_oscillations::Function`: A function to check for oscillations. - `params::NamedTuple`: Optimisation parameters. """ -struct HilbertianProjection{T,N} <: Optimiser - problem :: PDEConstrainedFunctionals{N} - stencil :: AdvectionStencil +struct HilbertianProjection <: Optimiser + problem :: PDEConstrainedFunctionals + ls_evolver :: LevelSetEvolution vel_ext :: VelocityExtension projector :: HilbertianProjectionMap history :: OptimiserHistory{Float64} @@ -156,7 +155,7 @@ struct HilbertianProjection{T,N} <: Optimiser @doc """ HilbertianProjection( problem :: PDEConstrainedFunctionals{N}, - stencil :: AdvectionStencil, + ls_evolver :: LevelSetEvolution, vel_ext :: VelocityExtension, φ0; orthog = HPModifiedGramSchmidt(), @@ -174,8 +173,7 @@ struct HilbertianProjection{T,N} <: Optimiser # Required - `problem::PDEConstrainedFunctionals{N}`: The objective and constraint setup. - - `stencil::AdvectionStencil{O}`: The finite difference stencil for - solving the evolution and reinitisation equations. + - `ls_evolver::LevelSetEvolution`: Solver for the evolution and reinitisation equations. - `vel_ext::VelocityExtension`: The velocity-extension method for extending shape sensitivities onto the computational domain. - `φ0`: An initial level-set function defined as a FEFunction or GridapDistributed equivilent. @@ -221,9 +219,9 @@ struct HilbertianProjection{T,N} <: Optimiser disabling the line search via `ls_enabled = false` will enable oscillation detection. """ function HilbertianProjection( - problem :: PDEConstrainedFunctionals{N}, - stencil :: AdvectionStencil, - vel_ext :: VelocityExtension, + problem :: PDEConstrainedFunctionals{N}, + ls_evolver :: LevelSetEvolution, + vel_ext :: VelocityExtension, φ0; orthog = HPModifiedGramSchmidt(), λ=0.5, α_min=0.1, α_max=1.0, γ=0.1, γ_reinit=0.5, reinit_mod = 1, @@ -234,7 +232,9 @@ struct HilbertianProjection{T,N} <: Optimiser converged::Function = default_hp_converged, debug = false, has_oscillations::Function = (ls_enabled ? (args...)->false : default_has_oscillations), os_γ_mult = 0.5 - ) where {N} + ) where N + + @assert α_min <= α_max "We require α_min <= α_max" constraint_names = map(Symbol,constraint_names) al_keys = [:J,constraint_names...,:γ] @@ -244,8 +244,7 @@ struct HilbertianProjection{T,N} <: Optimiser projector = HilbertianProjectionMap(N,orthog,vel_ext;λ,α_min,α_max,debug) params = (;debug,γ,γ_reinit,reinit_mod,ls_enabled,ls_max_iters,ls_δ_inc,ls_δ_dec,ls_ξ, ls_ξ_reduce_coef,ls_ξ_reduce_abs_tol,ls_γ_min,ls_γ_max,os_γ_mult) - T = typeof(orthog) - new{T,N}(problem,stencil,vel_ext,projector,history,converged, + new(problem,ls_evolver,vel_ext,projector,history,converged, has_oscillations,params,φ0) end end @@ -258,7 +257,7 @@ end function default_hp_converged( m::HilbertianProjection; - J_tol = 0.2*maximum(m.stencil.params.Δ), + J_tol = 0.2*maximum(get_dof_Δ(m.ls_evolver)), C_tol = 0.001 ) h = m.history @@ -296,7 +295,7 @@ function Base.iterate(m::HilbertianProjection) φh = m.φ0 ## Reinitialise as SDF - reinit!(m.stencil,φh,params.γ_reinit) + reinit!(m.ls_evolver,φh,params.γ_reinit) ## Compute FE problem and shape derivatives J, C, dJ, dC = Gridap.evaluate!(m.problem,φh) @@ -351,8 +350,8 @@ function Base.iterate(m::HilbertianProjection,state) φ = get_free_dof_values(φh); copy!(φ_tmp,φ) while !done && (ls_it <= ls_max_iters) # Advect & Reinitialise - advect!(m.stencil,φ,vel,γ) - iszero(it % reinit_mod) && reinit!(m.stencil,φ,params.γ_reinit) + evolve!(m.ls_evolver,φ,vel,γ) + iszero(it % reinit_mod) && reinit!(m.ls_evolver,φ,params.γ_reinit) ~ls_enabled && break diff --git a/test/seq/ThermalComplianceALMTests.jl b/test/seq/ThermalComplianceALMTests.jl index 1d227f11..a5891838 100644 --- a/test/seq/ThermalComplianceALMTests.jl +++ b/test/seq/ThermalComplianceALMTests.jl @@ -71,7 +71,7 @@ function main() dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ ## Finite difference solver and level set function - stencil = AdvectionStencil(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) + ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) ## Setup solver and FE operators state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N) @@ -83,7 +83,7 @@ function main() vel_ext = VelocityExtension(a_hilb,U_reg,V_reg) ## Optimiser - optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh; + optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh; γ,γ_reinit,verbose=true,constraint_names=[:Vol]) # Do a few iterations