Added custom body force.

src/Flow.jl

@@ -71,6 +71,13 @@ end
 accelerate!(r,dt,g::Nothing,U::Function) = accelerate!(r,dt,(i,t)->ForwardDiff.derivative(τ->U(i,τ),t),())
 accelerate!(r,dt,g::Function,U::Function) = accelerate!(r,dt,(i,t)->g(i,t)+ForwardDiff.derivative(τ->U(i,τ),t),())
 accelerate!(r,dt,::Nothing,::Tuple) = nothing
+"""
+    body_force!(r,force)
+
+Adds a body force to the RHS
+"""
+body_force!(_,::Nothing) = nothing
+body_force!(r,force::AbstractArray) = r .+= force
 """
     BCTuple(U,dt,N)
 
@@ -150,18 +157,20 @@ end
 Integrate the `Flow` one time step using the [Boundary Data Immersion Method](https://eprints.soton.ac.uk/369635/)
 and the `AbstractPoisson` pressure solver to project the velocity onto an incompressible flow.
 """
-@fastmath function mom_step!(a::Flow{N},b::AbstractPoisson; CFL_f=CFL) where N
+@fastmath function mom_step!(a::Flow{N},b::AbstractPoisson; bf=nothing, CFL_f=CFL) where N
     a.u⁰ .= a.u; scale_u!(a,0)
     # predictor u → u'
     U = BCTuple(a.U,@view(a.Δt[1:end-1]),N)
     conv_diff!(a.f,a.u⁰,a.σ,ν=a.ν,perdir=a.perdir)
+    body_force!(a.f,bf)
     accelerate!(a.f,@view(a.Δt[1:end-1]),a.g,a.U)
     BDIM!(a); BC!(a.u,U,a.exitBC,a.perdir)
     a.exitBC && exitBC!(a.u,a.u⁰,U,a.Δt[end]) # convective exit
     project!(a,b); BC!(a.u,U,a.exitBC,a.perdir)
     # corrector u → u¹
     U = BCTuple(a.U,a.Δt,N)
     conv_diff!(a.f,a.u,a.σ,ν=a.ν,perdir=a.perdir)
+    body_force!(a.f,bf)
     accelerate!(a.f,a.Δt,a.g,a.U)
     BDIM!(a); scale_u!(a,0.5); BC!(a.u,U,a.exitBC,a.perdir)
     project!(a,b,0.5); BC!(a.u,U,a.exitBC,a.perdir)

src/WaterLily.jl

@@ -94,17 +94,17 @@ Integrate the simulation `sim` up to dimensionless time `t_end`.
 If `remeasure=true`, the body is remeasured at every time step.
 Can be set to `false` for static geometries to speed up simulation.
 """
-function sim_step!(sim::Simulation,t_end;remeasure=true,max_steps=typemax(Int),verbose=false,CFL_f=CFL)
+function sim_step!(sim::Simulation,t_end;remeasure=true,max_steps=typemax(Int),bf=nothing,CFL_f=CFL,verbose=false)
     steps₀ = length(sim.flow.Δt)
     while sim_time(sim) < t_end && length(sim.flow.Δt) - steps₀ < max_steps
         sim_step!(sim; remeasure, CFL_f)
         verbose && println("tU/L=",round(sim_time(sim),digits=4),
             ", Δt=",round(sim.flow.Δt[end],digits=3))
     end
 end
-function sim_step!(sim::Simulation;remeasure=true,CFL_f=CFL)
+function sim_step!(sim::Simulation;remeasure=true,bf=nothing,CFL_f=CFL)
     remeasure && measure!(sim)
-    mom_step!(sim.flow,sim.pois; CFL_f)
+    mom_step!(sim.flow,sim.pois; bf, CFL_f)
 end
 
 """