From 4b94785e030ac2bc69ac58700b517ed5ea938ec6 Mon Sep 17 00:00:00 2001 From: "Z J Wegert (Workstation)" Date: Wed, 10 Apr 2024 10:43:19 +1000 Subject: [PATCH 01/12] Create LevelSetEvolution abstract type --- .../HamiltonJacobiEvolution.jl | 351 ++++++++++++++++++ src/LevelSetEvolution/LevelSetEvolution.jl | 39 ++ src/LevelSetEvolution/SpatialStencil.jl | 175 +++++++++ 3 files changed, 565 insertions(+) create mode 100644 src/LevelSetEvolution/HamiltonJacobiEvolution.jl create mode 100644 src/LevelSetEvolution/LevelSetEvolution.jl create mode 100644 src/LevelSetEvolution/SpatialStencil.jl diff --git a/src/LevelSetEvolution/HamiltonJacobiEvolution.jl b/src/LevelSetEvolution/HamiltonJacobiEvolution.jl new file mode 100644 index 00000000..8e86f07a --- /dev/null +++ b/src/LevelSetEvolution/HamiltonJacobiEvolution.jl @@ -0,0 +1,351 @@ +""" + struct HamiltonJacobiEvolution{O} + +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. + +Based on the scheme by Osher and Fedkiw ([link](https://doi.org/10.1007/b98879)). + +# Parameters + +- `stencil::SpatialStencil`: Spatial finite difference stencil for a single step HJ + equation and reinitialisation equation. +- `model`: A `CartesianDiscreteModel`. +- `space`: FE space for level-set function +- `perm`: A permutation vector +- `params`: Tuple of additional params +- `cache`: SpatialStencil cache +""" +struct HamiltonJacobiEvolution{O} <: LevelSetEvolution + stencil :: SpatialStencil + model + space + perm + params + cache +end + +""" + HamiltonJacobiEvolution(stencil::SpatialStencil,model,space,tol,max_steps,max_steps_reinit) + +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 HamiltonJacobiEvolution( + stencil::SpatialStencil, + model, + space, + tol=1.e-3, + max_steps=100, + max_steps_reinit=2000 +) + # Parameters + 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) + + # Caches + φ, vel = zero_free_values(space), zero_free_values(space) + cache = allocate_caches(stencil,φ,vel,perm,order,ndof) + + return HamiltonJacobiEvolution{order}(stencil,model,space,perm,params,cache) +end + +""" + get_dof_Δ(::HamiltonJacobiEvolution) + +Return the distance betweem degree of freedom +""" +function get_dof_Δ(m::HamiltonJacobiEvolution) + return m.params.Δ +end + +""" + compute_Δt(::HamiltonJacobiEvolution,φ,vel) + +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 + +""" + evolve!(s::HamiltonJacobiEvolution{O},φ,vel,γ) where O + +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,`` + +with ``\\phi(0,\\boldsymbol{x})=\\phi_0(\\boldsymbol{x})`` and ``\\boldsymbol{x}\\in D,~t\\in(0,T)``. + +# Arguments + +- `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 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 + + _φ = (O >= 2) ? permute!(perm_caches[1],φ,s.perm) : φ + _vel = (O >= 2) ? permute!(perm_caches[2],vel,s.perm) : vel + + ## CFL Condition (requires γ≤1.0) + Δ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) + evolve!(s.stencil,φ_mat,vel_mat,Δt,Δ,isperiodic,stencil_cache) + end + # Update ghost nodes + consistent!(_φ) |> fetch + end + φ = (O >= 2) ? permute_inv!(φ,_φ,s.perm) : _φ + return φ +end + +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 + + _φ = (O >= 2) ? permute!(perm_caches[1],φ,s.perm) : φ + _vel = (O >= 2) ? permute!(perm_caches[2],vel,s.perm) : vel + + ## CFL Condition (requires γ≤1.0) + Δt = compute_Δt(s,Δ,γ,φ,vel) + for _ in 1:max_steps + φ_mat = reshape(_φ,ndof) + vel_mat = reshape(_vel,ndof) + evolve!(s.stencil,φ_mat,vel_mat,Δt,Δ,isperiodic,stencil_cache) + end + φ = (O >= 2) ? permute_inv!(φ,_φ,s.perm) : _φ + return φ +end + +function evolve!(s::HamiltonJacobiEvolution,φh,args...) + evolve!(s,get_free_dof_values(φh),args...) +end + +""" + reinit!(s::HamiltonJacobiEvolution{O},φ,γ) where O + +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,`` + +with ``\\phi(0,\\boldsymbol{x})=\\phi_0(\\boldsymbol{x})`` and ``\\boldsymbol{x}\\in D,~t\\in(0,T)``. + +# Arguments + +- `s::HamiltonJacobiEvolution{O}`: Method +- `φ`: level set function as a vector of degrees of freedom +- `γ`: coeffient on the time step size. +""" +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 infnorm(S) = 1.0 + Δt = compute_Δt(s,Δ,γ,_φ,1.0) + + # Apply operations across partitions + step = 1; err = maximum(abs,φ); fill!(φ_tmp,0.0) + while (err > tol) && (step <= max_steps) + # Step of 1st order upwind reinitialisation equation + map(local_views(φ_tmp),local_views(_φ),local_views(vel_tmp),stencil_cache,ndof) do φ_tmp,_φ,vel_tmp,stencil_cache,S + φ_tmp_mat = reshape(φ_tmp,S) + φ_mat = reshape(_φ,S) + vel_tmp_mat = reshape(vel_tmp,S) + reinit!(s.stencil,φ_tmp_mat,φ_mat,vel_tmp_mat,Δt,Δ,isperiodic,stencil_cache) + end + + # Compute error + _φ .-= φ_tmp # φ - φ_tmp + err = maximum(abs,_φ) + step += 1 + + # Update φ + copy!(_φ,φ_tmp) + consistent!(_φ) |> fetch # We exchange ghosts here! + end + φ = (O >= 2) ? permute_inv!(φ,_φ,s.perm) : _φ + return φ +end + +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 + + _φ = (O >= 2) ? permute!(perm_caches[1],φ,s.perm) : φ + + ## CFL Condition (requires γ≤0.5) + Δt = compute_Δt(s,Δ,γ,_φ,1.0) + + # Apply operations across partitions + step = 1; err = maximum(abs,φ); fill!(φ_tmp,0.0) + while (err > tol) && (step <= max_steps) + # Step of 1st order upwind reinitialisation equation + φ_tmp_mat = reshape(φ_tmp,ndof) + φ_mat = reshape(_φ,ndof) + vel_tmp_mat = reshape(vel_tmp,ndof) + reinit!(s.stencil,φ_tmp_mat,φ_mat,vel_tmp_mat,Δt,Δ,isperiodic,stencil_cache) + + # Compute error + _φ .-= φ_tmp + err = maximum(abs,_φ) + step += 1 + + # Update φ + copy!(_φ,φ_tmp) + end + φ = (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::SpatialStencil,φ::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::SpatialStencil,φ::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..d6a743f4 --- /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("SpatialStencil.jl") +include("HamiltonJacobiEvolution.jl") \ No newline at end of file diff --git a/src/LevelSetEvolution/SpatialStencil.jl b/src/LevelSetEvolution/SpatialStencil.jl new file mode 100644 index 00000000..12dd6469 --- /dev/null +++ b/src/LevelSetEvolution/SpatialStencil.jl @@ -0,0 +1,175 @@ +""" + abstract type SpatialStencil + +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 SpatialStencil end + +""" + allocate_caches(::SpatialStencil,φ,vel) + +Allocate caches for a given `SpatialStencil`. +""" +function allocate_caches(::SpatialStencil,φ,vel) + nothing # By default, no caches are required. +end + +""" + check_order(::SpatialStencil,order) + +Throw error if insufficient reference element order +to implement stencil in parallel. +""" +function check_order(::SpatialStencil,order) + @abstractmethod +end + +""" + reinit!(::SpatialStencil,φ_new,φ,vel,Δt,Δx,isperiodic,caches) -> φ + +Single finite difference step of the reinitialisation equation for a given `SpatialStencil`. +""" +function reinit!(::SpatialStencil,φ_new,φ,vel,Δt,Δx,isperiodic,caches) + @abstractmethod +end + +""" + evolve!(::SpatialStencil,φ,vel,Δt,Δx,isperiodic,caches) -> φ + +Single finite difference step of the Hamilton-Jacobi evoluation equation for a given +`SpatialStencil`. +""" +function evolve!(::SpatialStencil,φ,vel,Δt,Δx,isperiodic,caches) + @abstractmethod +end + +""" + struct FirstOrderStencil{D,T} <: SpatialStencil end + +A first order Godunov upwind difference scheme based on Osher and Fedkiw +([link](https://doi.org/10.1007/b98879)). +""" +struct FirstOrderStencil{D,T} <: SpatialStencil + 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 From 3d0539e083f46e4f4f56263cc896688d0669c962 Mon Sep 17 00:00:00 2001 From: "Z J Wegert (Workstation)" Date: Wed, 10 Apr 2024 10:43:45 +1000 Subject: [PATCH 02/12] Update includes and exports --- src/LevelSetTopOpt.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/LevelSetTopOpt.jl b/src/LevelSetTopOpt.jl index 7b6ebc4d..cc30a752 100644 --- a/src/LevelSetTopOpt.jl +++ b/src/LevelSetTopOpt.jl @@ -52,8 +52,8 @@ 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 reinit! From 8850dfa052c64bdf4178f469203b44873c7b8444 Mon Sep 17 00:00:00 2001 From: "Z J Wegert (Workstation)" Date: Wed, 10 Apr 2024 10:46:14 +1000 Subject: [PATCH 03/12] Update Benchmark.jl --- src/Benchmarks.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/Benchmarks.jl b/src/Benchmarks.jl index 87cc4a6e..afffef2d 100644 --- a/src/Benchmarks.jl +++ b/src/Benchmarks.jl @@ -116,13 +116,13 @@ 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. """ -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,γ) end @@ -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 From 8b60442857aa7e78aeddd64f47ba0ff2d14b4609 Mon Sep 17 00:00:00 2001 From: "Z J Wegert (Workstation)" Date: Wed, 10 Apr 2024 11:08:07 +1000 Subject: [PATCH 04/12] Rename advect! --- src/Benchmarks.jl | 4 ++-- src/LevelSetTopOpt.jl | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/Benchmarks.jl b/src/Benchmarks.jl index afffef2d..0c052a44 100644 --- a/src/Benchmarks.jl +++ b/src/Benchmarks.jl @@ -120,11 +120,11 @@ end 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::LevelSetEvolution, φ0, v0, γ, ranks; nreps = 10) function f(stencil,φ,v,γ) - advect!(stencil,φ,v,γ) + evolve!(stencil,φ,v,γ) end function reset!(stencil,φ,v,γ) copy!(φ,φ0) diff --git a/src/LevelSetTopOpt.jl b/src/LevelSetTopOpt.jl index cc30a752..5de30155 100644 --- a/src/LevelSetTopOpt.jl +++ b/src/LevelSetTopOpt.jl @@ -55,7 +55,7 @@ export isotropic_elast_tensor include("LevelSetEvolution/LevelSetEvolution.jl") export HamiltonJacobiEvolution export FirstOrderStencil -export advect! +export evolve! export reinit! include("Solvers.jl") From 489f886c96d485e9070e6221bbc7ecab97a7751a Mon Sep 17 00:00:00 2001 From: "Z J Wegert (Workstation)" Date: Wed, 10 Apr 2024 11:22:10 +1000 Subject: [PATCH 05/12] Update src/Optimisers/.. --- src/Optimisers/AugmentedLagrangian.jl | 45 ++++++++++++-------------- src/Optimisers/HilbertianProjection.jl | 37 +++++++++++---------- 2 files changed, 39 insertions(+), 43 deletions(-) diff --git a/src/Optimisers/AugmentedLagrangian.jl b/src/Optimisers/AugmentedLagrangian.jl index 3a44ea75..1c879508 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 = 5.0, ζ = 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.Δ), + L_tol = 0.01*maximum(get_dof_Δ(m.ls_evolver)), 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 From bea8cf2be9c0e989b4d711e1f674668912f45982 Mon Sep 17 00:00:00 2001 From: "Z J Wegert (Workstation)" Date: Wed, 10 Apr 2024 11:25:22 +1000 Subject: [PATCH 06/12] Update scripts --- scripts/Benchmarks/benchmark.jl | 8 ++++---- scripts/MPI/3d_elastic_compliance_ALM.jl | 2 +- scripts/MPI/3d_hyperelastic_compliance_ALM.jl | 2 +- scripts/MPI/3d_hyperelastic_compliance_neohook_ALM.jl | 2 +- scripts/MPI/3d_inverse_homenisation_ALM.jl | 2 +- scripts/MPI/3d_inverter_ALM.jl | 2 +- scripts/MPI/3d_inverter_HPM.jl | 2 +- scripts/MPI/3d_nonlinear_thermal_compliance_ALM.jl | 2 +- scripts/MPI/3d_thermal_compliance_ALM.jl | 2 +- scripts/MPI/elastic_compliance_ALM.jl | 2 +- scripts/MPI/inverse_homenisation_ALM.jl | 2 +- scripts/MPI/thermal_compliance_ALM.jl | 2 +- scripts/MPI/thermal_compliance_ALM_PETSc.jl | 2 +- scripts/Serial/elastic_compliance_ALM.jl | 2 +- scripts/Serial/elastic_compliance_HPM.jl | 2 +- scripts/Serial/elastic_compliance_LM.jl | 2 +- scripts/Serial/hyperelastic_compliance_ALM.jl | 2 +- scripts/Serial/hyperelastic_compliance_neohook_ALM.jl | 2 +- scripts/Serial/inverse_homenisation_ALM.jl | 2 +- scripts/Serial/inverter_ALM.jl | 2 +- scripts/Serial/inverter_HPM.jl | 2 +- scripts/Serial/nonlinear_thermal_compliance_ALM.jl | 2 +- scripts/Serial/thermal_compliance_ALM.jl | 2 +- scripts/Serial/thermal_compliance_ALM_higherorderlsf.jl | 2 +- scripts/_dev/3d_hyperelastic_compliance_ALM.jl | 2 +- scripts/_dev/Paper_scripts/elast_comp_MPI_3D.jl | 2 +- scripts/_dev/Paper_scripts/hyperelast_comp_MPI_3D.jl | 2 +- scripts/_dev/Paper_scripts/inverse_hom_MPI_3D.jl | 2 +- scripts/_dev/Paper_scripts/inverter_MPI_3D.jl | 2 +- scripts/_dev/Paper_scripts/therm_comp_MPI.jl | 2 +- scripts/_dev/Paper_scripts/therm_comp_MPI_3D.jl | 2 +- scripts/_dev/Paper_scripts/therm_comp_serial.jl | 2 +- .../_dev/benchmark_tests/MPI_benchmark_test_2D_themal.jl | 2 +- .../_dev/benchmark_tests/MPI_benchmark_test_3D_themal.jl | 2 +- scripts/_dev/benchmark_tests/serial_benchmark_test.jl | 2 +- scripts/_dev/bug_issue_46/thermal_MPI.jl | 2 +- scripts/_dev/bug_issue_46/thermal_serial.jl | 2 +- scripts/_dev/full_piezo_script.jl | 2 +- scripts/_dev/inv_hom_block_assem_testing.jl | 2 +- scripts/_dev/inv_hom_block_assem_testing_3D_MPI.jl | 2 +- scripts/_dev/inv_hom_block_assem_testing_full_script.jl | 2 +- scripts/_dev/mem_leak.jl | 2 +- scripts/_dev/nonlinear_adjoint_MWE.jl | 2 +- .../minimum_thermal_compliance_implementation.jl | 2 +- .../minimum_thermal_compliance_implementation_3d.jl | 2 +- .../minimum_thermal_compliance_implementation_3d_petsc.jl | 2 +- ...imum_thermal_compliance_implementation_3d_petsc_mpi.jl | 2 +- ...minimum_thermal_compliance_implementation_nonlinear.jl | 2 +- 48 files changed, 51 insertions(+), 51 deletions(-) diff --git a/scripts/Benchmarks/benchmark.jl b/scripts/Benchmarks/benchmark.jl index 3c0fe455..759bf6e0 100644 --- a/scripts/Benchmarks/benchmark.jl +++ b/scripts/Benchmarks/benchmark.jl @@ -102,7 +102,7 @@ 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) + stencil = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps) # Optimiser return AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit,verbose) end @@ -180,7 +180,7 @@ 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) + stencil = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps) # Optimiser return AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit,verbose) end @@ -256,7 +256,7 @@ 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) + stencil = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps) # Optimiser return AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit,verbose) end @@ -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) + stencil = HamiltonJacobiEvolution(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps) ## Setup solver and FE operators Tm = SparseMatrixCSR{0,PetscScalar,PetscInt} diff --git a/scripts/MPI/3d_elastic_compliance_ALM.jl b/scripts/MPI/3d_elastic_compliance_ALM.jl index 58836523..d44fdff3 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) + stencil = HamiltonJacobiEvolution(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps) ## Setup solver and FE operators Tm = SparseMatrixCSR{0,PetscScalar,PetscInt} diff --git a/scripts/MPI/3d_hyperelastic_compliance_ALM.jl b/scripts/MPI/3d_hyperelastic_compliance_ALM.jl index 6064ae61..c103beab 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) + stencil = HamiltonJacobiEvolution(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps) ## Setup solver and FE operators Tm = SparseMatrixCSR{0,PetscScalar,PetscInt} diff --git a/scripts/MPI/3d_hyperelastic_compliance_neohook_ALM.jl b/scripts/MPI/3d_hyperelastic_compliance_neohook_ALM.jl index 7b02dd0c..1a9fe02e 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) + stencil = HamiltonJacobiEvolution(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps) ## Setup solver and FE operators Tm = SparseMatrixCSR{0,PetscScalar,PetscInt} diff --git a/scripts/MPI/3d_inverse_homenisation_ALM.jl b/scripts/MPI/3d_inverse_homenisation_ALM.jl index ec6f7872..2fed57d8 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) + stencil = HamiltonJacobiEvolution(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps) ## Setup solver and FE operators Tm = SparseMatrixCSR{0,PetscScalar,PetscInt} diff --git a/scripts/MPI/3d_inverter_ALM.jl b/scripts/MPI/3d_inverter_ALM.jl index 0b600003..5d270952 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) + stencil = HamiltonJacobiEvolution(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps) ## Setup solver and FE operators Tm = SparseMatrixCSR{0,PetscScalar,PetscInt} diff --git a/scripts/MPI/3d_inverter_HPM.jl b/scripts/MPI/3d_inverter_HPM.jl index bd6d16e9..f1684cf0 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) + stencil = HamiltonJacobiEvolution(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps) ## Setup solver and FE operators Tm = SparseMatrixCSR{0,PetscScalar,PetscInt} diff --git a/scripts/MPI/3d_nonlinear_thermal_compliance_ALM.jl b/scripts/MPI/3d_nonlinear_thermal_compliance_ALM.jl index f4f110d6..126309d5 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) + stencil = HamiltonJacobiEvolution(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps) ## Setup solver and FE operators Tm = SparseMatrixCSR{0,PetscScalar,PetscInt} diff --git a/scripts/MPI/3d_thermal_compliance_ALM.jl b/scripts/MPI/3d_thermal_compliance_ALM.jl index 0fa01deb..e0cea609 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) + stencil = HamiltonJacobiEvolution(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps) ## Setup solver and FE operators Tm = SparseMatrixCSR{0,PetscScalar,PetscInt} diff --git a/scripts/MPI/elastic_compliance_ALM.jl b/scripts/MPI/elastic_compliance_ALM.jl index fcfb451f..6fa7a17c 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) + stencil = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) ## Setup solver and FE operators Tm = SparseMatrixCSR{0,PetscScalar,PetscInt} diff --git a/scripts/MPI/inverse_homenisation_ALM.jl b/scripts/MPI/inverse_homenisation_ALM.jl index 7a068021..ec4d8886 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) + stencil = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) ## Setup solver and FE operators Tm = SparseMatrixCSR{0,PetscScalar,PetscInt} diff --git a/scripts/MPI/thermal_compliance_ALM.jl b/scripts/MPI/thermal_compliance_ALM.jl index 8bb1e9aa..bce212f1 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) + stencil = 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) diff --git a/scripts/MPI/thermal_compliance_ALM_PETSc.jl b/scripts/MPI/thermal_compliance_ALM_PETSc.jl index 39bd0271..cfe2acd6 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) + stencil = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) ## Setup solver and FE operators Tm = SparseMatrixCSR{0,PetscScalar,PetscInt} diff --git a/scripts/Serial/elastic_compliance_ALM.jl b/scripts/Serial/elastic_compliance_ALM.jl index c313ef68..7cd20e1a 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) + stencil = 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) diff --git a/scripts/Serial/elastic_compliance_HPM.jl b/scripts/Serial/elastic_compliance_HPM.jl index e7fe84c8..50964928 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) + stencil = 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) diff --git a/scripts/Serial/elastic_compliance_LM.jl b/scripts/Serial/elastic_compliance_LM.jl index 004e2119..0517ef1b 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) + stencil = 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) diff --git a/scripts/Serial/hyperelastic_compliance_ALM.jl b/scripts/Serial/hyperelastic_compliance_ALM.jl index fe45ee8e..c455de7b 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) + stencil = 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) diff --git a/scripts/Serial/hyperelastic_compliance_neohook_ALM.jl b/scripts/Serial/hyperelastic_compliance_neohook_ALM.jl index c212f0ef..03ea9f55 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) + stencil = 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) diff --git a/scripts/Serial/inverse_homenisation_ALM.jl b/scripts/Serial/inverse_homenisation_ALM.jl index 7779b851..c079c08c 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) + stencil = 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Ω) diff --git a/scripts/Serial/inverter_ALM.jl b/scripts/Serial/inverter_ALM.jl index d1144e2f..2774a9ea 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) + stencil = 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) diff --git a/scripts/Serial/inverter_HPM.jl b/scripts/Serial/inverter_HPM.jl index ae21a683..f2587eec 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) + stencil = 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) diff --git a/scripts/Serial/nonlinear_thermal_compliance_ALM.jl b/scripts/Serial/nonlinear_thermal_compliance_ALM.jl index edbc0b1b..51597244 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) + stencil = 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) diff --git a/scripts/Serial/thermal_compliance_ALM.jl b/scripts/Serial/thermal_compliance_ALM.jl index 44903c27..6ddf66cc 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) + stencil = 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) diff --git a/scripts/Serial/thermal_compliance_ALM_higherorderlsf.jl b/scripts/Serial/thermal_compliance_ALM_higherorderlsf.jl index 9123a832..0b076760 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) + stencil = 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) diff --git a/scripts/_dev/3d_hyperelastic_compliance_ALM.jl b/scripts/_dev/3d_hyperelastic_compliance_ALM.jl index 2565536b..de79e6fc 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) + # stencil = 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..c903951a 100644 --- a/scripts/_dev/Paper_scripts/elast_comp_MPI_3D.jl +++ b/scripts/_dev/Paper_scripts/elast_comp_MPI_3D.jl @@ -78,7 +78,7 @@ function main(mesh_partition,distribute) ) # Finite difference scheme scheme = FirstOrderStencil(length(el_size),Float64) - stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps) + stencil = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps) # Optimiser optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit, verbose=i_am_main(ranks)) diff --git a/scripts/_dev/Paper_scripts/hyperelast_comp_MPI_3D.jl b/scripts/_dev/Paper_scripts/hyperelast_comp_MPI_3D.jl index 5b4730c5..6f463318 100644 --- a/scripts/_dev/Paper_scripts/hyperelast_comp_MPI_3D.jl +++ b/scripts/_dev/Paper_scripts/hyperelast_comp_MPI_3D.jl @@ -92,7 +92,7 @@ function main(mesh_partition,distribute) ) # Finite difference scheme scheme = FirstOrderStencil(length(el_size),Float64) - stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps) + stencil = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps) # Optimiser optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit, verbose=i_am_main(ranks)) diff --git a/scripts/_dev/Paper_scripts/inverse_hom_MPI_3D.jl b/scripts/_dev/Paper_scripts/inverse_hom_MPI_3D.jl index d15c424c..657b35c4 100644 --- a/scripts/_dev/Paper_scripts/inverse_hom_MPI_3D.jl +++ b/scripts/_dev/Paper_scripts/inverse_hom_MPI_3D.jl @@ -82,7 +82,7 @@ function main(mesh_partition,distribute) ) # Finite difference scheme scheme = FirstOrderStencil(length(el_size),Float64) - stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps) + stencil = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps) # Optimiser optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit, verbose=i_am_main(ranks)) diff --git a/scripts/_dev/Paper_scripts/inverter_MPI_3D.jl b/scripts/_dev/Paper_scripts/inverter_MPI_3D.jl index b9f3269c..75163cb0 100644 --- a/scripts/_dev/Paper_scripts/inverter_MPI_3D.jl +++ b/scripts/_dev/Paper_scripts/inverter_MPI_3D.jl @@ -88,7 +88,7 @@ function main(mesh_partition,distribute) ) # Finite difference scheme scheme = FirstOrderStencil(length(el_size),Float64) - stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps) + stencil = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps) # Optimiser optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit, verbose=i_am_main(ranks)) diff --git a/scripts/_dev/Paper_scripts/therm_comp_MPI.jl b/scripts/_dev/Paper_scripts/therm_comp_MPI.jl index 9b76c455..8320ff76 100644 --- a/scripts/_dev/Paper_scripts/therm_comp_MPI.jl +++ b/scripts/_dev/Paper_scripts/therm_comp_MPI.jl @@ -78,7 +78,7 @@ function main(mesh_partition,distribute) ) # Finite difference scheme scheme = FirstOrderStencil(length(el_size),Float64) - stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps) + stencil = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps) # Optimiser optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit,verbose=i_am_main(ranks)) # Solve diff --git a/scripts/_dev/Paper_scripts/therm_comp_MPI_3D.jl b/scripts/_dev/Paper_scripts/therm_comp_MPI_3D.jl index 99be6ec8..3cb4f8ba 100644 --- a/scripts/_dev/Paper_scripts/therm_comp_MPI_3D.jl +++ b/scripts/_dev/Paper_scripts/therm_comp_MPI_3D.jl @@ -80,7 +80,7 @@ function main(mesh_partition,distribute) ) # Finite difference scheme scheme = FirstOrderStencil(length(el_size),Float64) - stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps) + stencil = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps) # Optimiser optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit, verbose=i_am_main(ranks)) diff --git a/scripts/_dev/Paper_scripts/therm_comp_serial.jl b/scripts/_dev/Paper_scripts/therm_comp_serial.jl index 78cd651c..9e7e2342 100644 --- a/scripts/_dev/Paper_scripts/therm_comp_serial.jl +++ b/scripts/_dev/Paper_scripts/therm_comp_serial.jl @@ -63,7 +63,7 @@ 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) + stencil = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps) # Optimiser optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit, verbose=true) 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..2d09b546 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) + stencil = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) reinit!(stencil,φh,γ_reinit) ## Setup solver and FE operators 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..5842ead6 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) + stencil = HamiltonJacobiEvolution(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps) reinit!(stencil,φh,γ_reinit) ## Setup solver and FE operators diff --git a/scripts/_dev/benchmark_tests/serial_benchmark_test.jl b/scripts/_dev/benchmark_tests/serial_benchmark_test.jl index 8de90168..0f8a784a 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) + stencil = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) reinit!(stencil,φh,γ_reinit) ## Setup solver and FE operators diff --git a/scripts/_dev/bug_issue_46/thermal_MPI.jl b/scripts/_dev/bug_issue_46/thermal_MPI.jl index 2e4596ee..176e7b05 100644 --- a/scripts/_dev/bug_issue_46/thermal_MPI.jl +++ b/scripts/_dev/bug_issue_46/thermal_MPI.jl @@ -68,7 +68,7 @@ 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) + stencil = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps) # Optimiser optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh; γ,γ_reinit,verbose=i_am_main(ranks),constraint_names=[:Vol]) diff --git a/scripts/_dev/bug_issue_46/thermal_serial.jl b/scripts/_dev/bug_issue_46/thermal_serial.jl index e9c29a6c..b94a4d09 100644 --- a/scripts/_dev/bug_issue_46/thermal_serial.jl +++ b/scripts/_dev/bug_issue_46/thermal_serial.jl @@ -67,7 +67,7 @@ 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) + stencil = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps) # Optimiser optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh; γ,γ_reinit,verbose=true,constraint_names=[:Vol]) diff --git a/scripts/_dev/full_piezo_script.jl b/scripts/_dev/full_piezo_script.jl index 7cb7c149..49dc507a 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) + stencil = HamiltonJacobiEvolution(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps) reinit!(stencil,φh,γ_reinit) ## Setup solver and FE operators diff --git a/scripts/_dev/inv_hom_block_assem_testing.jl b/scripts/_dev/inv_hom_block_assem_testing.jl index acd4c897..d8623a05 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); +stencil = 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..d0317914 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) + stencil = HamiltonJacobiEvolution(FirstOrderStencil(3,Float64),model,V_φ,el_size./order,max_steps,tol) reinit!(stencil,φ,γ_reinit) ## Setup solver and FE operators 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..13da1a00 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) + stencil = 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/mem_leak.jl b/scripts/_dev/mem_leak.jl index af9a2b28..2865affb 100644 --- a/scripts/_dev/mem_leak.jl +++ b/scripts/_dev/mem_leak.jl @@ -62,7 +62,7 @@ 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) +stencil = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps) # Optimiser optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit,verbose=true,constraint_names=[:Vol]) # Solve diff --git a/scripts/_dev/nonlinear_adjoint_MWE.jl b/scripts/_dev/nonlinear_adjoint_MWE.jl index a84705a2..082e78e0 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) + stencil = 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..c9f05eca 100644 --- a/scripts/_dev/tutorial_test/minimum_thermal_compliance_implementation.jl +++ b/scripts/_dev/tutorial_test/minimum_thermal_compliance_implementation.jl @@ -61,7 +61,7 @@ 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) +stencil = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps) # Optimiser optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit,verbose=true,constraint_names=[:Vol]) # Solve 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..f6db4233 100644 --- a/scripts/_dev/tutorial_test/minimum_thermal_compliance_implementation_3d.jl +++ b/scripts/_dev/tutorial_test/minimum_thermal_compliance_implementation_3d.jl @@ -63,7 +63,7 @@ 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) +stencil = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps) # Optimiser optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit,verbose=true,constraint_names=[:Vol]) # Solve 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..b41760df 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,7 +77,7 @@ function main() ) # Finite difference scheme scheme = FirstOrderStencil(length(el_size),Float64) - stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps) + stencil = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps) # Optimiser optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit,verbose=true,constraint_names=[:Vol]) # Solve 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..7ea63770 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,7 +78,7 @@ function main(mesh_partition,distribute) ) # Finite difference scheme scheme = FirstOrderStencil(length(el_size),Float64) - stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps) + stencil = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps) # Optimiser optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh; γ,γ_reinit,verbose=i_am_main(ranks),constraint_names=[:Vol]) 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..69c03024 100644 --- a/scripts/_dev/tutorial_test/minimum_thermal_compliance_implementation_nonlinear.jl +++ b/scripts/_dev/tutorial_test/minimum_thermal_compliance_implementation_nonlinear.jl @@ -59,7 +59,7 @@ 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) +stencil = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps) # Optimiser optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit,verbose=true,constraint_names=[:Vol]) # Solve From a8d2e265f673d33db34ead5aaa07a8e71f50e2dd Mon Sep 17 00:00:00 2001 From: "Z J Wegert (Workstation)" Date: Wed, 10 Apr 2024 11:27:40 +1000 Subject: [PATCH 07/12] Update warmup script --- compile/warmup.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/compile/warmup.jl b/compile/warmup.jl index c79ecbd1..70adbc03 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) + stencil = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,Δ./order,max_steps,tol) reinit!(stencil,φ,γ_reinit) ## Setup solver and FE operators From 6358a6febfb2ff65730a388640b5e8ca45dacd61 Mon Sep 17 00:00:00 2001 From: "Z J Wegert (Workstation)" Date: Wed, 10 Apr 2024 11:34:29 +1000 Subject: [PATCH 08/12] Started updating docs --- docs/make.jl | 2 +- .../{advection.md => levelsetevolution.md} | 8 ++++---- docs/src/tutorials/minimum_thermal_compliance.md | 14 +++++++------- 3 files changed, 12 insertions(+), 12 deletions(-) rename docs/src/reference/{advection.md => levelsetevolution.md} (64%) 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/levelsetevolution.md similarity index 64% rename from docs/src/reference/advection.md rename to docs/src/reference/levelsetevolution.md index 46166d58..a9ba77e4 100644 --- a/docs/src/reference/advection.md +++ b/docs/src/reference/levelsetevolution.md @@ -1,9 +1,9 @@ -# Advection +# LevelSetEvolution (NEED TO UPDATE THIS) -## `AdvectionStencil` +## `HamiltonJacobiEvolution` ```@docs -LevelSetTopOpt.AdvectionStencil -LevelSetTopOpt.AdvectionStencil(stencil::LevelSetTopOpt.Stencil,model,space,tol=1.e-3,max_steps=100,max_steps_reinit=2000) +LevelSetTopOpt.HamiltonJacobiEvolution +LevelSetTopOpt.HamiltonJacobiEvolution(stencil::LevelSetTopOpt.Stencil,model,space,tol=1.e-3,max_steps=100,max_steps_reinit=2000) LevelSetTopOpt.advect! LevelSetTopOpt.reinit! ``` diff --git a/docs/src/tutorials/minimum_thermal_compliance.md b/docs/src/tutorials/minimum_thermal_compliance.md index ddbb4337..5d1f4452 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) +stencil = 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 [`advect!`](@ref) and [`reinit!`](@ref) that correspond to solving the Hamilton-Jacobi evolution equation and reinitialisation equation, respectively. ### Optimiser, visualisation and IO @@ -461,7 +461,7 @@ 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) +stencil = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps) # Optimiser optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit,verbose=true,constraint_names=[:Vol]) # Solve @@ -604,7 +604,7 @@ 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) +stencil = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps) # Optimiser optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit,verbose=true,constraint_names=[:Vol]) # Solve @@ -754,7 +754,7 @@ function main() ) # Finite difference scheme scheme = FirstOrderStencil(length(el_size),Float64) - stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps) + stencil = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps) # Optimiser optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit,verbose=true,constraint_names=[:Vol]) # Solve @@ -918,7 +918,7 @@ function main(mesh_partition,distribute) ) # Finite difference scheme scheme = FirstOrderStencil(length(el_size),Float64) - stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps) + stencil = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps) # Optimiser optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh; γ,γ_reinit,verbose=i_am_main(ranks),constraint_names=[:Vol]) @@ -1070,7 +1070,7 @@ 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) +stencil = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps) # Optimiser optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit,verbose=true,constraint_names=[:Vol]) # Solve From 29e628029cbd79d5b34152bfd6a000f56adac3d8 Mon Sep 17 00:00:00 2001 From: "Z J Wegert (Workstation)" Date: Wed, 10 Apr 2024 11:50:08 +1000 Subject: [PATCH 09/12] Remove src/Advection.jl and ALM converge change --- src/Advection.jl | 501 -------------------------- src/Optimisers/AugmentedLagrangian.jl | 4 +- 2 files changed, 2 insertions(+), 503 deletions(-) delete mode 100644 src/Advection.jl diff --git a/src/Advection.jl b/src/Advection.jl deleted file mode 100644 index 241af032..00000000 --- a/src/Advection.jl +++ /dev/null @@ -1,501 +0,0 @@ -""" - abstract type Stencil - -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 -Hamilton-Jacobi evolution equation and reinitialisation equation -on order `O` finite elements in serial or parallel. - -# Parameters - -- `stencil::Stencil`: Finite difference stencil for a single step HJ - equation and reinitialisation equation. -- `model`: A `CartesianDiscreteModel`. -- `space`: FE space for level-set function -- `perm`: A permutation vector -- `params`: Tuple of additional params -- `cache`: Stencil cache -""" -struct AdvectionStencil{O} - stencil :: Stencil - model - space - perm - params - cache -end - -""" - AdvectionStencil(stencil::Stencil,model,space,tol,max_steps,max_steps_reinit) - -Create an instance of `AdvectionStencil` given a stencil, model, FE space, and -additional optional arguments. This automatically creates the DoF permutation -to handle high-order finite elements. -""" -function AdvectionStencil( - stencil::Stencil, - model, - space, - tol=1.e-3, - max_steps=100, - max_steps_reinit=2000 -) - # Parameters - order, isperiodic, Δ, ndof = get_stencil_params(model,space) - params = (;isperiodic,Δ,ndof,max_steps,max_steps_reinit,tol) - - # Dof permutation - perm = create_dof_permutation(model,space,order) - - # Caches - φ, 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 -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 - -function advect!(s::AdvectionStencil,φh,args...) - advect!(s,get_free_dof_values(φh),args...) -end - -""" - advect!(s::AdvectionStencil{O},φ,vel,γ) where O - -Solve the Hamilton-Jacobi evolution equation using the `AdvectionStencil`. - -# Hamilton-Jacobi evolution equation -``\\frac{\\partial\\phi}{\\partial t} + V(\\boldsymbol{x})\\lVert\\boldsymbol{\\nabla}\\phi\\rVert = 0,`` - -with ``\\phi(0,\\boldsymbol{x})=\\phi_0(\\boldsymbol{x})`` and ``\\boldsymbol{x}\\in D,~t\\in(0,T)``. - -# Arguments - -- `s::AdvectionStencil{O}`: Stencil for computation -- `φ`: 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 - _, _, perm_caches, stencil_cache = s.cache - Δ, isperiodic, = s.params.Δ, s.params.isperiodic - ndof, max_steps = s.params.ndof, s.params.max_steps - - _φ = (O >= 2) ? permute!(perm_caches[1],φ,s.perm) : φ - _vel = (O >= 2) ? permute!(perm_caches[2],vel,s.perm) : vel - - ## CFL Condition (requires γ≤1.0) - Δt = compute_Δt(s.stencil,Δ,γ,φ,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) - end - # Update ghost nodes - consistent!(_φ) |> fetch - end - φ = (O >= 2) ? permute_inv!(φ,_φ,s.perm) : _φ - return φ -end - -function advect!(s::AdvectionStencil{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 - - _φ = (O >= 2) ? permute!(perm_caches[1],φ,s.perm) : φ - _vel = (O >= 2) ? permute!(perm_caches[2],vel,s.perm) : vel - - ## CFL Condition (requires γ≤1.0) - Δt = compute_Δt(s.stencil,Δ,γ,φ,vel) - for _ in 1:max_steps - φ_mat = reshape(_φ,ndof) - vel_mat = reshape(_vel,ndof) - advect!(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...) -end - -""" - reinit!(s::AdvectionStencil{O},φ,γ) where O - -Solve the reinitialisation equation using the `AdvectionStencil`. - -# Reinitialisation equation -``\\frac{\\partial\\phi}{\\partial t} + \\mathrm{sign}(\\phi_0)(\\lVert\\boldsymbol{\\nabla}\\phi\\rVert-1) = 0,`` - -with ``\\phi(0,\\boldsymbol{x})=\\phi_0(\\boldsymbol{x})`` and ``\\boldsymbol{x}\\in D,~t\\in(0,T)``. - -# Arguments - -- `s::AdvectionStencil{O}`: Stencil for computation -- `φ`: level set function as a vector of degrees of freedom -- `γ`: coeffient on the time step size. -""" -function reinit!(s::AdvectionStencil{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) - - # Apply operations across partitions - step = 1; err = maximum(abs,φ); fill!(φ_tmp,0.0) - while (err > tol) && (step <= max_steps) - # Step of 1st order upwind reinitialisation equation - map(local_views(φ_tmp),local_views(_φ),local_views(vel_tmp),stencil_cache,ndof) do φ_tmp,_φ,vel_tmp,stencil_cache,S - φ_tmp_mat = reshape(φ_tmp,S) - φ_mat = reshape(_φ,S) - vel_tmp_mat = reshape(vel_tmp,S) - reinit!(s.stencil,φ_tmp_mat,φ_mat,vel_tmp_mat,Δt,Δ,isperiodic,stencil_cache) - end - - # Compute error - _φ .-= φ_tmp # φ - φ_tmp - err = maximum(abs,_φ) - step += 1 - - # Update φ - copy!(_φ,φ_tmp) - consistent!(_φ) |> fetch # We exchange ghosts here! - end - φ = (O >= 2) ? permute_inv!(φ,_φ,s.perm) : _φ - return φ -end - -function reinit!(s::AdvectionStencil{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 - - _φ = (O >= 2) ? permute!(perm_caches[1],φ,s.perm) : φ - - ## CFL Condition (requires γ≤0.5) - Δt = compute_Δt(s.stencil,Δ,γ,_φ,1.0) - - # Apply operations across partitions - step = 1; err = maximum(abs,φ); fill!(φ_tmp,0.0) - while (err > tol) && (step <= max_steps) - # Step of 1st order upwind reinitialisation equation - φ_tmp_mat = reshape(φ_tmp,ndof) - φ_mat = reshape(_φ,ndof) - vel_tmp_mat = reshape(vel_tmp,ndof) - reinit!(s.stencil,φ_tmp_mat,φ_mat,vel_tmp_mat,Δt,Δ,isperiodic,stencil_cache) - - # Compute error - _φ .-= φ_tmp - err = maximum(abs,_φ) - step += 1 - - # Update φ - copy!(_φ,φ_tmp) - end - φ = (O >= 2) ? permute_inv!(φ,_φ,s.perm) : _φ - return φ -end diff --git a/src/Optimisers/AugmentedLagrangian.jl b/src/Optimisers/AugmentedLagrangian.jl index 1c879508..6c2d48b5 100644 --- a/src/Optimisers/AugmentedLagrangian.jl +++ b/src/Optimisers/AugmentedLagrangian.jl @@ -37,7 +37,7 @@ struct AugmentedLagrangian <: Optimiser ls_evolver :: LevelSetEvolution, vel_ext :: VelocityExtension, φ0; - Λ_max = 5.0, ζ = 1.1, update_mod = 5, γ = 0.1, γ_reinit = 0.5, os_γ_mult = 0.75, + Λ_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), converged::Function = default_al_converged, debug = false, has_oscillations::Function = default_has_oscillations @@ -123,7 +123,7 @@ end function default_al_converged( m::AugmentedLagrangian; - L_tol = 0.01*maximum(get_dof_Δ(m.ls_evolver)), + L_tol = 0.01*maximum(get_dof_Δ(m.ls_evolver))/(length(get_dof_Δ(m.ls_evolver))-1), C_tol = 0.01 ) h = m.history From f142cd5e5b9b4b8d506ababf381815b5229018c9 Mon Sep 17 00:00:00 2001 From: "Z J Wegert (Workstation)" Date: Wed, 10 Apr 2024 14:29:56 +1000 Subject: [PATCH 10/12] Update docs --- docs/src/reference/levelsetevolution.md | 34 +++++++++++++------ docs/src/reference/optimisers.md | 2 +- .../tutorials/minimum_thermal_compliance.md | 2 +- .../HamiltonJacobiEvolution.jl | 8 ++--- src/LevelSetEvolution/SpatialStencil.jl | 6 ++-- 5 files changed, 31 insertions(+), 21 deletions(-) diff --git a/docs/src/reference/levelsetevolution.md b/docs/src/reference/levelsetevolution.md index a9ba77e4..3fd8ebf4 100644 --- a/docs/src/reference/levelsetevolution.md +++ b/docs/src/reference/levelsetevolution.md @@ -1,25 +1,39 @@ -# LevelSetEvolution (NEED TO UPDATE THIS) +# 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 `SpatialStencil` 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.advect! +LevelSetTopOpt.HamiltonJacobiEvolution(stencil::LevelSetTopOpt.SpatialStencil,model,space,tol=1.e-3,max_steps=100,max_steps_reinit=2000) +LevelSetTopOpt.evolve! LevelSetTopOpt.reinit! +LevelSetTopOpt.get_dof_Δ(m::HamiltonJacobiEvolution) ``` -## Stencils +## Spatial stencils for `HamiltonJacobiEvolution` ```@docs LevelSetTopOpt.FirstOrderStencil ``` -## Custom `Stencil` +## Custom `SpatialStencil` + +```@docs +LevelSetTopOpt.SpatialStencil +LevelSetTopOpt.evolve!(::LevelSetTopOpt.SpatialStencil,φ,vel,Δt,Δx,isperiodic,caches) +LevelSetTopOpt.reinit!(::LevelSetTopOpt.SpatialStencil,φ_new,φ,vel,Δt,Δx,isperiodic,caches) +LevelSetTopOpt.allocate_caches(::LevelSetTopOpt.SpatialStencil,φ,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 `SpatialStencil` types. ```@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) +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 5d1f4452..0839f6e7 100644 --- a/docs/src/tutorials/minimum_thermal_compliance.md +++ b/docs/src/tutorials/minimum_thermal_compliance.md @@ -349,7 +349,7 @@ scheme = FirstOrderStencil(length(el_size),Float64) stencil = 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 [`HamiltonJacobiEvolution`](@ref) which enables finite differencing on order `O` finite elements in serial or parallel. The [`HamiltonJacobiEvolution`](@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 diff --git a/src/LevelSetEvolution/HamiltonJacobiEvolution.jl b/src/LevelSetEvolution/HamiltonJacobiEvolution.jl index 8e86f07a..f8052495 100644 --- a/src/LevelSetEvolution/HamiltonJacobiEvolution.jl +++ b/src/LevelSetEvolution/HamiltonJacobiEvolution.jl @@ -59,7 +59,7 @@ function HamiltonJacobiEvolution( end """ - get_dof_Δ(::HamiltonJacobiEvolution) + get_dof_Δ(m::HamiltonJacobiEvolution) Return the distance betweem degree of freedom """ @@ -67,11 +67,7 @@ function get_dof_Δ(m::HamiltonJacobiEvolution) return m.params.Δ end -""" - compute_Δt(::HamiltonJacobiEvolution,φ,vel) - -Compute the time step for the `HamiltonJacobiEvolution`. -""" +# Compute the time step for the `HamiltonJacobiEvolution`. function compute_Δt(::HamiltonJacobiEvolution,Δ,γ,φ,vel) T = eltype(γ) v_norm = maximum(abs,vel) diff --git a/src/LevelSetEvolution/SpatialStencil.jl b/src/LevelSetEvolution/SpatialStencil.jl index 12dd6469..f306442f 100644 --- a/src/LevelSetEvolution/SpatialStencil.jl +++ b/src/LevelSetEvolution/SpatialStencil.jl @@ -1,10 +1,10 @@ """ abstract type SpatialStencil -Finite difference stencil for a single step of the Hamilton-Jacobi +Spatial 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. +Your own spatial stencil can be implemented by extending the methods below. """ abstract type SpatialStencil end @@ -49,7 +49,7 @@ end """ struct FirstOrderStencil{D,T} <: SpatialStencil end -A first order Godunov upwind difference scheme based on Osher and Fedkiw +A first order upwind difference scheme based on Osher and Fedkiw ([link](https://doi.org/10.1007/b98879)). """ struct FirstOrderStencil{D,T} <: SpatialStencil From ae39c6cb382b8fd7b067cd2fb503411f2fdf3cb1 Mon Sep 17 00:00:00 2001 From: zjwegert <60646897+zjwegert@users.noreply.github.com> Date: Wed, 10 Apr 2024 16:51:56 +1000 Subject: [PATCH 11/12] Rename stencil in scripts --- compile/warmup.jl | 4 +-- .../tutorials/minimum_thermal_compliance.md | 26 +++++++++---------- scripts/Benchmarks/benchmark.jl | 16 ++++++------ scripts/MPI/3d_elastic_compliance_ALM.jl | 4 +-- scripts/MPI/3d_hyperelastic_compliance_ALM.jl | 4 +-- .../3d_hyperelastic_compliance_neohook_ALM.jl | 4 +-- scripts/MPI/3d_inverse_homenisation_ALM.jl | 4 +-- scripts/MPI/3d_inverter_ALM.jl | 4 +-- scripts/MPI/3d_inverter_HPM.jl | 4 +-- .../3d_nonlinear_thermal_compliance_ALM.jl | 4 +-- scripts/MPI/3d_thermal_compliance_ALM.jl | 4 +-- scripts/MPI/elastic_compliance_ALM.jl | 4 +-- scripts/MPI/inverse_homenisation_ALM.jl | 4 +-- scripts/MPI/thermal_compliance_ALM.jl | 4 +-- scripts/MPI/thermal_compliance_ALM_PETSc.jl | 4 +-- scripts/Serial/elastic_compliance_ALM.jl | 4 +-- scripts/Serial/elastic_compliance_HPM.jl | 4 +-- scripts/Serial/elastic_compliance_LM.jl | 4 +-- scripts/Serial/hyperelastic_compliance_ALM.jl | 4 +-- .../hyperelastic_compliance_neohook_ALM.jl | 4 +-- scripts/Serial/inverse_homenisation_ALM.jl | 4 +-- scripts/Serial/inverter_ALM.jl | 4 +-- scripts/Serial/inverter_HPM.jl | 4 +-- .../nonlinear_thermal_compliance_ALM.jl | 4 +-- scripts/Serial/thermal_compliance_ALM.jl | 4 +-- .../thermal_compliance_ALM_higherorderlsf.jl | 4 +-- .../_dev/3d_hyperelastic_compliance_ALM.jl | 2 +- .../_dev/Paper_scripts/elast_comp_MPI_3D.jl | 4 +-- .../Paper_scripts/hyperelast_comp_MPI_3D.jl | 4 +-- .../_dev/Paper_scripts/inverse_hom_MPI_3D.jl | 4 +-- scripts/_dev/Paper_scripts/inverter_MPI_3D.jl | 4 +-- scripts/_dev/Paper_scripts/therm_comp_MPI.jl | 4 +-- .../_dev/Paper_scripts/therm_comp_MPI_3D.jl | 4 +-- .../_dev/Paper_scripts/therm_comp_serial.jl | 4 +-- .../MPI_benchmark_test_2D_themal.jl | 4 +-- .../MPI_benchmark_test_3D_themal.jl | 4 +-- .../benchmark_tests/serial_benchmark_test.jl | 4 +-- scripts/_dev/bug_issue_46/thermal_MPI.jl | 4 +-- scripts/_dev/bug_issue_46/thermal_serial.jl | 4 +-- scripts/_dev/full_piezo_script.jl | 4 +-- scripts/_dev/inv_hom_block_assem_testing.jl | 2 +- .../inv_hom_block_assem_testing_3D_MPI.jl | 4 +-- ...inv_hom_block_assem_testing_full_script.jl | 4 +-- scripts/_dev/mem_leak.jl | 4 +-- scripts/_dev/nonlinear_adjoint_MWE.jl | 2 +- ...nimum_thermal_compliance_implementation.jl | 4 +-- ...um_thermal_compliance_implementation_3d.jl | 4 +-- ...rmal_compliance_implementation_3d_petsc.jl | 4 +-- ..._compliance_implementation_3d_petsc_mpi.jl | 4 +-- ...mal_compliance_implementation_nonlinear.jl | 4 +-- test/seq/ThermalComplianceALMTests.jl | 4 +-- 51 files changed, 116 insertions(+), 116 deletions(-) diff --git a/compile/warmup.jl b/compile/warmup.jl index 70adbc03..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 = HamiltonJacobiEvolution(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/src/tutorials/minimum_thermal_compliance.md b/docs/src/tutorials/minimum_thermal_compliance.md index 0839f6e7..3d52db3c 100644 --- a/docs/src/tutorials/minimum_thermal_compliance.md +++ b/docs/src/tutorials/minimum_thermal_compliance.md @@ -346,7 +346,7 @@ 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 = HamiltonJacobiEvolution(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 [`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. @@ -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 = HamiltonJacobiEvolution(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 = HamiltonJacobiEvolution(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 = HamiltonJacobiEvolution(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 = HamiltonJacobiEvolution(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 = HamiltonJacobiEvolution(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 759bf6e0..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 = HamiltonJacobiEvolution(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 = HamiltonJacobiEvolution(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 = HamiltonJacobiEvolution(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 = HamiltonJacobiEvolution(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 d44fdff3..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 = HamiltonJacobiEvolution(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 c103beab..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 = HamiltonJacobiEvolution(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 1a9fe02e..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 = HamiltonJacobiEvolution(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 2fed57d8..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 = HamiltonJacobiEvolution(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 5d270952..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 = HamiltonJacobiEvolution(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 f1684cf0..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 = HamiltonJacobiEvolution(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 126309d5..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 = HamiltonJacobiEvolution(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 e0cea609..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 = HamiltonJacobiEvolution(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 6fa7a17c..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 = HamiltonJacobiEvolution(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 ec4d8886..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 = HamiltonJacobiEvolution(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 bce212f1..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 = HamiltonJacobiEvolution(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 cfe2acd6..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 = HamiltonJacobiEvolution(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 7cd20e1a..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 = HamiltonJacobiEvolution(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 50964928..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 = HamiltonJacobiEvolution(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 0517ef1b..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 = HamiltonJacobiEvolution(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 c455de7b..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 = HamiltonJacobiEvolution(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 03ea9f55..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 = HamiltonJacobiEvolution(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 c079c08c..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 = HamiltonJacobiEvolution(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 2774a9ea..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 = HamiltonJacobiEvolution(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 f2587eec..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 = HamiltonJacobiEvolution(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 51597244..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 = HamiltonJacobiEvolution(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 6ddf66cc..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 = HamiltonJacobiEvolution(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 0b076760..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 = HamiltonJacobiEvolution(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 de79e6fc..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 = HamiltonJacobiEvolution(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 c903951a..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 = HamiltonJacobiEvolution(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 6f463318..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 = HamiltonJacobiEvolution(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 657b35c4..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 = HamiltonJacobiEvolution(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 75163cb0..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 = HamiltonJacobiEvolution(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 8320ff76..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 = HamiltonJacobiEvolution(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 3cb4f8ba..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 = HamiltonJacobiEvolution(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 9e7e2342..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 = HamiltonJacobiEvolution(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 2d09b546..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 = HamiltonJacobiEvolution(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 5842ead6..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 = HamiltonJacobiEvolution(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 0f8a784a..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 = HamiltonJacobiEvolution(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 176e7b05..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 = HamiltonJacobiEvolution(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 b94a4d09..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 = HamiltonJacobiEvolution(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 49dc507a..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 = HamiltonJacobiEvolution(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 d8623a05..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 = HamiltonJacobiEvolution(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 d0317914..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 = HamiltonJacobiEvolution(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 13da1a00..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 = HamiltonJacobiEvolution(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 2865affb..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 = HamiltonJacobiEvolution(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 082e78e0..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 = HamiltonJacobiEvolution(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 c9f05eca..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 = HamiltonJacobiEvolution(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 f6db4233..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 = HamiltonJacobiEvolution(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 b41760df..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 = HamiltonJacobiEvolution(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 7ea63770..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 = HamiltonJacobiEvolution(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 69c03024..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 = HamiltonJacobiEvolution(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/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 From 7603a8934548635565f36748b90e7f0a9d064564 Mon Sep 17 00:00:00 2001 From: zjwegert <60646897+zjwegert@users.noreply.github.com> Date: Wed, 10 Apr 2024 17:08:32 +1000 Subject: [PATCH 12/12] Undo Stencil naming --- docs/src/reference/levelsetevolution.md | 16 +++++----- .../HamiltonJacobiEvolution.jl | 14 ++++---- src/LevelSetEvolution/LevelSetEvolution.jl | 2 +- .../{SpatialStencil.jl => Stencil.jl} | 32 +++++++++---------- 4 files changed, 32 insertions(+), 32 deletions(-) rename src/LevelSetEvolution/{SpatialStencil.jl => Stencil.jl} (89%) diff --git a/docs/src/reference/levelsetevolution.md b/docs/src/reference/levelsetevolution.md index 3fd8ebf4..3f41a926 100644 --- a/docs/src/reference/levelsetevolution.md +++ b/docs/src/reference/levelsetevolution.md @@ -1,12 +1,12 @@ # 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 `SpatialStencil` 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. +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.SpatialStencil,model,space,tol=1.e-3,max_steps=100,max_steps_reinit=2000) +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) @@ -18,18 +18,18 @@ LevelSetTopOpt.get_dof_Δ(m::HamiltonJacobiEvolution) LevelSetTopOpt.FirstOrderStencil ``` -## Custom `SpatialStencil` +## Custom `Stencil` ```@docs -LevelSetTopOpt.SpatialStencil -LevelSetTopOpt.evolve!(::LevelSetTopOpt.SpatialStencil,φ,vel,Δt,Δx,isperiodic,caches) -LevelSetTopOpt.reinit!(::LevelSetTopOpt.SpatialStencil,φ_new,φ,vel,Δt,Δx,isperiodic,caches) -LevelSetTopOpt.allocate_caches(::LevelSetTopOpt.SpatialStencil,φ,vel) +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 `SpatialStencil` types. +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 diff --git a/src/LevelSetEvolution/HamiltonJacobiEvolution.jl b/src/LevelSetEvolution/HamiltonJacobiEvolution.jl index f8052495..06dba88b 100644 --- a/src/LevelSetEvolution/HamiltonJacobiEvolution.jl +++ b/src/LevelSetEvolution/HamiltonJacobiEvolution.jl @@ -9,16 +9,16 @@ Based on the scheme by Osher and Fedkiw ([link](https://doi.org/10.1007/b98879)) # Parameters -- `stencil::SpatialStencil`: Spatial 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 - `perm`: A permutation vector - `params`: Tuple of additional params -- `cache`: SpatialStencil cache +- `cache`: Stencil cache """ struct HamiltonJacobiEvolution{O} <: LevelSetEvolution - stencil :: SpatialStencil + stencil :: Stencil model space perm @@ -27,14 +27,14 @@ struct HamiltonJacobiEvolution{O} <: LevelSetEvolution end """ - HamiltonJacobiEvolution(stencil::SpatialStencil,model,space,tol,max_steps,max_steps_reinit) + HamiltonJacobiEvolution(stencil::Stencil,model,space,tol,max_steps,max_steps_reinit) 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 HamiltonJacobiEvolution( - stencil::SpatialStencil, + stencil::Stencil, model, space, tol=1.e-3, @@ -302,7 +302,7 @@ function PartitionedArrays.permute_indices(indices::LocalIndices,perm) return LocalIndices(n_glob,id,l2g,l2o) end -function allocate_caches(s::SpatialStencil,φ::Vector,vel::Vector,perm,order,ndofs) +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) @@ -310,7 +310,7 @@ function allocate_caches(s::SpatialStencil,φ::Vector,vel::Vector,perm,order,ndo return φ_tmp, vel_tmp, perm_caches, stencil_caches end -function allocate_caches(s::SpatialStencil,φ::PVector,vel::PVector,perm,order,local_ndofs) +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 diff --git a/src/LevelSetEvolution/LevelSetEvolution.jl b/src/LevelSetEvolution/LevelSetEvolution.jl index d6a743f4..6f9c4f32 100644 --- a/src/LevelSetEvolution/LevelSetEvolution.jl +++ b/src/LevelSetEvolution/LevelSetEvolution.jl @@ -35,5 +35,5 @@ function get_dof_Δ(::LevelSetEvolution) @abstractmethod end -include("SpatialStencil.jl") +include("Stencil.jl") include("HamiltonJacobiEvolution.jl") \ No newline at end of file diff --git a/src/LevelSetEvolution/SpatialStencil.jl b/src/LevelSetEvolution/Stencil.jl similarity index 89% rename from src/LevelSetEvolution/SpatialStencil.jl rename to src/LevelSetEvolution/Stencil.jl index f306442f..4b15e82e 100644 --- a/src/LevelSetEvolution/SpatialStencil.jl +++ b/src/LevelSetEvolution/Stencil.jl @@ -1,58 +1,58 @@ """ - abstract type SpatialStencil + abstract type Stencil -Spatial finite difference stencil for a single step of the Hamilton-Jacobi +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 SpatialStencil end +abstract type Stencil end """ - allocate_caches(::SpatialStencil,φ,vel) + allocate_caches(::Stencil,φ,vel) -Allocate caches for a given `SpatialStencil`. +Allocate caches for a given `Stencil`. """ -function allocate_caches(::SpatialStencil,φ,vel) +function allocate_caches(::Stencil,φ,vel) nothing # By default, no caches are required. end """ - check_order(::SpatialStencil,order) + check_order(::Stencil,order) Throw error if insufficient reference element order to implement stencil in parallel. """ -function check_order(::SpatialStencil,order) +function check_order(::Stencil,order) @abstractmethod end """ - reinit!(::SpatialStencil,φ_new,φ,vel,Δt,Δx,isperiodic,caches) -> φ + reinit!(::Stencil,φ_new,φ,vel,Δt,Δx,isperiodic,caches) -> φ -Single finite difference step of the reinitialisation equation for a given `SpatialStencil`. +Single finite difference step of the reinitialisation equation for a given `Stencil`. """ -function reinit!(::SpatialStencil,φ_new,φ,vel,Δt,Δx,isperiodic,caches) +function reinit!(::Stencil,φ_new,φ,vel,Δt,Δx,isperiodic,caches) @abstractmethod end """ - evolve!(::SpatialStencil,φ,vel,Δt,Δx,isperiodic,caches) -> φ + evolve!(::Stencil,φ,vel,Δt,Δx,isperiodic,caches) -> φ Single finite difference step of the Hamilton-Jacobi evoluation equation for a given -`SpatialStencil`. +`Stencil`. """ -function evolve!(::SpatialStencil,φ,vel,Δt,Δx,isperiodic,caches) +function evolve!(::Stencil,φ,vel,Δt,Δx,isperiodic,caches) @abstractmethod end """ - struct FirstOrderStencil{D,T} <: SpatialStencil 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} <: SpatialStencil +struct FirstOrderStencil{D,T} <: Stencil function FirstOrderStencil(D::Integer,::Type{T}) where T<:Real new{D,T}() end