Understand apply! of JuAFEM

Hey there,

I’m new to Julia and I wanted to transport my existing FEM Code written with FEniCS in C++ to JuAFEM code. First of all thanks to the developers of the package!

My problem is a SUPG stabilized advection diffusion system
\int c_{n+1} v \ dΩ + \int c_{n+1} \cdot (\delta_T \mathbf{u} \cdot \nabla v)+ \\ \Delta t\int D \cdot (∇v ⋅ ∇c_{n+1}) dΩ + \Delta t\int \mathbf{u} ⋅ ∇c_{n+1} v \ dΩ + \Delta t \int \mathbf{u} ⋅ ∇c_{n+1} ⋅ (δ_T ⋅ \mathbf{u} ⋅ ∇v) dΩ
= \int R \cdot v dΩ + \int R \cdot (δ_T \mathbf{u} ⋅ ∇v) dΩ + \int c_n\cdot v dΩ + \int c_n \cdot (\delta_T \mathbf{u} \cdot \nabla v)

sorry for the mess of equation. \mathbf{u} denotes the velocity, v is the test function and c the field of interest. R is in this particular case 0, so the first two terms of the right hand side can be dropped. Anyways, my main problem is that I don’t fully understand the differences of the apply! function. The matrix entries (of the stiffness matrix) are correct, I double checked them with hand derived values.
no_f
The gif aboves shows the evolution of c over time. At the left part of the domain a smooth boundary condition is applied which reaches its maximum at approximately t=350, c=1. The problem is the “wall” between the first and second finite element. Since I tried several possibilities, I can guarantee that I mess up something when I apply the boundary conditions.

using JuAFEM, SparseArrays, Plots
import CSV

theme(:gruvbox_dark)
#plotly()
pyplot()
#gr()
N = (100,)
h = 0.05 / N[1]
L = 5e-2
T = 1000
Δt = 1
w = 1.9128e-4 * (1/0.37)
D = 1e-9
left = zero(Vec{1})
right = L * ones(Vec{1})
grid = generate_grid(Line, N, left, right)
input_exp = []
output_exp = []
for row in CSV.File("new-data.csv"; delim=" ")
    push!(input_exp, row.I)
    push!(output_exp, row.O)
end

ip = Lagrange{1,RefCube,1}() #Interpolation
qr = QuadratureRule{1,RefCube}(3) #QuadratureRule
cv = CellScalarValues(qr, ip) #FEValues
dh = DofHandler(grid)
push!(dh, :c, 1) #add concentration field
close!(dh)

K = create_sparsity_pattern(dh)
M = create_sparsity_pattern(dh)
fill!(K.nzval, 1.0)
fill!(M.nzval, 1.0)

ch = ConstraintHandler(dh)

function inputExp(t)
    t_int = convert(Int, t)
    return input_exp[t_int]
end
dbc = Dirichlet(:c, getfaceset(grid, "left"), (x, t) -> inputExp(t))
add!(ch, dbc)
close!(ch)
update!(ch,1)

function doassemble(
    w,
    δT,
    cellvalues::CellScalarValues{dim},
    M::SparseMatrixCSC,
    dh::DofHandler,
) where {dim}
    n_basefuncs = getnbasefunctions(cellvalues)
    Me = zeros(n_basefuncs, n_basefuncs)
    assembler = start_assemble(M)
    @inbounds for cell in CellIterator(dh)
        fill!(Me, 0.0)
        reinit!(cellvalues, cell)
        for q_point = 1:getnquadpoints(cellvalues)
            dΩ = getdetJdV(cellvalues, q_point)
            for i = 1:n_basefuncs
                v = shape_value(cellvalues, q_point, i)
                ∇v = shape_gradient(cellvalues, q_point, i)
                for j = 1:n_basefuncs
                    c = shape_value(cellvalues, q_point, j)
                    Me[i, j] += c ⋅ v * dΩ  + c ⋅w ⋅ ∇v *δT * dΩ
                end
            end
        end
        assemble!(assembler, celldofs(cell), Me)
    end
    return M
end

function doassemble(
    D::Float64,
    w::Float64,
    δT::Float64,
    cellvalues::CellScalarValues{dim},
    K::SparseMatrixCSC,
    dh::DofHandler,
) where {dim}
    R = 0
    n_basefuncs = getnbasefunctions(cellvalues)
    Ke = zeros(n_basefuncs, n_basefuncs)
    fe = zeros(n_basefuncs)
    f = zeros(ndofs(dh))
    assembler = start_assemble(K, f)
    @inbounds for cell in CellIterator(dh)
        fill!(Ke, 0.0)
        fill!(fe, 0.0)
        reinit!(cellvalues, cell)
        for q_point = 1:getnquadpoints(cellvalues)
            dΩ = getdetJdV(cellvalues, q_point)
            for i = 1:n_basefuncs
                v = shape_value(cellvalues, q_point, i)
                ∇v = shape_gradient(cellvalues, q_point, i)
                fe[i] += R * v * dΩ + R * (δT * w ⋅ ∇v) * dΩ
                for j = 1:n_basefuncs
                    c = shape_value(cellvalues, q_point, j)
                    ∇c = shape_gradient(cellvalues, q_point, j)
                    Ke[i, j] +=
                        D * (∇v ⋅ ∇c) * dΩ +
                        w ⋅ ∇c ⋅ v * dΩ +
                        w ⋅ ∇c ⋅ (δT ⋅ w ⋅ ∇v) * dΩ
                end
            end
        end
        assemble!(assembler, celldofs(cell), fe, Ke)
    end
    return K, f
end


#δT = h/(2*abs(w)) # stabilization parameter
δT = 0. # stabilization parameter

K, f = doassemble(D, w, δT, cv, K, dh);
M = doassemble(w, δT, cv, M, dh);
A = M+(Δt*K)

c_0 = zeros(ndofs(dh))
c_n = copy(c_0)

store = []


for t in 1:Δt:T
    update!(ch, t) # load current dbc values from input_exp

    global b = f + M * c_n # get the discrete rhs
    apply!(A, b, ch) # apply time-dependent dbc
    #apply!(A, ch) #
    #apply!(b, ch) # this yields a completely different result
    global c = A \ b # solve the current time step

    push!(store, c) # store current solution
    global c_n = c # update time step
end

is the whole program, where the input is some experiment measurements and can be changed with anything you want. When reading this post, you should focus anyway on the time loop, i.e. this part:

K, f = doassemble(D, w, δT, cv, K, dh);
M = doassemble(w, δT, cv, M, dh);
A = M+(Δt*K)

c_0 = zeros(ndofs(dh))
c_n = copy(c_0)

store = []


for t in 1:Δt:T
    update!(ch, t) # load current dbc values from input_exp

    global b = f + M * c_n # get the discrete rhs
    apply!(A, b, ch) # apply time-dependent dbc
    #apply!(A, ch) #
    #apply!(b, ch) # this yields a completely different result
    global c = A \ b # solve the current time step

    push!(store, c) # store current solution
    global c_n = c # update time step
end

The stationary version works just fine. The problem is coming from the time stepping.

I want to know what the purpose of apply! is when either one matrix or vector is passed to the function. I guess, it’s something obvious I just don’t see. Maybe someone knows what goes wrong in the time loop or can explain the difference between

apply!(A, ch), apply!(b, ch) and apply(A, b, ch)

Thanks in advance!

1 Like

I think if you inspect the code, it’s pretty obvious what the apply! methods do. https://github.com/KristofferC/JuAFEM.jl/blob/master/src/Dofs/ConstraintHandler.jl#L307

So the apply! function “applies” the boundary condition values in some way. The vector-only version just writes the dof values to the corresponding elements in the vector. This is useful for inspection for example. The main method is the one that takes a matrix and a vector in. This method applies the boundary conditions by changing A and b to A' and b' such that the solution c = A' \ b' satisfies the boundary conditions. And for every unconstrained dof i, (A * c)[i] will be equal to b[i]. Assuming Dirichlet boundary conditions (bcs), this is done in a few steps:

  1. Add -v[i] * A[:,i] to b for each constrained dof i where v[i] is the value of the constrained dof.
  2. Find the mean diagonal value of the matrix m.
  3. Zero out the columns and rows of A corresponding to the bcs
  4. For each diagonal element A[i,i] corresponding to a bc, assign the value m to A[i,i] and the value m * v[i] to b[i] where v[i] is the value of the constrained dof i. Solving c = A \ b will now guarantee that c[i] == v[i] if i is a constrained dof hence the bc is satisfied. The c update part is ignored if only the matrix is passed to apply!.

Edit: edited to reflect @PetrKryslUCSD 's comment. Thanks!

3 Likes

The described procedure I think drops the interaction of the prescribed degrees of freedom with the free degrees of freedom. Normally this interaction results in a contribution to the right-hand side. Perhaps it is done this way in the code (I don’t have time to check), but it may be worthwhile to look into this.

1 Like

Yes my bad, I missed the step of adding the interaction. It is done correctly in the code.

Thank you guys!
I know its not that hard to understand the procedure from the source code, but I started two days ago with Julia. I appreciate the time you’ve taken to explain it. It seems that I have to revisit my code to get rid of the weird wall between the first two elements.
Thanks once again :slight_smile:

Can you post the CSV file somewhere? And maybe also the plotting code?

I suspect that the problem is that you modify A in each step. In Julia the convention for function names ending with ! is that they mutate their arguments. Maybe this diff is all you need to fix the problem:

--- a/main.jl
+++ b/main.jl
@@ -131,7 +131,7 @@ for t in 1:Δt:T
     update!(ch, t) # load current dbc values from input_exp
 
     global b = f + M * c_n # get the discrete rhs
-    apply!(A, b, ch) # apply time-dependent dbc
+    apply!(copy(A), b, ch) # apply time-dependent dbc
     #apply!(A, ch) #
     #apply!(b, ch) # this yields a completely different result
     global c = A \ b # solve the current time step
1 Like

Hey @fredrikekre
I’ll try that out! Here is the whole code with the plotting procedure:

using JuAFEM, SparseArrays, UnicodePlots, Plots
import Gadfly, CSV

theme(:gruvbox_dark)
plotly()
#pyplot()
#gr()
N = (100,)
h = 0.05 / N[1]
L = 5e-2
T = 1000
Δt = 1
w = 1.9128e-4 * (1/0.37)
D = 1e-9
left = zero(Vec{1})
right = L * ones(Vec{1})
grid = generate_grid(Line, N, left, right)
input_exp = []
output_exp = []
for row in CSV.File("new-data.csv"; delim=" ")
    push!(input_exp, row.I)
    push!(output_exp, row.O)
end

ip = Lagrange{1,RefCube,1}() #Interpolation
qr = QuadratureRule{1,RefCube}(3) #QuadratureRule
cv = CellScalarValues(qr, ip) #FEValues
dh = DofHandler(grid)
push!(dh, :c, 1) #add concentration field
close!(dh)

K = create_sparsity_pattern(dh)
M = create_sparsity_pattern(dh)
fill!(K.nzval, 1.0)
fill!(M.nzval, 1.0)
UnicodePlots.spy(K; height = 15)
UnicodePlots.spy(M; height=15)

ch = ConstraintHandler(dh)

function inputExp(x,t)
    t_int = convert(Int, t)
    return input_exp[t_int]
end

dbc = Dirichlet(:c, getfaceset(grid, "left"), (x, t) -> inputExp(x,t))
add!(ch, dbc)
close!(ch)
#update!(ch,1)

function doassemble(
    w,
    δT,
    cellvalues::CellScalarValues{dim},
    M::SparseMatrixCSC,
    dh::DofHandler,
) where {dim}
    n_basefuncs = getnbasefunctions(cellvalues)
    Me = zeros(n_basefuncs, n_basefuncs)
    assembler = start_assemble(M)
    @inbounds for cell in CellIterator(dh)
        fill!(Me, 0.0)
        reinit!(cellvalues, cell)
        for q_point = 1:getnquadpoints(cellvalues)
            dΩ = getdetJdV(cellvalues, q_point)
            for i = 1:n_basefuncs
                v = shape_value(cellvalues, q_point, i)
                ∇v = shape_gradient(cellvalues, q_point, i)
                for j = 1:n_basefuncs
                    c = shape_value(cellvalues, q_point, j)
                    Me[i, j] += c ⋅ v * dΩ  + c ⋅ w ⋅ ∇v ⋅ δT * dΩ
                end
            end
        end
        assemble!(assembler, celldofs(cell), Me)
    end
    return M
end

function doassemble(
    D::Float64,
    w::Float64,
    δT::Float64,
    cellvalues::CellScalarValues{dim},
    K::SparseMatrixCSC,
    dh::DofHandler,
) where {dim}
    R = 0
    n_basefuncs = getnbasefunctions(cellvalues)
    Ke = zeros(n_basefuncs, n_basefuncs)
    fe = zeros(n_basefuncs)
    f = zeros(ndofs(dh))
    assembler = start_assemble(K, f)
    @inbounds for cell in CellIterator(dh)
        fill!(Ke, 0.0)
        fill!(fe, 0.0)
        reinit!(cellvalues, cell)
        for q_point = 1:getnquadpoints(cellvalues)
            dΩ = getdetJdV(cellvalues, q_point)
            for i = 1:n_basefuncs
                v = shape_value(cellvalues, q_point, i)
                ∇v = shape_gradient(cellvalues, q_point, i)
                fe[i] += R * v * dΩ + R * (δT * w ⋅ ∇v) * dΩ
                for j = 1:n_basefuncs
                    c = shape_value(cellvalues, q_point, j)
                    ∇c = shape_gradient(cellvalues, q_point, j)
                    Ke[i, j] +=
                        D * (∇v ⋅ ∇c) * dΩ +
                        w ⋅ ∇c ⋅ v * dΩ +
                        w ⋅ ∇c ⋅ (δT ⋅ w ⋅ ∇v) * dΩ
                end
            end
        end
        assemble!(assembler, celldofs(cell), fe, Ke)
    end
    return K, f
end


δT = h/(2*abs(w)) # stabilization parameter
#δT = 0. # stabilization parameter

K, f = doassemble(D, w, δT, cv, K, dh);
M = doassemble(w, δT, cv, M, dh);
A = M+(Δt*K)

c_0 = zeros(ndofs(dh))
c_n = copy(c_0)

store = []
#b = f + M * c_n
#apply!(A, ch)

for t in 1:Δt:T
    update!(ch, t) # load current dbc values from input_exp

    global b = Δt*f + M * c_n # get the discrete rhs
    apply!(A, b, ch) # apply time-dependent dbc
    global c = A \ b # solve the current time step

    push!(store, c) # store current solution
    global c_n = copy(c) # update time step
end

function plotAnimation(storage::Array, gifname::String)
    t = 0
    anim = @animate for field in storage
        plot(field, ylim=(0,1), label="time=$t")
        t += 1
    end

    gif(anim, gifname, fps = 30)
end

plot(c, ylims=(0,1))


plotAnimation(store, "no_f.gif")

Unfortunately I’m not aware of how to post files in discourse (except pictures), so here is a dropbox upload of the csv: https://www.dropbox.com/s/cxalbxi6xmpjjdc/new-data.csv?dl=0
If you prefer any other hoster or a direct mail, just let me know. I think the bug is something obvious, due to my lacking julia experience.

BTW: since the developers are part of the conversation: Is it actually dumb to wrap my DirichletBC Array in a function? The !update can only handle Real types for time (which makes total sense), but I’m curious if you would’ve done it the same way :slight_smile:

This video has a bit more detailed description of the algorithm:

FWIW, I would highly recommend the whole lecture series by Wolfgang: https://www.math.colostate.edu/~bangerth/videos.html

3 Likes

Sorry, my patch above is not quite right, you of course need to use the copied A for the solution too:

--- a/main.jl
+++ b/main.jl
@@ -135,8 +135,9 @@ for t in 1:Δt:T
     update!(ch, t) # load current dbc values from input_exp
 
     global b = Δt*f + M * c_n # get the discrete rhs
-    apply!(A, b, ch) # apply time-dependent dbc
-    global c = A \ b # solve the current time step
+    copyA = copy(A)
+    apply!(copyA, b, ch) # apply time-dependent dbc
+    global c = copyA \ b # solve the current time step
 
     push!(store, c) # store current solution
     global c_n = copy(c) # update time step

which results in this:
no_f_fix

2 Likes

Damn I owe you guys a couple of beers.

Thanks for solving the problem @fredrikekre!