How to code initial condition in the "Time Dependant Problem" example ?

[Drachta]

Hello,

Sorry if this is a silly question (I'm currently still learning FEM techniques and Julia), but how do you add initial condition in the time dependant problem example ?

I guess it would be something like this, but I miss the details :

# upper half initial temperature is 100°C, bottom half is 0°C
@inbounds for cell in CellIterator(dh)
    if cell.pos[2] > 0
        cell.u=100.0
    end   
end

Thank you for your help !

[fredrikekre]

Since the initial conditions is usually on the whole domain (not just the boundary) you would have to perform an initial solve to compute u_0(x). If you have a known function u_0(x) and need to get that into the "FE-space" you can solve something like:

which, I think, is what the weak form of the heat equation reduced to at time t = 0. Solving this would be the correct way I think.

If you just want to compute the values for the dofs exactly, then there is no easy way to do this right now in Ferrite, but it would be good to have a way to interpolate like this, I guess.

The best way I can think of now is this "hack" of creating a constraint for all nodes and applying to your dof vector:

julia> using Ferrite

julia> grid = generate_grid(Line, (4,));

julia> dh = DofHandler(grid); push!(dh, :u, 1); close!(dh);

julia> ch = ConstraintHandler(dh);

julia> d = Dirichlet(:u, Set(1:getnnodes(grid)), (x, t) -> x[1]^2); # Initial condition that u(x, 0) = x^2

julia> add!(ch, d); close!(ch); update!(ch, 0);

julia> u0 = zeros(ndofs(dh))
5-element Vector{Float64}:
 0.0
 0.0
 0.0
 0.0
 0.0

julia> apply!(u0, ch)
5-element Vector{Float64}:
 1.0
 0.25
 0.0
 0.25
 1.0

Alternatively, you can compute the values in quadrature points and do an L2 interpolation, which should work nicely if the function is smooth at least. See this example. This would correspond to solving the initial equation (at least for the heat equation).

[Drachta]

Thank you for your suggestion !

Here is my code modification to the example:

uₙ = zeros(length(f));

# Initial condition such that T=100°C on upper half plate and T=0°C on lower half plate
initial_condition = ConstraintHandler(dh);

d = Dirichlet(:u, Set(1:getnnodes(grid)), (x, t) -> x[2] > 0 ? 100. : 0.); 
add!(initial_condition, d);

close!(initial_condition)
update!(initial_condition, 0.0);
apply!(uₙ, initial_condition);

apply!(A, ch);

pvd = paraview_collection("transient-heat.pvd");

vtk_grid("transient-heat-0", dh) do vtk
    vtk_point_data(vtk, dh, uₙ)
    vtk_save(vtk)
    pvd[0] = vtk
end

for t in Δt:Δt:t_max

I also commented the Dirichlet boundary condition at the beginning of the file, and replace T with t_max so I won't mistake it with the temperature.

The initial condition seems to be working, however the average temperature in the plate keeps getting higher as the time pass. I expected the temperature to stay at about 50°C once equilibrium is reached. What did I missed ?

[termi-official]

I guess you do not modify the RHS of your problem. Note that the line apply!(A, ch); just modified the system matrix A, but not your right hand side b. Does this example help?

[fredrikekre]

The initial condition seems to be working, however the average temperature in the plate keeps getting higher as the time pass. I expected the temperature to stay at about 50°C once equilibrium is reached. What did I missed ?

What boundary conditions do you use in the time interval?

[Drachta]

@termi-official I started from the example you propose.
@fredrikekre Basically, I just removed the Dirichlet boundary condition (by commenting lines 24 to 30 of the example, I guess I could just remove the apply!(A, ch); and apply_rhs!(rhsdata, b, ch) lines for the same effect), and replaced it with the initial condition (by inserting code between the line uₙ = zeros(length(f)); and the line for t in Δt:Δt:T as shown in my previous message).

What I try to achieve is to simulate the temperature evolution of a metal plate with adiabatic condition (so no heat leaves or enters the plate). I start with one half of the plate at 100°C, and the other half at 0°C. The temperature in the plate should rise or decrease to the average of 50°C depending on which node you are looking at.

However, what I get with my current code, is that the temperature of the whole plate keeps increasing, even after equilibrium is reached. The increase is linear with time (for delta_t = 0.1, the temperature increase is 0.1°C per step everywhere)

The way I understant FEM, is that I should not have any Dirichlet boundary condition, and instead have Neumann condition where \nabla u = 0 on the domain boundary: \delta\Omega.
If my understanding of FEM is correct, I don't need to add any code to achieve those Neumann condition ?

P.S. how do you use LateX in your message ?

[fredrikekre]

Is the system really solvable without any Dirichlet BC? I think you need it somewhere to avoid "rigid body movement", no?

P.S. how do you use LateX in your message ?

I used some online tex generator and copied the resulting png into the message :)

[fredrikekre]

Is the system really solvable without any Dirichlet BC?

I guess it is due to the mass matrix. With this diff:

diff --git a/docs/src/literate/transient_heat_equation.jl b/docs/src/literate/transient_heat_equation.jl
index eecb94c6..419b2503 100644
--- a/docs/src/literate/transient_heat_equation.jl
+++ b/docs/src/literate/transient_heat_equation.jl
@@ -94,13 +94,13 @@ T = 200
 ch = ConstraintHandler(dh);
 
 # Here, we define the boundary condition related to $\partial \Omega_1$.
-∂Ω₁ = union(getfaceset.((grid, ), ["left", "right"])...)
-dbc = Dirichlet(:u, ∂Ω₁, (x, t) -> 0)
-add!(ch, dbc);
+# ∂Ω₁ = union(getfaceset.((grid, ), ["left", "right"])...)
+# dbc = Dirichlet(:u, ∂Ω₁, (x, t) -> 0)
+# add!(ch, dbc);
 # While the next code block corresponds to the linearly increasing temperature description on $\partial \Omega_2$.
-∂Ω₂ = union(getfaceset.((grid, ), ["top", "bottom"])...)
-dbc = Dirichlet(:u, ∂Ω₂, (x, t) -> t*(max_temp/T))
-add!(ch, dbc)
+# ∂Ω₂ = union(getfaceset.((grid, ), ["top", "bottom"])...)
+# dbc = Dirichlet(:u, ∂Ω₂, (x, t) -> t*(max_temp/T))
+# add!(ch, dbc)
 close!(ch)
 update!(ch, 0.0);
 
@@ -127,7 +127,7 @@ function doassemble_K!(K::SparseMatrixCSC, f::Vector, cellvalues::CellScalarValu
             for i in 1:n_basefuncs
                 v  = shape_value(cellvalues, q_point, i)
                 ∇v = shape_gradient(cellvalues, q_point, i)
-                fe[i] += v * dΩ
+                # fe[i] += v * dΩ
                 for j in 1:n_basefuncs
                     ∇u = shape_gradient(cellvalues, q_point, j)
                     Ke[i, j] += (∇v ⋅ ∇u) * dΩ
@@ -182,6 +182,14 @@ A = (Δt .* K) + M;
 rhsdata = get_rhs_data(ch, A);
 # We set the initial time step, denoted by uₙ,  to $\mathbf{0}$.
 uₙ = zeros(length(f));
+
+initial_condition = ConstraintHandler(dh);
+d = Dirichlet(:u, Set(1:getnnodes(grid)), (x, t) -> x[2] > 0 ? 100. : 0.);
+add!(initial_condition, d);
+close!(initial_condition)
+update!(initial_condition, 0.0);
+apply!(uₙ, initial_condition);
+
 # Here, we apply **once** the boundary conditions to the system matrix `A`.
 apply!(A, ch);

it works for me and reaches a steady state of ~50.

The increase is linear with time (for delta_t = 0.1, the temperature increase is 0.1°C per step everywhere)

Perhaps you forgot to comment out the internal heat source added to f?

[Drachta]

Is the system really solvable without any Dirichlet BC? I think you need it somewhere to avoid "rigid body movement", no?

That would be true for a steady-state problem. For a transient problem, the rigid body movement should be handled by initial condition.

Perhaps you forgot to comment out the internal heat source added to f?

Yes, that was exactly it. Thank you very much !

[fredrikekre]

Yea, I realized that. Great that you got it working :)