Setting up a time-varying PDE

In this example we will demonstrate the use of the package on a time-dependent PDE, a problem of unsteady heat conduction. We will use the package to solve for the interior diffusion of temperature from a circle held at constant temperature.

We seek to solve the heat conduction equation with Dirichlet boundary conditions

\[\dfrac{\partial T}{\partial t} = \kappa \nabla^2 T + q + \delta(\chi) \sigma - \kappa \nabla\cdot \left( \delta(\chi) \mathbf{n} [T] \right)\]

subject to $T = T_b$ on the immersed surface. We might be solving this external to a surface, or it might be internal. The quantity $\sigma$ is the Lagrange multiplier. In this context, it is the heat flux through the surface.

In the spatially discrete formulation, the problem takes the form

\[\begin{bmatrix} \mathcal{L}_C^\kappa & R_C \\ R_C^T & 0 \end{bmatrix}\begin{pmatrix} T \\ -\sigma \end{pmatrix} = \begin{pmatrix} q - \kappa D_s [T] \\ (T^+_b + T^-_b)/2 \end{pmatrix}\]

where $\mathcal{L}_C^\kappa = \mathrm{d}/\mathrm{d}t - \kappa L_C$, where $[T] = T_b^+ - T_b^-$ is the jump in temperature across the surface. As in the time-independent problems, we can specify whether we are solving it external or internal to a surface by setting the boundary value to zero in the other region. However, in contrast to the time-independent problems, we have to advance this problem in time. The system above has the form of a constrained ODE system, which the ConstrainedSystems.jl package treats. We will make use of this package in the example below.

To support this, there are a few additional steps in our setup of the problem:

• we (as the implementers of the PDE) need to specify the functions that calculate the various parts of this constrained ODE system.
• we (as the users of this implementation) need to specify the time step size, the initial conditions, the time integration range, and create the integrator to advance the solution.

The latter of these is very easy, as we'll find. Most of our attention will be on the first part: how to set up the constrained ODE system. For this, we will make use of the ODEFunctionList, which assembles the various functions and operators into a ConstrainedODEFunction, to be used by the ConstrainedSystems.jl package.

using ImmersedLayers
using Plots
using UnPack

Set up the constrained ODE system operators

The problem type is generated with the usual macro call. In this example, we will make use of more of the capabilities of the resulting problem constructor for "packing" it with information about the problem.

@ilmproblem DirichletHeatConduction scalar

The constrained ODE system requires us to provide functions that calculate the RHS of the ODE, the RHS of the constraint equation, the Lagrange multiplier force term in the ODE, and the action of the boundary operator on the state vector. (You can see the generic form of the system by typing ?ConstrainedODEFunction) As you will see, in this example these are in-place operators: their first argument holds the result, which is changed (i.e., mutated) by the function. Below, we construct the function that calculates the RHS of the heat conduction ODE. We have omitted the volumetric heat flux here, supplying only the double-layer term. Note how this makes use of the physical parameters in phys_params and the boundary data via functions in bc. The functions for the boundary data supply the boundary values. Also, note that the function returns dT in the first argument. This represents this function's contribution to $dT/dt$. The argument x isn't used here, but would generally hold information about the body state.

function heatconduction_ode_rhs!(dT,T,x,sys::ILMSystem,t)
    @unpack bc, forcing, phys_params, extra_cache, base_cache = sys
    @unpack dTb = extra_cache

    κ = phys_params["diffusivity"]

    # Calculate the double-layer term
    prescribed_surface_jump!(dTb,x,t,sys)
    surface_divergence!(dT,-κ*dTb,sys)

    return dT
end

heatconduction_ode_rhs! (generic function with 1 method)

Now, we create the function that calculates the RHS of the boundary condition. For this Dirichlet condition, we simply take the average of the interior and exterior prescribed values. The first argument dTb holds the result. Again, x isn't used here.

function heatconduction_bc_rhs!(dTb,x,sys::ILMSystem,t)
    prescribed_surface_average!(dTb,x,t,sys)
    return dTb
end

heatconduction_bc_rhs! (generic function with 1 method)

This function calculates the contribution to $dT/dt$ from the Lagrange multiplier (the input σ). Here, we simply regularize the negative of this to the grid.

function heatconduction_constraint_force!(dT,σ,x,sys::ILMSystem)
    regularize!(dT,-σ,sys)
    return dT
end

heatconduction_constraint_force! (generic function with 1 method)

Now, we provide the transpose term of the previous function: a function that interpolates the temperature (state vector) onto the boundary. The first argument dTb holds the result.

function heatconduction_bc_op!(dTb,T,x,sys::ILMSystem)
    interpolate!(dTb,T,sys)
    return dTb
end

heatconduction_bc_op! (generic function with 1 method)

Set up the extra cache and extend prob_cache

Here, we construct an extra cache that holds a few extra intermediate variables, used in the routines above. But this cache also, crucially, holds the functions and operators of the constrained ODE function. We call the function ODEFunctionList to assemble these together.

The prob_cache function creates this ODE function, supplying the functions that we just defined. We also create a Laplacian operator with the heat diffusivity built into it. (This operator is singled out from the other terms in the heat conduction equation, because we account for it separately in the time marching using a matrix exponential.) We also create prototypes of the state and constraint force vectors. Here, the state is the grid temperature data and the constraint is the Lagrange multipliers on the boundary.

struct DirichletHeatConductionCache{DTT,FT} <: AbstractExtraILMCache
   dTb :: DTT
   f :: FT
end

function ImmersedLayers.prob_cache(prob::DirichletHeatConductionProblem,
                                   base_cache::BasicILMCache{N,scaling}) where {N,scaling}
    @unpack phys_params = prob
    @unpack gdata_cache, g = base_cache

    dTb = zeros_surface(base_cache)

    # Construct a Lapacian outfitted with the diffusivity
    κ = phys_params["diffusivity"]
    heat_L = Laplacian(base_cache,κ)

    # State (grid temperature data) and constraint (surface Lagrange multipliers)
    f = ODEFunctionList(state = zeros_grid(base_cache),
                        constraint = zeros_surface(base_cache),
                        ode_rhs=heatconduction_ode_rhs!,
                        lin_op=heat_L,
                        bc_rhs=heatconduction_bc_rhs!,
                        constraint_force = heatconduction_constraint_force!,
                        bc_op = heatconduction_bc_op!)

    DirichletHeatConductionCache(dTb,f)
end

Before we move on to solving the problem, we need to set up a function that will calculate the time step size. The time marching algorithm will call this function. Of course, this could just be used to specify a time step directly, e.g., by supplying it in phys_params. But it is better to use a stability condition (a Fourier condition) to determine it based on the other data.

function timestep_fourier(u,sys)
    @unpack phys_params = sys
    g = get_grid(sys)
    κ = phys_params["diffusivity"]
    Fo = phys_params["Fourier"]
    Δt = Fo*cellsize(g)^2/κ
    return Δt
end

timestep_fourier (generic function with 1 method)

Solve the problem

We will solve heat conduction inside a circular region with uniform temperature, with thermal diffusivity equal to 1.

Set up the grid

Δx = 0.01
Lx = 4.0
xlim = (-Lx/2,Lx/2)
ylim = (-Lx/2,Lx/2)
g = PhysicalGrid(xlim,ylim,Δx);

Set up the body shape.

Here, we will demonstrate the solution on a circular shape of radius 1.

Δs = 1.4*cellsize(g)
body = Circle(1.0,Δs);

Though the body is stationary, we still need to provide some minimal information about its placement and (lack of) motion. We do this with the help of the RigidBodyMotion structure. Technically, the placement of a body constitutes a basic joint with the inertial coordinate system. the MotionTransform below places the joint at the origin of the inertial coordinate system (the first argument), with no relative rotation (the second argument).

X = MotionTransform([0,0],0)
joint = Joint(X)
m = RigidBodyMotion(joint,body)

1 linked system(s) of bodies
   1 bodies
   1 joints

We don't have to do anything more here because the placement of the body is trivial. However, to demonstrate how we might do it in other problems,

x = zero_motion_state(body,m)
update_body!(body,x,m)

Circular body with 448 points and radius 1.0
   Current position: (0.0,0.0)
   Current angle (rad): 0.0

Specify the physical parameters, data, etc.

These can be changed later without having to regenerate the system.

Here, we create a dict with physical parameters to be passed in.

phys_params = Dict("diffusivity" => 1.0, "Fourier" => 1.0)

Dict{String, Float64} with 2 entries:
  "Fourier"     => 1.0
  "diffusivity" => 1.0

The temperature boundary functions on the exterior and interior are defined here and assembled into a dict. Note that these functions must have a slightly more complex signature than in time-invariant problems: for generality, they must accept the time argument and another argument accepting possible motions of the surfaces.

get_Tbplus(t,x,base_cache,phys_params,motions) = zeros_surface(base_cache)
get_Tbminus(t,x,base_cache,phys_params,motions) = ones_surface(base_cache)
bcdict = Dict("exterior" => get_Tbplus,"interior" => get_Tbminus)

Dict{String, Function} with 2 entries:
  "interior" => get_Tbminus
  "exterior" => get_Tbplus

Construct the problem, passing in the data and functions we've just created. We pass in the body's motion (however trivial) via the motions keyword.

prob = DirichletHeatConductionProblem(g,body,scaling=GridScaling,
                                             phys_params=phys_params,
                                             bc=bcdict,
                                             motions=m,
                                             timestep_func=timestep_fourier);

Construct the system

sys = construct_system(prob);

Solving the problem

In contrast to the previous (time-independent) example, we have not extended the solve function here to serve us in solving this problem. Instead, we rely on the tools in ConstrainedSystems.jl to advance the solution forward in time. This package builds from the OrdinaryDiffEq.jl package, and leverages most of the tools of that package.

Set an initial condition. Here, we just get an initial (zeroed) copy of the solution prototype that we have stored in the extra cache. We also get the time step size for our own inspection.

u0 = init_sol(sys)
Δt = timestep_fourier(u0,sys)

0.0001

It is instructive to note that u0 has two parts: a state and a constraint, each obtained respectively with a convenience function. The state in this case is the temperature; the constraint is the Lagrange multiplier.

state(u0)

Nodes{Primal, 405, 405, Float64, Matrix{Float64}} data
Printing in grid orientation (lower left is (1,1))
404×404 Matrix{Float64}:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  …  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  …  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 ⋮                        ⋮              ⋱                 ⋮              
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  …  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  …  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0

constraint(u0)

448 points of scalar-valued Float64 data
448-element Vector{Float64}:
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 ⋮
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0

Now, create the integrator, with a time interval of 0 to 1. We have not specified the algorithm here explicitly; it defaults to the LiskaIFHERK() time-marching algorithm, which is a second-order algorithm for constrained ODE systems that utilizes a matrix exponential (i.e., integrating factor) for the linear part of the problem (the Laplacian), and a half-explicit Runge-Kutta method for the constrained part. This method is most suitable for problems in which there is no dependence on the Lagrange multipliers in the constraint. (Such a problem is an index-2 differential-algebraic equation.) Another possible choice for this problem is the first-order Euler method, IFHEEuler(), which can be specified with the keyword alg=IFHEEuler().

For problems that do have a constraint that depends on the Lagrange multipliers such as the Neumann problem (in an upcoming example), then the default method switches to HETrapezoidalAB2(), which uses a half-explicit trapezoidal method for the constrained and linear parts, and 2nd-order Adams-Bashforth for the explicit part.

tspan = (0.0,1.0)
integrator = init(u0,tspan,sys)

t: 0.0
u: (Primal nodes in a (nx = 405, ny = 405) cell grid of type Float64 data
  Number of Primal nodes: (nx = 404, ny = 404), [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0])

Now advance the solution by 0.01 convective time units, by using the step! function, which steps through the solution.

step!(integrator,0.01)

Plot the solution

The integrator holds the most recent solution in the field u, which has the same type as our initial condition u0. Here, we plot the state of the system at the end of the interval.

plot(state(integrator.u),sys)

It would be nice to just define a function called temperature to get this more explicitly. We will do that here, and also apply a macro @snapshotoutput that automatically extends this function with some convenient interfaces. For example, if we simply pass in the integrator to temperature, it will pick off the u field for us.

temperature(T,σ,x,sys::ILMSystem,t) = T
@snapshotoutput temperature

Now we can write

plot(temperature(integrator),sys)

The solution history is in the field integrator.sol. The macro we called earlier enables temperature to work for this, as well, and we can obtain the temperature at any time in the interval of our solution. For example, to get the solution at time 0.51:

sol = integrator.sol
plot(temperature(sol,sys,0.0051),sys)

We can also get it for an array of times, e.g.,

temperature(sol,sys,0.0051:0.0001:0.0061);

This also gives us a spatially-interpolated version of the temperature field. To get this interpolated version at $t = 0.0051$, for example,

Tfcn = temperature_xy(sol,sys,0.0051)

404×404 extrapolate(scale(interpolate(OffsetArray(::Matrix{Float64}, 0:405, 0:405), BSpline(Cubic(Line(Interpolations.OnGrid())))), (-2.0100000000000002:0.01:2.02, -2.0100000000000002:0.01:2.02)), (Flat(), Flat())) with element type Float64:
 -1.77286e-17  -3.78596e-17  -6.35158e-17  …  -7.54779e-17  -6.05771e-17
 -3.40437e-17  -6.33286e-17  -1.04037e-16     -1.08414e-16  -8.0363e-17
 -4.7266e-17   -8.31912e-17  -1.37305e-16     -1.2848e-16   -9.10788e-17
 -5.68022e-17  -9.95877e-17  -1.63764e-16     -1.39737e-16  -9.66344e-17
 -6.29781e-17  -1.12677e-16  -1.83598e-16     -1.45121e-16  -9.7715e-17
 -6.79236e-17  -1.20193e-16  -1.96147e-16  …  -1.44664e-16  -9.43561e-17
 -7.19671e-17  -1.22719e-16  -2.00434e-16     -1.38065e-16  -8.67141e-17
 -7.40692e-17  -1.24228e-16  -2.00505e-16     -1.2541e-16   -7.48548e-17
 -7.64705e-17  -1.26257e-16  -1.99634e-16     -1.10926e-16  -6.39332e-17
 -8.14516e-17  -1.28433e-16  -2.00185e-16     -9.92013e-17  -5.90073e-17
  ⋮                                        ⋱                
  7.14172e-17   1.17779e-16   1.82534e-16  …  -7.30695e-18   1.58588e-17
  5.60095e-17   9.56898e-17   1.58566e-16     -2.88939e-17  -7.24523e-19
  3.81243e-17   7.34378e-17   1.34498e-16     -4.62025e-17  -1.43791e-17
  2.05563e-17   5.19896e-17   1.08213e-16     -5.42875e-17  -2.0337e-17
  3.99212e-18   3.05975e-17   8.03739e-17     -5.40531e-17  -2.02717e-17
 -1.03685e-17   9.48817e-18   5.1529e-17   …  -4.70992e-17  -1.67412e-17
 -2.1125e-17   -1.03597e-17   2.35658e-17     -3.55238e-17  -1.09431e-17
 -2.78918e-17  -2.57481e-17   2.94888e-18     -2.01651e-17  -1.65294e-18
 -2.75553e-17  -3.09548e-17  -4.88058e-18     -1.39145e-18   9.35426e-18

And then, to get the temperature at location $(x,y) = (-0.9,0)$, for example,

Tfcn(-0.9,0)

0.3584147231951137

We can also get a time array of these spatially-interpolated fields by supplying temperature_xy with an array of times.

Tfcn_array = temperature_xy(sol,sys,0.0051:0.0001:0.0061);

Then, each element in Tfcn_array corresponds to an entry in the given time array, and can be accessed at an $(x,y)$ point.

Tfcn_array[4](-0.9,0)

0.3728574046514253

Motions

It is straightforward to make bodies move in time-varying problems, using the RigidBodyMotion function. In the previous example we used this in a trivial fashion. Here, we will use it in a non-trivial example. The only caveat is that the time-stepping becomes slower in such problems, since the system operators must be regenerated at every time step.

The RigidBodyTools.jl package provides a versatile set of motions, both rigid-body and deforming, and associated tools.

Rigid body motion

For example, to simply make the body move at constant velocity 1 in the x direction.

Xp_to_jp = MotionTransform([0,0],0)
Xc_to_jc = MotionTransform([0,0],0)
dofs = [ConstantVelocityDOF(0),ConstantVelocityDOF(1),ConstantVelocityDOF(0)]
joint = Joint(FreeJoint2d,0,Xp_to_jp,1,Xc_to_jc,dofs)
m = RigidBodyMotion(joint,body)

1 linked system(s) of bodies
   1 bodies
   1 joints

The MotionTransform operator Xp_to_jp places the parent side of the joint at the origin of the inertial coordinate system (body 0) with no rotation, and Xc_to_jc places the child side of the joint at the origin of the body's system (body 1), also with no rotation. We are creating a free joint, with all degrees of freedom possibly in motion.

These three degrees of freedom are all assigned a prescribed behavior in the dofs vector. They are ordered as [rotation, x position, y position]. The rotational and y degrees of freedom are all set to zero velocity and the x degree is set to velocity of 1.

Surface deformation

Deformations to the body surface can be superposed on rigid body motions. Here's an example of a deforming motion on a stationary body

ufcn(x,y,t) = 0.25*x*y*cos(t)
vfcn(x,y,t) = 0.25*(x^2-y^2)*cos(t)
def = DeformationMotion(ufcn,vfcn)
X = MotionTransform([0,0],0)
joint = Joint(X)
m = RigidBodyMotion(joint,body,def)

1 linked system(s) of bodies
   1 bodies
   1 joints

Either of these would be provided in the motions keyword of the problem construction. Consult the documentation of RigidBodyTools.jl to learn more about these. However, for time-marching purposes, it is helpful to know that the maximum surface velocity is provided by the maxvelocity function:

x = zero_motion_state(body,m)
maxvelocity(body,x,m)

(0.24998770586539437, 225, 0.0, 1)

The first element of the output is the maximum velocity, the second is the index on which it occurs, the third is the time at which it occurs, and the fourth is the body on which it occurs.

Time-varying PDE functions

ODEFunctionList

zeros_sol(sys::ILMSystem)

init_sol(sys::ILMSystem)

init_sol(s::AbstractSpatialField,sys::ILMSystem)

ConstrainedSystems.init(u0,tspan,sys::ILMSystem,[alg=LiskaIFHERK()/HETrapezoidalAB2()])

timestep(u,sys::ILMSystem) -> Float64

evaluate_motion!(motion::ILMMotion,x::AbstractVector,t::Real)

evaluate_motion!(sys::ILMSystem,x::AbstractVector,t::Real)

surface_velocity_in_translating_frame!(vel,x,base_cache,motions,t)

This page was generated using Literate.jl.