Skip to content

Commit

Permalink
Fixing GalerkinToolkitExamples
Browse files Browse the repository at this point in the history
  • Loading branch information
@fverdugo
fverdugo committed Dec 23, 2024
1 parent d9f8acb commit c6d5e45
Show file tree
Hide file tree
Showing 6 changed files with 83 additions and 42 deletions.
10 changes: 5 additions & 5 deletions GalerkinToolkitExamples/Project.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,8 @@ authors = ["Francesc Verdugo <[email protected]>"]
version = "0.1.0"

[deps]
AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
GalerkinToolkit = "5e3ba9c4-dd81-444d-b69a-0e7bd7bf60a4"
IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
Expand All @@ -20,14 +22,12 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"

[compat]
PartitionedSolvers = "0.3"

[extras]
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[targets]
test = ["Test"]

[compat]
PartitionedSolvers = "0.3"

31 changes: 15 additions & 16 deletions GalerkinToolkitExamples/src/poisson.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,8 @@ using SparseArrays
using PartitionedArrays
using TimerOutputs
using WriteVTK
using AbstractTrees


function main(params_in)
params_default = default_params()
Expand All @@ -26,18 +28,15 @@ function main(params_in)
end
end

gradient(u) = x->ForwardDiff.gradient(u,x)
jacobian(u) = x->ForwardDiff.jacobian(u,x)
laplacian(u,x) = tr(ForwardDiff.jacobian(y->ForwardDiff.gradient(u,y),x))
laplacian(u) = x->laplacian(u,x)
# TODO hide this
laplacian(u::GT.AbstractQuantity,x::GT.AbstractQuantity) = GT.call(laplacian,u,x)

Δ(u) = GT.call(laplacian,u)
∇(u) = GT.call(gradient,u)
Δ(u,x) = Δ(u)(x)
∇(u,x) = ∇(u)(x)
const ∇ = ForwardDiff.gradient
const Δ = laplacian

mean(u,x) = 0.5*(u(x)[+]+u(x)[-])
jump(u,n,x) = u(x)[+]*n[+](x) + u(x)[-]*n[-](x)
mean(u) = 0.5*(u[1]+u[2])
jump(u,n) = u[2]*n[2] + u[1]*n[1]

function main_automatic(params)
timer = params[:timer]
Expand All @@ -54,7 +53,7 @@ function main_automatic(params)

u = GT.analytical_field(params[:u],Ω)
f(x) = -Δ(u,x)
n_Γn = GT.unit_normal(Γn,Ω)
n_Γn = GT.unit_normal(mesh,D-1)
g(x) = n_Γn(x)⋅∇(u,x)

interpolation_degree = params[:interpolation_degree]
Expand All @@ -68,7 +67,7 @@ function main_automatic(params)
GT.label_interior_faces!(mesh;physical_name="__INTERIOR_FACES__")
Λ = GT.skeleton(mesh;physical_names=["__INTERIOR_FACES__"])
dΛ = GT.measure(Λ,integration_degree)
n_Λ = GT.unit_normal(Λ,Ω)
n_Λ = GT.unit_normal(mesh,D-1)
h_Λ = GT.face_diameter_field(Λ)
else
conformity = :default
Expand All @@ -79,7 +78,7 @@ function main_automatic(params)
uhd = GT.dirichlet_field(Float64,V)
GT.interpolate_dirichlet!(u,uhd)
else
n_Γd = GT.unit_normal(Γd,Ω)
n_Γd = GT.unit_normal(mesh,D-1)
h_Γd = GT.face_diameter_field(Γd)
dΓd = GT.measure(Γd,integration_degree)
V = GT.lagrange_space(Ω,interpolation_degree;conformity)
Expand All @@ -98,7 +97,7 @@ function main_automatic(params)
end
if params[:discretization_method] === :interior_penalty
r += ∫( x->
(γ/h_Λ(x))*jump(v,n_Λ,x)⋅jump(u,n_Λ,x)-jump(v,n_Λ,x)⋅mean(∇(u),x)-mean(∇(v),x)⋅jump(u,n_Λ,x), dΛ)
(γ/h_Λ(x))*jump(v(x),n_Λ(x))⋅jump(u(x),n_Λ(x))-jump(v(x),n_Λ(x))⋅mean(∇(u,x))-mean(∇(v,x))⋅jump(u(x),n_Λ(x)), dΛ)
end
r
end
Expand Down Expand Up @@ -202,7 +201,7 @@ function main_hand_written(params)
dΩ = GT.measure(Ω,integration_degree)

u = params[:u]
f(x) = -laplacian(u,x)
f(x) = -Δ(u,x)

# Set interpolation space -> Lagrange polynomial basis function, the element type and geometry
V = GT.lagrange_space(Ω,interpolation_degree;dirichlet_boundary=Γd)
Expand Down Expand Up @@ -564,8 +563,8 @@ function default_params()
params[:dirichlet_tags] = ["boundary"]
params[:neumann_tags] = String[]
params[:dirichlet_method] = :strong
params[:integration_degree] = 1
params[:interpolation_degree] = 2*params[:integration_degree]
params[:interpolation_degree] = 1
params[:integration_degree] = 2*params[:interpolation_degree]
params[:solver] = PS.LinearAlgebra_lu
params[:example_path] = joinpath(mkpath(joinpath(@__DIR__,"..","output")),"example_001")
params[:export_vtu] = true
Expand Down
4 changes: 2 additions & 2 deletions GalerkinToolkitExamples/test/poisson_tests.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@ for discretization_method in (:continuous_galerkin,)
params[:dirichlet_tags] = ["1-face-1","1-face-2","1-face-3","1-face-4"]
params[:discretization_method] = discretization_method
params[:dirichlet_method] = dirichlet_method
params[:integration_degree] = 2
params[:interpolation_degree] = 2
results = Poisson.main(params)
@test results[:error_l2_norm] < tol
@test results[:error_h1_norm] < tol
Expand All @@ -35,7 +35,7 @@ for discretization_method in (:interior_penalty,:continuous_galerkin)
params[:neumann_tags] = ["1-face-2"]
params[:discretization_method] = discretization_method
params[:dirichlet_method] = dirichlet_method
params[:integration_degree] = 2
params[:interpolation_degree] = 2
results = Poisson.main(params)
@test results[:error_l2_norm] < tol
@test results[:error_h1_norm] < tol
Expand Down
60 changes: 45 additions & 15 deletions src/assembly.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -192,11 +192,17 @@ function assemble_vector_count(integral,state)
for field in 1:nfields
Dface_to_sDface = field_to_Dface_to_sDface[field]
D = field_to_D[field]
if D < d
continue
end
face_to_Dfaces = face_incidence(topo,d,D)
sDface_to_dofs = field_to_sDface_to_dofs[field]
Dfaces = face_to_Dfaces[face]
for Dface in Dfaces
sDface = Dface_to_sDface[Dface]
if sDface == 0
continue
end
dofs = sDface_to_dofs[sDface]
for dof in dofs
counter = vector_strategy.count(counter,setup,dof,field)
Expand Down Expand Up @@ -230,7 +236,7 @@ function assemble_vector_fill!(integral,state)
nfields = GT.num_fields(space)
field_to_domain = map(field->domain(space,field),1:nfields)
field_to_sDface_to_dofs = map(field->face_dofs(space,field),1:nfields)
field_to_Dface_to_sDface = map(faces,field_to_domain)
field_to_Dface_to_sDface = map(inverse_faces,field_to_domain)
field_to_D = map(num_dims,field_to_domain)
counter = vector_strategy.counter(setup)
T = vector_strategy.scalar_type(setup)
Expand Down Expand Up @@ -270,13 +276,19 @@ function assemble_vector_fill!(integral,state)
$(s_qty[3])
Dface_to_sDface = args.field_to_Dface_to_sDface[$field]
D = args.field_to_D[$field]
if D < args.d
continue
end
face_to_Dfaces = face_incidence(args.topo,args.d,D)
sDface_to_dofs = args.field_to_sDface_to_dofs[$field]
Dfaces = face_to_Dfaces[$face]
for ($Dface_around,Dface) in enumerate(Dfaces)
$(s_qty[4])
fill!(b,zero(eltype(b)))
sDface = Dface_to_sDface[Dface]
if sDface == 0
continue
end
dofs = sDface_to_dofs[sDface]
for $point in 1:npoints
$(s_qty[5])
Expand Down Expand Up @@ -386,20 +398,32 @@ function assemble_matrix_count(integral,state)
for field_test in 1:axis_to_nfields[test]
test_Dface_to_sDface = axis_to_field_to_Dface_to_sDface[test][field_test]
test_D = axis_to_field_to_D[test][field_test]
if test_D < d
continue
end
test_face_to_Dfaces = face_incidence(topo,d,test_D)
test_sDface_to_dofs = axis_to_field_to_sDface_to_dofs[test][field_test]
test_Dfaces = test_face_to_Dfaces[face]
for field_trial in 1:axis_to_nfields[trial]
trial_Dface_to_sDface = axis_to_field_to_Dface_to_sDface[trial][field_trial]
trial_D = axis_to_field_to_D[trial][field_trial]
if trial_D < d
continue
end
trial_face_to_Dfaces = face_incidence(topo,d,trial_D)
trial_sDface_to_dofs = axis_to_field_to_sDface_to_dofs[trial][field_trial]
trial_Dfaces = trial_face_to_Dfaces[face]
for test_Dface in test_Dfaces
test_sDface = test_Dface_to_sDface[test_Dface]
if test_sDface == 0
continue
end
test_dofs = test_sDface_to_dofs[test_sDface]
for trial_Dface in trial_Dfaces
trial_sDface = trial_Dface_to_sDface[trial_Dface]
if trial_sDface == 0
continue
end
trial_dofs = trial_sDface_to_dofs[trial_sDface]
for test_dof in test_dofs
for trial_dof in trial_dofs
Expand Down Expand Up @@ -453,10 +477,10 @@ function assemble_matrix_fill!(integral,state)
axes = (test, trial)

idof_test, idof_trial = map(axis -> dof_index(index,axis), axes) # TODO: test, trial
test_Dface, trial_Dface = map(axis -> face_around_index(index,axis), axes)
test_Dface_around, trial_Dface_around = map(axis -> face_around_index(index,axis), axes)
field_test, field_trial = map(axis -> field_index(index,axis), axes)

s_qty = GT.topological_sort(expr_qty,(face, field_test, field_trial, test_Dface, trial_Dface, point, idof_test, idof_trial)) # 8 deps
s_qty = GT.topological_sort(expr_qty,(face, field_test, field_trial, test_Dface_around, trial_Dface_around, point, idof_test, idof_trial)) # 8 deps
s_npoints = GT.topological_sort(expr_npoints,(face,))

expr = quote
Expand All @@ -466,51 +490,57 @@ function assemble_matrix_fill!(integral,state)
$(s_npoints[1])
b = args.b
nsfaces = length(args.sface_to_face)

for sface in 1:nsfaces
$face = args.sface_to_face[sface]
$(s_qty[2])
npoints = $(s_npoints[2])

for $field_test in 1:args.axis_to_nfields[$test]
$(s_qty[3])
test_Dface_to_sDface = args.axis_to_field_to_Dface_to_sDface[$test][$field_test]
test_D = args.axis_to_field_to_D[$test][$field_test]
if test_D < args.d
continue
end
test_face_to_Dfaces = face_incidence(args.topo,args.d,test_D)
test_sDface_to_dofs = args.axis_to_field_to_sDface_to_dofs[$test][$field_test]
test_Dfaces = test_face_to_Dfaces[$face]
for $field_trial in 1:args.axis_to_nfields[$trial]
$(s_qty[4])
trial_Dface_to_sDface = args.axis_to_field_to_Dface_to_sDface[$trial][$field_trial]
trial_D = args.axis_to_field_to_D[$trial][$field_trial]
if trial_D < args.d
continue
end
trial_face_to_Dfaces = face_incidence(args.topo,args.d,trial_D)
trial_sDface_to_dofs = args.axis_to_field_to_sDface_to_dofs[$trial][$field_trial]
trial_Dfaces = trial_face_to_Dfaces[$face]

for $test_Dface in test_Dfaces
for ($test_Dface_around,test_Dface) in enumerate(test_Dfaces)
$(s_qty[5])
test_sDface = test_Dface_to_sDface[$test_Dface]
test_sDface = test_Dface_to_sDface[test_Dface]
if test_sDface == 0
continue
end
test_dofs = test_sDface_to_dofs[test_sDface]
for $trial_Dface in trial_Dfaces
for ($trial_Dface_around,trial_Dface) in enumerate(trial_Dfaces)
$(s_qty[6])
trial_sDface = trial_Dface_to_sDface[$trial_Dface]
trial_sDface = trial_Dface_to_sDface[trial_Dface]
if trial_sDface == 0
continue
end
trial_dofs = trial_sDface_to_dofs[trial_sDface]
fill!(b,zero(eltype(b)))


for $point in 1:npoints
$(s_qty[7])
for ($idof_test, dof_test) in enumerate(test_dofs)
$(s_qty[8])
for ($idof_trial, dof_trial) in enumerate(trial_dofs)
b[$idof_trial, $idof_test] += $(s_qty[9])
b[$idof_test, $idof_trial] += $(s_qty[9])
end
end
end

for ($idof_test, dof_test) in enumerate(test_dofs)
for ($idof_trial, dof_trial) in enumerate(trial_dofs)
counter = args.matrix_strategy.set!(args.alloc,counter,args.setup,b[$idof_trial, $idof_test],dof_test,dof_trial,$field_test,$field_trial)
counter = args.matrix_strategy.set!(args.alloc,counter,args.setup,b[$idof_test, $idof_trial],dof_test,dof_trial,$field_test,$field_trial)
end
end
end
Expand Down
2 changes: 1 addition & 1 deletion src/interpolation.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1004,7 +1004,7 @@ function form_argument(a::AbstractSpace,axis,field=1)
refid_to_reffes = reference_fes(a)
refid_to_funs = map(GT.shape_functions,refid_to_reffes)
domain = GT.domain(a)
qty = form_argument(axis,field,refid_to_funs,domain;reference=true,face_reference_id=face_to_refid)
qty = form_argument(axis,field,refid_to_funs,domain;reference=true)# TODO,face_reference_id=face_to_refid)
if is_reference_domain(GT.domain(a))
return qty
end
Expand Down
18 changes: 15 additions & 3 deletions src/symbolics.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -332,6 +332,15 @@ function binary_call_term(f,a,b::BoundaryTerm)
boundary_term(b.dim,b.dim_around,t,b.face_around)
end

function binary_call_term(f,a::BoundaryTerm,b)
t = binary_call_term(f,a.term,b)
boundary_term(a.dim,a.dim_around,t,a.face_around)
end

function binary_call_term(f,a::BoundaryTerm,b::BoundaryTerm)
BinaryCallTerm(f,a,b)
end

skeleton_term(args...) = SkeletonTerm(args...)

function skeleton_term(dim,dim_around,term::BoundaryTerm)
Expand Down Expand Up @@ -825,7 +834,10 @@ end
free_dims(a::PhysicalMapTerm) = [a.dim]


prototype(a::PhysicalMapTerm) = x-> zero(SVector{val_parameter(a.dim),Float64})
function prototype(a::PhysicalMapTerm)
D = num_ambient_dims(mesh(index(a)))
x-> zero(SVector{D,Float64})
end


function Base.:(==)(a::PhysicalMapTerm,b::PhysicalMapTerm)
Expand Down Expand Up @@ -976,7 +988,7 @@ function binary_call_term_physical_maps(::typeof(call),a,b,phi1::PhysicalMapTerm
f(phi(x))
end

function binary_call_term_physical_maps(::typeof(ForwardDiff.gradient),a,b,phi1::PhysicalMapTerm,phi2::PhysicalMapTerm)
function binary_call_term_physical_maps(d::typeof(ForwardDiff.gradient),a,b,phi1::PhysicalMapTerm,phi2::PhysicalMapTerm)
@assert phi1.dim != phi2.dim
phi = reference_map_term(phi2.dim,phi1.dim,index(a))
f = a.arg1
Expand Down Expand Up @@ -1237,7 +1249,7 @@ struct UnitNormalTerm{A} <: AbstractTerm
end

free_dims(a::UnitNormalTerm) = free_dims(a.outwards_normals)
prototype(a::UnitNormalTerm) = prototype(a.outwards_normals)
prototype(a::UnitNormalTerm) = x -> prototype(a.outwards_normals)
index(a::UnitNormalTerm) = index(a.outwards_normals)

function expression(c::BinaryCallTerm{typeof(call),<:UnitNormalTerm,<:Any})
Expand Down

0 comments on commit c6d5e45

Please sign in to comment.