Linear Time Dependent Problem

This example is the same example as FESolvers.jl's transient heat flow which was taken from Ferrite.jl's transient heat flow. Please see the theoretical derivations in those examples, with the specific formulation used here in the former.

Commented Program

Now we solve the problem by using FerriteProblems. The full program, without comments, can be found in the next section.

First we load required packages

using Ferrite, FerriteProblems, FerriteAssembly, FESolvers
import FerriteProblems as FP
import FerriteAssembly as FA
import FerriteAssembly.ExampleElements: TransientFourier

Physics

The transient heat element, TransientFourier, is available from FerriteAssembly.ExampleElements, noting that in the original example, the heat capacity, c=1, is implicitly assumed.

material() = TransientFourier(#=k=#1.0-3, #=c=#1.0)

material (generic function with 1 method)

The material doesn't include the heat source: We later add this with FerriteAssembly's LoadHandler, which works similar to the ConstraintHandler.

Problem setup

We start by a function that will create the problem definition

function create_definition()
    grid = generate_grid(Quadrilateral, (100, 100));
    ip = Lagrange{RefQuadrilateral,1}()
    cellvalues = CellValues(
        QuadratureRule{RefQuadrilateral}(2), ip);

    dh = DofHandler(grid); add!(dh, :u, ip); close!(dh)

    # Boundary conditions
    # Zero pressure on ∂Ω₁ and linear ramp followed by constant pressure on ∂Ω₂
    max_temp = 100; t_rise = 100
    ch = ConstraintHandler(dh);
    ∂Ω₁ = union(getfaceset.((grid,), ["left", "right"])...)
    add!(ch, Dirichlet(:u, ∂Ω₁, (x, t) -> 0));
    ∂Ω₂ = union(getfaceset.((grid,), ["top", "bottom"])...)
    add!(ch, Dirichlet(:u, ∂Ω₂, (x, t) -> max_temp * clamp(t / t_rise, 0, 1)))
    close!(ch)

    # Body load: constant heat source, f=5.0e-1
    lh = LoadHandler(dh)
    add!(lh, BodyLoad(:u, 2, Returns(5.0e-1)))

    # Create and return the `FEDefinition`
    domainspec = DomainSpec(dh, material(), cellvalues)
    return FEDefinition(domainspec; ch, lh)
end;

Postprocessing

After defining all the physics and problem setup, we must decide what data to save. In this example, we use the vtk-file exports as in the original example. To this end, we define the custom postprocessing struct

struct TH_PostProcessing{PVD}
    pvd::PVD
end
TH_PostProcessing() = TH_PostProcessing(paraview_collection("transient-heat.pvd"));

And the postprocessing function that is called after each time step

function FESolvers.postprocess!(post::TH_PostProcessing, p, solver)
    step = FESolvers.get_step(solver)
    if step < 5 || mod(step, 20) == 0
        @info "postprocessing step $step"
        dh = FP.get_dofhandler(p)
        vtk_grid("transient-heat-$step", dh) do vtk
            vtk_point_data(vtk, dh, FP.getunknowns(p))
            vtk_save(vtk)
            post.pvd[step] = vtk
        end
    end
end;

At the end of the simulation, we want to finish all IO operations. We can then define the function close_postprocessing which will be called even in the case that an error is thrown during the simulation

function FP.close_postprocessing(post::TH_PostProcessing, p)
    vtk_save(post.pvd)
end;

And now we create the problem type, and define the QuasiStaticSolver with the LinearProblemSolver as well as fixed time steps

def = create_definition()
post = TH_PostProcessing()
problem = FerriteProblem(def, post)
solver = QuasiStaticSolver(;nlsolver=LinearProblemSolver(), timestepper=FixedTimeStepper(collect(0.0:1.0:200)));

Finally, we can solve the problem

solve_problem!(problem, solver);

[ Info: postprocessing step 1
[ Info: postprocessing step 2
[ Info: postprocessing step 3
[ Info: postprocessing step 4
[ Info: postprocessing step 20
[ Info: postprocessing step 40
[ Info: postprocessing step 60
[ Info: postprocessing step 80
[ Info: postprocessing step 100
[ Info: postprocessing step 120
[ Info: postprocessing step 140
[ Info: postprocessing step 160
[ Info: postprocessing step 180
[ Info: postprocessing step 200

Plain program

Here follows a version of the program without any comments. The file is also available here: transient_heat.jl.

using Ferrite, FerriteProblems, FerriteAssembly, FESolvers
import FerriteProblems as FP
import FerriteAssembly as FA
import FerriteAssembly.ExampleElements: TransientFourier

material() = TransientFourier(#=k=#1.0-3, #=c=#1.0)

function create_definition()
    grid = generate_grid(Quadrilateral, (100, 100));
    ip = Lagrange{RefQuadrilateral,1}()
    cellvalues = CellValues(
        QuadratureRule{RefQuadrilateral}(2), ip);

    dh = DofHandler(grid); add!(dh, :u, ip); close!(dh)

    # Boundary conditions
    # Zero pressure on ∂Ω₁ and linear ramp followed by constant pressure on ∂Ω₂
    max_temp = 100; t_rise = 100
    ch = ConstraintHandler(dh);
    ∂Ω₁ = union(getfaceset.((grid,), ["left", "right"])...)
    add!(ch, Dirichlet(:u, ∂Ω₁, (x, t) -> 0));
    ∂Ω₂ = union(getfaceset.((grid,), ["top", "bottom"])...)
    add!(ch, Dirichlet(:u, ∂Ω₂, (x, t) -> max_temp * clamp(t / t_rise, 0, 1)))
    close!(ch)

    # Body load: constant heat source, f=5.0e-1
    lh = LoadHandler(dh)
    add!(lh, BodyLoad(:u, 2, Returns(5.0e-1)))

    # Create and return the `FEDefinition`
    domainspec = DomainSpec(dh, material(), cellvalues)
    return FEDefinition(domainspec; ch, lh)
end;

struct TH_PostProcessing{PVD}
    pvd::PVD
end
TH_PostProcessing() = TH_PostProcessing(paraview_collection("transient-heat.pvd"));

function FESolvers.postprocess!(post::TH_PostProcessing, p, solver)
    step = FESolvers.get_step(solver)
    if step < 5 || mod(step, 20) == 0
        @info "postprocessing step $step"
        dh = FP.get_dofhandler(p)
        vtk_grid("transient-heat-$step", dh) do vtk
            vtk_point_data(vtk, dh, FP.getunknowns(p))
            vtk_save(vtk)
            post.pvd[step] = vtk
        end
    end
end;

function FP.close_postprocessing(post::TH_PostProcessing, p)
    vtk_save(post.pvd)
end;

def = create_definition()
post = TH_PostProcessing()
problem = FerriteProblem(def, post)
solver = QuasiStaticSolver(;nlsolver=LinearProblemSolver(), timestepper=FixedTimeStepper(collect(0.0:1.0:200)));

solve_problem!(problem, solver);

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl

This page was generated using Literate.jl.