Merge pull request #55 from zjwegert/advection-refactor

Advection refactor

compile/warmup.jl

@@ -70,7 +70,7 @@ function main(mesh_partition,distribute)
   dJ = (q,u,φ,dΩ,dΓ_N) -> ∫((ξ-D*∇(u)⋅∇(u))*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ;
 
   ## Finite difference solver and level set function
-  stencil = AdvectionStencil(FirstOrderStencil(2,Float64),model,V_φ,Δ./order,max_steps,tol)
+  ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,Δ./order,max_steps,tol)
   reinit!(stencil,φ,γ_reinit)
 
   ## Setup solver and FE operators
@@ -84,7 +84,7 @@ function main(mesh_partition,distribute)
 
   ## Optimiser
   _conv_cond = t->LevelSetTopOpt.conv_cond(t;coef=1/50);
-  optimiser = AugmentedLagrangian(φ,pcfs,stencil,vel_ext,interp,el_size,γ,γ_reinit,conv_criterion=_conv_cond);
+  optimiser = AugmentedLagrangian(φ,pcfs,ls_evo,vel_ext,interp,el_size,γ,γ_reinit,conv_criterion=_conv_cond);
   for history in optimiser
     it,Ji,_,_ = last(history)
     print_history(it,["J"=>Ji];ranks=ranks)

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",

docs/src/reference/levelsetevolution.md

@@ -0,0 +1,39 @@
+# LevelSetEvolution
+In LevelSetTopOpt, the level set is evolved and reinitialised using a `LevelSetEvolution` method. The most standard of these is the Hamilton-Jacobi evolution equation solved using a first order upwind finite difference scheme. A forward Euler in time method is provided below via `HamiltonJacobiEvolution <: LevelSetEvolution` along with an upwind finite difference stencil for the spatial discretisation via `FirstOrderStencil`.
+
+This can be extended in several ways. For example, higher order spatial stencils can be implemented by extending the `Stencil` interface below. In addition, more advanced ODE solvers could be implemented (e.g., Runge–Kutta methods) or entirely different level set evolution methods by extending the `LevelSetEvolution` interface below.
+
+## `HamiltonJacobiEvolution`
+```@docs
+LevelSetTopOpt.HamiltonJacobiEvolution
+LevelSetTopOpt.HamiltonJacobiEvolution(stencil::LevelSetTopOpt.Stencil,model,space,tol=1.e-3,max_steps=100,max_steps_reinit=2000)
+LevelSetTopOpt.evolve!
+LevelSetTopOpt.reinit!
+LevelSetTopOpt.get_dof_Δ(m::HamiltonJacobiEvolution)
+```
+
+## Spatial stencils for `HamiltonJacobiEvolution`
+
+```@docs
+LevelSetTopOpt.FirstOrderStencil
+```
+
+## Custom `Stencil`
+
+```@docs
+LevelSetTopOpt.Stencil
+LevelSetTopOpt.evolve!(::LevelSetTopOpt.Stencil,φ,vel,Δt,Δx,isperiodic,caches)
+LevelSetTopOpt.reinit!(::LevelSetTopOpt.Stencil,φ_new,φ,vel,Δt,Δx,isperiodic,caches)
+LevelSetTopOpt.allocate_caches(::LevelSetTopOpt.Stencil,φ,vel)
+LevelSetTopOpt.check_order
+```
+
+## Custom `LevelSetEvolution`
+To implement a custom level set evolution method, we can extend the methods below. For example, one could consider Reaction-Diffusion-based evolution of the level set function. This can be solved with a finite element method and so we can implement a new type that inherits from `LevelSetEvolution` independently of the `Stencil` types.
+
+```@docs
+LevelSetTopOpt.LevelSetEvolution
+LevelSetTopOpt.evolve!(::LevelSetTopOpt.LevelSetEvolution,φ,args...)
+LevelSetTopOpt.reinit!(::LevelSetTopOpt.LevelSetEvolution,φ,args...)
+LevelSetTopOpt.get_dof_Δ(::LevelSetTopOpt.LevelSetEvolution)
+```

docs/src/reference/optimisers.md

@@ -27,7 +27,7 @@ LevelSetTopOpt.OptimiserHistory
 LevelSetTopOpt.OptimiserHistorySlice
 ```
 
-## Custom optimiser
+## Custom `Optimiser`
 ```@docs
 LevelSetTopOpt.Optimiser
 LevelSetTopOpt.iterate(::LevelSetTopOpt.Optimiser)

docs/src/tutorials/minimum_thermal_compliance.md

@@ -346,18 +346,18 @@ Both of these equations can be solved numerically on a Cartesian mesh using a fi
 ```julia
 # Finite difference scheme
 scheme = FirstOrderStencil(length(el_size),Float64)
-stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps)
+ls_evo = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps)
 ```
 
-In the above we first build an object [`FirstOrderStencil`](@ref) that represents a finite difference stencil for a single step of the Hamilton-Jacobi evolution equation and reinitialisation equation. We use `length(el_size)` to indicate the dimension of the problem. We then create an [`AdvectionStencil`](@ref) which enables finite differencing on order `O` finite elements in serial or parallel. The [`AdvectionStencil`](@ref) object provides two important methods [`advect!`](@ref) and [`reinit!`](@ref) that correspond to solving the Hamilton-Jacobi evolution equation and reinitialisation equation, respectively.
+In the above we first build an object [`FirstOrderStencil`](@ref) that represents a finite difference stencil for a single step of the Hamilton-Jacobi evolution equation and reinitialisation equation. We use `length(el_size)` to indicate the dimension of the problem. We then create an [`HamiltonJacobiEvolution`](@ref) which enables finite differencing on order `O` finite elements in serial or parallel. The [`HamiltonJacobiEvolution`](@ref) object provides two important methods [`evolve!`](@ref) and [`reinit!`](@ref) that correspond to solving the Hamilton-Jacobi evolution equation and reinitialisation equation, respectively.
 
 ### Optimiser, visualisation and IO
 
 We may now create the optimiser object. This structure holds all information regarding the optimisation problem that we wish to solve and implements an optimisation algorithm as a Julia [iterator](https://docs.julialang.org/en/v1/manual/interfaces/). For the purpose of this tutorial we use a standard augmented Lagrangian method based on [6]. In our script, we create an instance of the [`AugmentedLagrangian`](@ref) via
 
 ```julia
 # Optimiser
-optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit,verbose=true,constraint_names=[:Vol])
+optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;γ,γ_reinit,verbose=true,constraint_names=[:Vol])
 ```
 
 As optimisers inheriting from [`LevelSetTopOpt.Optimiser`](@ref) implement Julia's iterator functionality, we can solve the optimisation problem to convergence by iterating over the optimiser:
@@ -461,9 +461,9 @@ a_hilb(p,q) =∫(α^2*∇(p)⋅∇(q) + p*q)dΩ
 vel_ext = VelocityExtension(a_hilb,U_reg,V_reg)
 # Finite difference scheme
 scheme = FirstOrderStencil(length(el_size),Float64)
-stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps)
+ls_evo = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps)
 # Optimiser
-optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit,verbose=true,constraint_names=[:Vol])
+optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;γ,γ_reinit,verbose=true,constraint_names=[:Vol])
 # Solve
 for (it,uh,φh) in optimiser
   data = ["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh]
@@ -604,9 +604,9 @@ a_hilb(p,q) =∫(α^2*∇(p)⋅∇(q) + p*q)dΩ
 vel_ext = VelocityExtension(a_hilb,U_reg,V_reg)
 # Finite difference scheme
 scheme = FirstOrderStencil(length(el_size),Float64)
-stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps)
+ls_evo = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps)
 # Optimiser
-optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit,verbose=true,constraint_names=[:Vol])
+optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;γ,γ_reinit,verbose=true,constraint_names=[:Vol])
 # Solve
 for (it,uh,φh) in optimiser
   data = ["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh]
@@ -754,9 +754,9 @@ function main()
   )
   # Finite difference scheme
   scheme = FirstOrderStencil(length(el_size),Float64)
-  stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps)
+  ls_evo = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps)
   # Optimiser
-  optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit,verbose=true,constraint_names=[:Vol])
+  optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;γ,γ_reinit,verbose=true,constraint_names=[:Vol])
   # Solve
   for (it,uh,φh) in optimiser
     data = ["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh]
@@ -827,7 +827,7 @@ We then adjust lines 28-31 as follows:
 The function `i_am_main` returns true only on the first processor. This function is useful for ensuring certain operations only happen once instead of several times across each executable. In addition, we now create a partitioned Cartesian model using `CartesianDiscreteModel(ranks,mesh_partition,dom,el_size)`. Finally, we adjust line 82 to ensure that verbosity only happens on the first processors:
 
 ```julia
-optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;
+optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;
   γ,γ_reinit,verbose=i_am_main(ranks),constraint_names=[:Vol])
 ```
 
@@ -918,9 +918,9 @@ function main(mesh_partition,distribute)
   )
   # Finite difference scheme
   scheme = FirstOrderStencil(length(el_size),Float64)
-  stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps)
+  ls_evo = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps)
   # Optimiser
-  optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;
+  optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;
     γ,γ_reinit,verbose=i_am_main(ranks),constraint_names=[:Vol])
   # Solve
   for (it,uh,φh) in optimiser
@@ -1070,9 +1070,9 @@ a_hilb(p,q) =∫(α^2*∇(p)⋅∇(q) + p*q)dΩ
 vel_ext = VelocityExtension(a_hilb,U_reg,V_reg)
 # Finite difference scheme
 scheme = FirstOrderStencil(length(el_size),Float64)
-stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps)
+ls_evo = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps)
 # Optimiser
-optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit,verbose=true,constraint_names=[:Vol])
+optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;γ,γ_reinit,verbose=true,constraint_names=[:Vol])
 # Solve
 for (it,uh,φh) in optimiser
   data = ["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh]

scripts/Benchmarks/benchmark.jl

@@ -102,9 +102,9 @@ function nl_elast(mesh_partition,ranks,el_size,order,verbose)
   )
   # Finite difference scheme
   scheme = FirstOrderStencil(length(el_size),Float64)
-  stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps)
+  ls_evo = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps)
   # Optimiser
-  return AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit,verbose)
+  return AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;γ,γ_reinit,verbose)
 end
 
 function therm(mesh_partition,ranks,el_size,order,verbose)
@@ -180,9 +180,9 @@ function therm(mesh_partition,ranks,el_size,order,verbose)
   )
   # Finite difference scheme
   scheme = FirstOrderStencil(length(el_size),Float64)
-  stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps)
+  ls_evo = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps)
   # Optimiser
-  return AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit,verbose)
+  return AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;γ,γ_reinit,verbose)
 end
 
 function elast(mesh_partition,ranks,el_size,order,verbose)
@@ -256,9 +256,9 @@ function elast(mesh_partition,ranks,el_size,order,verbose)
   )
   # Finite difference scheme
   scheme = FirstOrderStencil(length(el_size),Float64)
-  stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps)
+  ls_evo = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps)
   # Optimiser
-  return AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit,verbose)
+  return AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;γ,γ_reinit,verbose)
 end
 
 function inverter_HPM(mesh_partition,ranks,el_size,order,verbose)
@@ -329,7 +329,7 @@ function inverter_HPM(mesh_partition,ranks,el_size,order,verbose)
   UΓ_out(u,φ,dΩ,dΓ_in,dΓ_out) = ∫((u⋅-e₁-δₓ)/vol_Γ_out)dΓ_out
 
   ## Finite difference solver
-  stencil = AdvectionStencil(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps)
+  ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps)
 
   ## Setup solver and FE operators
   Tm = SparseMatrixCSR{0,PetscScalar,PetscInt}
@@ -355,7 +355,7 @@ function inverter_HPM(mesh_partition,ranks,el_size,order,verbose)
   )
 
   ## Optimiser
-  return HilbertianProjection(pcfs,stencil,vel_ext,φh;γ,γ_reinit,verbose=verbose)
+  return HilbertianProjection(pcfs,ls_evo,vel_ext,φh;γ,γ_reinit,verbose=verbose)
 end
 
 with_mpi() do distribute

scripts/MPI/3d_elastic_compliance_ALM.jl

@@ -81,7 +81,7 @@ function main(mesh_partition,distribute,el_size)
   dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
 
   ## Finite difference solver and level set function
-  stencil = AdvectionStencil(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps)
+  ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps)
 
   ## Setup solver and FE operators
   Tm = SparseMatrixCSR{0,PetscScalar,PetscInt}
@@ -107,7 +107,7 @@ function main(mesh_partition,distribute,el_size)
   )
 
   ## Optimiser
-  optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;
+  optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;
     γ,γ_reinit,verbose=i_am_main(ranks),constraint_names=[:Vol])
   for (it, uh, φh) in optimiser
     data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh]

scripts/MPI/3d_hyperelastic_compliance_ALM.jl

@@ -83,7 +83,7 @@ function main(mesh_partition,distribute,el_size)
   dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
 
   ## Finite difference solver and level set function
-  stencil = AdvectionStencil(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps)
+  ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps)
 
   ## Setup solver and FE operators
   Tm = SparseMatrixCSR{0,PetscScalar,PetscInt}
@@ -110,7 +110,7 @@ function main(mesh_partition,distribute,el_size)
   )
 
   ## Optimiser
-  optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;
+  optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;
     γ,γ_reinit,verbose=i_am_main(ranks),constraint_names=[:Vol])
   for (it, uh, φh) in optimiser
     data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh]

scripts/MPI/3d_hyperelastic_compliance_neohook_ALM.jl

@@ -97,7 +97,7 @@ function main(mesh_partition,distribute,el_size)
   dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
 
   ## Finite difference solver and level set function
-  stencil = AdvectionStencil(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps)
+  ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps)
 
   ## Setup solver and FE operators
   Tm = SparseMatrixCSR{0,PetscScalar,PetscInt}
@@ -124,7 +124,7 @@ function main(mesh_partition,distribute,el_size)
   )
 
   ## Optimiser
-  optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;
+  optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;
     γ,γ_reinit,verbose=i_am_main(ranks),constraint_names=[:Vol])
   for (it, uh, φh) in optimiser
     data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh]

scripts/MPI/3d_inverse_homenisation_ALM.jl

@@ -85,7 +85,7 @@ function main(mesh_partition,distribute,el_size)
   dVol(q,u,φ,dΩ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
 
   ## Finite difference solver and level set function
-  stencil = AdvectionStencil(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps)
+  ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps)
 
   ## Setup solver and FE operators
   Tm = SparseMatrixCSR{0,PetscScalar,PetscInt}
@@ -111,7 +111,7 @@ function main(mesh_partition,distribute,el_size)
   )
 
   ## Optimiser
-  optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;
+  optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;
     γ,γ_reinit,verbose=i_am_main(ranks),constraint_names=[:Vol])
   for (it, uh, φh) in optimiser
     data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh))]

scripts/MPI/3d_inverter_ALM.jl

@@ -99,7 +99,7 @@ function main(mesh_partition,distribute,el_size)
   UΓ_out(u,φ,dΩ,dΓ_in,dΓ_out) = ∫((u⋅-e₁-δₓ)/vol_Γ_out)dΓ_out
 
   ## Finite difference solver and level set function
-  stencil = AdvectionStencil(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps)
+  ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps)
 
   ## Setup solver and FE operators
   Tm = SparseMatrixCSR{0,PetscScalar,PetscInt}
@@ -125,7 +125,7 @@ function main(mesh_partition,distribute,el_size)
   )
 
   ## Optimiser
-  optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;
+  optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;
     γ,γ_reinit,verbose=i_am_main(ranks),constraint_names=[:Vol,:UΓ_out])
   for (it, uh, φh) in optimiser
     data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh]

scripts/MPI/3d_inverter_HPM.jl

@@ -99,7 +99,7 @@ function main(mesh_partition,distribute,el_size)
   UΓ_out(u,φ,dΩ,dΓ_in,dΓ_out) = ∫((u⋅-e₁-δₓ)/vol_Γ_out)dΓ_out
 
   ## Finite difference solver and level set function
-  stencil = AdvectionStencil(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps)
+  ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps)
 
   ## Setup solver and FE operators
   Tm = SparseMatrixCSR{0,PetscScalar,PetscInt}
@@ -126,7 +126,7 @@ function main(mesh_partition,distribute,el_size)
 
   ## Optimiser
   ls_enabled=false # Setting to true will use a line search instead of oscillation detection
-  optimiser = HilbertianProjection(pcfs,stencil,vel_ext,φh;γ,γ_reinit,ls_enabled,α_min=0.7,
+  optimiser = HilbertianProjection(pcfs,ls_evo,vel_ext,φh;γ,γ_reinit,ls_enabled,α_min=0.7,
     ls_γ_max=0.05,verbose=i_am_main(ranks),constraint_names=[:Vol,:UΓ_out])
   for (it, uh, φh) in optimiser
     data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh]

scripts/MPI/3d_nonlinear_thermal_compliance_ALM.jl

@@ -85,7 +85,7 @@ function main(mesh_partition,distribute,el_size)
   dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
 
   ## Finite difference solver and level set function
-  stencil = AdvectionStencil(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps)
+  ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps)
 
   ## Setup solver and FE operators
   Tm = SparseMatrixCSR{0,PetscScalar,PetscInt}
@@ -113,7 +113,7 @@ function main(mesh_partition,distribute,el_size)
   )
 
   ## Optimiser
-  optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;
+  optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;
     γ,γ_reinit,verbose=i_am_main(ranks),constraint_names=[:Vol])
   for (it, uh, φh) in optimiser
     data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh]

scripts/MPI/3d_thermal_compliance_ALM.jl

@@ -81,7 +81,7 @@ function main(mesh_partition,distribute,el_size)
   dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
 
   ## Finite difference solver and level set function
-  stencil = AdvectionStencil(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps)
+  ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps)
 
   ## Setup solver and FE operators
   Tm = SparseMatrixCSR{0,PetscScalar,PetscInt}
@@ -107,7 +107,7 @@ function main(mesh_partition,distribute,el_size)
   )
 
   ## Optimiser
-  optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;
+  optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;
     γ,γ_reinit,verbose=i_am_main(ranks),constraint_names=[:Vol])
   for (it, uh, φh) in optimiser
     data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh]