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.
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!
