Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Enhance example002 #54

Merged
merged 4 commits into from
Feb 7, 2024
Merged

Enhance example002 #54

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
fc9e75b
cleanup in example002
fverdugo Feb 7, 2024
7c32c7e
Running example002 in parallel with a petsc solver
fverdugo Feb 7, 2024
9063a32
Adding missing file
fverdugo Feb 7, 2024
a6c26df
Fix test
fverdugo Feb 7, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions extensions/GalerkinToolkitExamples/Project.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,9 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
GalerkinToolkit = "5e3ba9c4-dd81-444d-b69a-0e7bd7bf60a4"
IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
PartitionedArrays = "5a9dfac6-5c52-46f7-8278-5e2210713be9"
PetscCall = "1194c000-87c4-4102-b4a0-a6217ec4849e"
Preconditioners = "af69fa37-3177-5a40-98ee-561f696e4fcd"
SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
Expand Down
211 changes: 2 additions & 209 deletions extensions/GalerkinToolkitExamples/src/example001.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -222,14 +222,14 @@ end

function assemble_system(state)
timer = state.timer
@timeit timer "assemble_sytem_coo" args = assemble_sytem_coo(state)
@timeit timer "assemble_system_coo" args = assemble_system_coo(state)
I,J,V,b = args
n_dofs = state.dofs.n_dofs
@timeit timer "sparse" A = sparse(I,J,V,n_dofs,n_dofs)
A,b
end

function assemble_sytem_coo(state)
function assemble_system_coo(state)
timer = state.timer
@timeit timer "assemble_sybmolic" args = assemble_sybmolic(state)
@timeit timer "assemble_numeric_cells!" assemble_numeric_cells!(args...,state)
Expand Down Expand Up @@ -548,211 +548,4 @@ function export_results(uh,params,state)
nothing
end

### From here, it is old stuff to to be deleted

function assemble_system_loop(state)

cell_to_nodes = state.cell_isomap.face_to_nodes
cell_to_rid = state.cell_isomap.face_to_rid
free_and_dirichlet_nodes = state.dirichlet_bcs.free_and_dirichlet_nodes
node_to_x = state.cell_isomap.node_to_coords
∇s = state.cell_isomap.rid_to_shape_grads
s = state.cell_isomap.rid_to_shape_vals
u_dirichlet = state.dirichlet_bcs.u_dirichlet
f = state.user_funs.f
w = state.cell_integration.rid_to_weights
d = state.cell_integration.d

# Count coo entries
n_coo = 0
ncells = length(cell_to_nodes)
node_to_free_node = gk.permutation(free_and_dirichlet_nodes)
nfree = length(first(free_and_dirichlet_nodes))
for cell in 1:ncells
nodes = cell_to_nodes[cell]
for nj in nodes
if node_to_free_node[nj] > nfree
continue
end
for ni in nodes
if node_to_free_node[ni] > nfree
continue
end
n_coo += 1
end
end
end

# Allocate coo values
I_coo = zeros(Int32,n_coo)
J_coo = zeros(Int32,n_coo)
V_coo = zeros(Float64,n_coo)
b = zeros(nfree)

# Allocate auxiliary buffers
∇x = map(i->similar(i,size(i,2)),∇s)
Aes = map(i->zeros(size(i,2),size(i,2)),∇s)
fes = map(i->zeros(size(i,2)),∇s)
ues = map(i->zeros(size(i,2)),∇s)

Tx = SVector{d,Float64}
TJ = typeof(zero(Tx)*zero(Tx)')

# Fill coo values
i_coo = 0
for cell in 1:ncells
rid = cell_to_rid[cell]
nodes = cell_to_nodes[cell]
Ae = Aes[rid]
ue = ues[rid]
fe = fes[rid]
nl = length(nodes)
∇se = ∇s[rid]
∇xe = ∇x[rid]
se = s[rid]
we = w[rid]
nq = length(we)
fill!(Ae,zero(eltype(Ae)))
fill!(fe,zero(eltype(Ae)))
fill!(ue,zero(eltype(Ae)))
for k in 1:nl
nk = nodes[k]
gk = node_to_free_node[nk]
if gk <= nfree
continue
end
ue[k] = u_dirichlet[gk-nfree]
end
for iq in 1:nq
Jt = zero(TJ)
xint = zero(Tx)
for k in 1:nl
x = node_to_x[nodes[k]]
∇sqx = ∇se[iq,k]
sqx = se[iq,k]
Jt += ∇sqx*x'
xint += sqx*x
end
detJt = det(Jt)
invJt = inv(Jt)
dV = abs(detJt)*we[iq]
for k in 1:nl
∇xe[k] = invJt*∇se[iq,k]
end
for j in 1:nl
for i in 1:nl
Ae[i,j] += ∇xe[i]⋅∇xe[j]*dV
end
end
fx = f(xint)
for k in 1:nl
fe[k] += fx*se[iq,k]*dV
end
end
for i in 1:nl
for j in 1:nl
fe[i] -= Ae[i,j]*ue[j]
end
end
for i in 1:nl
ni = nodes[i]
gi = node_to_free_node[ni]
if gi > nfree
continue
end
b[gi] += fe[i]
end
for j in 1:nl
nj = nodes[j]
gj = node_to_free_node[nj]
if gj > nfree
continue
end
for i in 1:nl
ni = nodes[i]
gi = node_to_free_node[ni]
if gi > nfree
continue
end
i_coo += 1
I_coo[i_coo] = gi
J_coo[i_coo] = gj
V_coo[i_coo] = Ae[i,j]
end
end
end
add_neumann_bcs!(b,state)
I_coo,J_coo,V_coo,b
end

function add_neumann_bcs!(b,state)

neum_face_to_face = state.neumann_bcs.neum_face_to_face
if length(neum_face_to_face) == 0
return nothing
end
face_to_rid = state.face_isomap.face_to_rid
face_to_nodes = state.face_isomap.face_to_nodes
∇s_f = state.face_isomap.rid_to_shape_grads
s_f = state.face_isomap.rid_to_shape_vals
node_to_x = state.face_isomap.node_to_coords
free_and_dirichlet_nodes = state.dirichlet_bcs.free_and_dirichlet_nodes
g = state.user_funs.g
w_f = state.face_integration.rid_to_weights
d = state.cell_integration.d

fes = map(i->zeros(size(i,2)),s_f)
Tx = SVector{d,Float64}
Ts = typeof(zero(SVector{d-1,Float64})*zero(Tx)')
node_to_free_node = gk.permutation(free_and_dirichlet_nodes)
nfree = length(first(free_and_dirichlet_nodes))

for face in neum_face_to_face
rid = face_to_rid[face]
nodes = face_to_nodes[face]
∇se = ∇s_f[rid]
se = s_f[rid]
we = w_f[rid]
nq = length(we)
fe = fes[rid]
fill!(fe,zero(eltype(fe)))
nl = length(nodes)
for iq in 1:nq
Jt = zero(Ts)
xint = zero(Tx)
for k in 1:nl
x = node_to_x[nodes[k]]
∇sqx = ∇se[iq,k]
sqx = se[iq,k]
Jt += ∇sqx*x'
xint += sqx*x
end
J = transpose(Jt)
dS = sqrt(det(Jt*J))*we[iq]
gx = g(xint)
for k in 1:nl
fe[k] += gx*se[iq,k]*dS
end
end
for i in 1:nl
ni = nodes[i]
gi = node_to_free_node[ni]
if gi > nfree
continue
end
b[gi] += fe[i]
end
end

nothing
end

function assemble_system_old(state)
I,J,V,b = assemble_system_loop(state)
nfree = length(b)
A = sparse(I,J,V,nfree,nfree)
A,b
end


end # module
Loading
Loading