Skip to content

Commit fed91ac

Browse files
authored
Merge pull request #81 from oriolcg/TransientDistributedCellField
Transient distributed cell field
2 parents 4036e66 + 5166991 commit fed91ac

17 files changed

+607
-16
lines changed

.gitignore

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -9,3 +9,4 @@
99
/tmp/
1010
*.vtu
1111
*.pvtu
12+
*.pvd

Project.toml

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -14,10 +14,10 @@ SparseMatricesCSR = "a0a7dd2c-ebf4-11e9-1f05-cf50bc540ca1"
1414
WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
1515

1616
[compat]
17-
Gridap = "0.17.8"
18-
PartitionedArrays = "0.2.10"
1917
FillArrays = "0.8.4, 0.9, 0.10, 0.11, 0.12"
18+
Gridap = "0.17.9"
2019
MPI = "0.16, 0.17, 0.18, 0.19"
20+
PartitionedArrays = "0.2.10"
2121
SparseMatricesCSR = "0.6.6"
2222
WriteVTK = "1.12.0"
2323
julia = "1.3"

src/GridapDistributed.jl

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -11,6 +11,8 @@ using Gridap.CellData
1111
using Gridap.Visualization
1212
using Gridap.FESpaces
1313
using Gridap.MultiField
14+
using Gridap.ODEs.TransientFETools
15+
using Gridap.ODEs.ODETools
1416

1517
using PartitionedArrays
1618
const PArrays = PartitionedArrays
@@ -23,6 +25,7 @@ import Gridap.TensorValues: inner, outer, double_contraction, symmetric_part
2325
import LinearAlgebra: det, tr, cross, dot,
2426
import Base: inv, abs, abs2, *, +, -, /, adjoint, transpose, real, imag, conj
2527
import Gridap.Fields: grad2curl
28+
import Gridap.ODEs.ODETools: ∂t, ∂tt
2629

2730
export FullyAssembledRows
2831
export SubAssembledRows
@@ -41,4 +44,10 @@ include("DivConformingFESpaces.jl")
4144

4245
include("MultiField.jl")
4346

47+
include("TransientDistributedCellField.jl")
48+
49+
include("TransientMultiFieldDistributedCellField.jl")
50+
51+
include("TransientFESpaces.jl")
52+
4453
end # module

src/TransientDistributedCellField.jl

Lines changed: 99 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,99 @@
1+
# Transient Distributed CellField
2+
abstract type TransientDistributedCellField <: DistributedCellDatum end
3+
4+
# Transient SingleField
5+
struct TransientSingleFieldDistributedCellField{A} <: TransientDistributedCellField
6+
cellfield::A
7+
derivatives::Tuple
8+
end
9+
10+
# Constructors
11+
function TransientFETools.TransientCellField(single_field::DistributedSingleFieldFEFunction,derivatives::Tuple)
12+
TransientSingleFieldDistributedCellField(single_field,derivatives)
13+
end
14+
15+
function TransientFETools.TransientCellField(single_field::DistributedCellField,derivatives::Tuple)
16+
TransientSingleFieldDistributedCellField(single_field,derivatives)
17+
end
18+
19+
# Time derivative
20+
function ∂t(f::TransientDistributedCellField)
21+
cellfield, derivatives = first_and_tail(f.derivatives)
22+
TransientCellField(cellfield,derivatives)
23+
end
24+
25+
∂tt(f::TransientDistributedCellField) = ∂t(∂t(f))
26+
27+
# Integration related
28+
function Fields.integrate(f::TransientDistributedCellField,b::DistributedMeasure)
29+
integrate(f.cellfield,b)
30+
end
31+
32+
# Differential Operations
33+
Fields.gradient(f::TransientDistributedCellField) = gradient(f.cellfield)
34+
Fields.∇∇(f::TransientDistributedCellField) = ∇∇(f.cellfield)
35+
36+
# Unary ops
37+
for op in (:symmetric_part,:inv,:det,:abs,:abs2,:+,:-,:tr,:transpose,:adjoint,:grad2curl,:real,:imag,:conj)
38+
@eval begin
39+
($op)(a::TransientDistributedCellField) = ($op)(a.cellfield)
40+
end
41+
end
42+
43+
# Binary ops
44+
for op in (:inner,:outer,:double_contraction,:+,:-,:*,:cross,:dot,:/)
45+
@eval begin
46+
($op)(a::TransientDistributedCellField,b::TransientDistributedCellField) = ($op)(a.cellfield,b.cellfield)
47+
($op)(a::TransientDistributedCellField,b::DistributedCellField) = ($op)(a.cellfield,b)
48+
($op)(a::DistributedCellField,b::TransientDistributedCellField) = ($op)(a,b.cellfield)
49+
($op)(a::TransientDistributedCellField,b::Number) = ($op)(a.cellfield,b)
50+
($op)(a::Number,b::TransientDistributedCellField) = ($op)(a,b.cellfield)
51+
($op)(a::TransientDistributedCellField,b::Function) = ($op)(a.cellfield,b)
52+
($op)(a::Function,b::TransientDistributedCellField) = ($op)(a,b.cellfield)
53+
end
54+
end
55+
56+
Base.broadcasted(f,a::TransientDistributedCellField,b::TransientDistributedCellField) = broadcasted(f,a.cellfield,b.cellfield)
57+
Base.broadcasted(f,a::TransientDistributedCellField,b::DistributedCellField) = broadcasted(f,a.cellfield,b)
58+
Base.broadcasted(f,a::DistributedCellField,b::TransientDistributedCellField) = broadcasted(f,a,b.cellfield)
59+
Base.broadcasted(f,a::Number,b::TransientDistributedCellField) = broadcasted(f,a,b.cellfield)
60+
Base.broadcasted(f,a::TransientDistributedCellField,b::Number) = broadcasted(f,a.cellfield,b)
61+
Base.broadcasted(f,a::Function,b::TransientDistributedCellField) = broadcasted(f,a,b.cellfield)
62+
Base.broadcasted(f,a::TransientDistributedCellField,b::Function) = broadcasted(f,a.cellfield,b)
63+
Base.broadcasted(a::typeof(*),b::typeof(∇),f::TransientDistributedCellField) = broadcasted(a,b,f.cellfield)
64+
Base.broadcasted(a::typeof(*),s::Fields.ShiftedNabla,f::TransientDistributedCellField) = broadcasted(a,s,f.cellfield)
65+
66+
dot(::typeof(∇),f::TransientDistributedCellField) = dot(∇,f.cellfield)
67+
outer(::typeof(∇),f::TransientDistributedCellField) = outer(∇,f.cellfield)
68+
outer(f::TransientDistributedCellField,::typeof(∇)) = outer(f.cellfield,∇)
69+
cross(::typeof(∇),f::TransientDistributedCellField) = cross(∇,f.cellfield)
70+
71+
# Skeleton related
72+
function Base.getproperty(f::TransientDistributedCellField, sym::Symbol)
73+
if sym in (:⁺,:plus,:⁻, :minus)
74+
derivatives = ()
75+
cellfield = DistributedCellField(f.cellfield,sym)
76+
for iderivative in f.derivatives
77+
derivatives = (derivatives...,DistributedCellField(iderivative))
78+
end
79+
return TransientSingleFieldCellField(cellfield,derivatives)
80+
else
81+
return getfield(f, sym)
82+
end
83+
end
84+
85+
Base.propertynames(x::TransientDistributedCellField, private::Bool=false) = propertynames(x.cellfield, private)
86+
87+
for op in (:outer,:*,:dot)
88+
@eval begin
89+
($op)(a::TransientDistributedCellField,b::SkeletonPair{<:DistributedCellField}) = ($op)(a.cellfield,b)
90+
($op)(a::SkeletonPair{<:DistributedCellField},b::TransientDistributedCellField) = ($op)(a,b.cellfield)
91+
end
92+
end
93+
94+
Arrays.evaluate!(cache,k::Operation,a::TransientDistributedCellField,b::SkeletonPair{<:DistributedCellField}) = evaluate!(cache,k,a.cellfield,b)
95+
96+
Arrays.evaluate!(cache,k::Operation,a::SkeletonPair{<:DistributedCellField},b::TransientDistributedCellField) = evaluate!(cache,k,a,b.cellfield)
97+
98+
CellData.jump(a::TransientDistributedCellField) = jump(a.cellfield)
99+
CellData.mean(a::TransientDistributedCellField) = mean(a.cellfield)

src/TransientFESpaces.jl

Lines changed: 66 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,66 @@
1+
# Functions for transient FE spaces
2+
3+
Fields.evaluate(U::DistributedSingleFieldFESpace,t::Nothing) = U
4+
5+
(U::DistributedSingleFieldFESpace)(t) = U
6+
7+
∂t(U::DistributedSingleFieldFESpace) = HomogeneousTrialFESpace(U)
8+
∂t(U::DistributedMultiFieldFESpace) = MultiFieldFESpace(∂t.(U.field_fe_space))
9+
∂tt(U::DistributedSingleFieldFESpace) = HomogeneousTrialFESpace(U)
10+
∂tt(U::DistributedMultiFieldFESpace) = MultiFieldFESpace(∂tt.(U.field_fe_spaces))
11+
12+
function TransientMultiFieldFESpace(spaces::Vector{<:DistributedSingleFieldFESpace})
13+
MultiFieldFESpace(spaces)
14+
end
15+
16+
# Functions for transient FE Functions
17+
18+
function ODETools.allocate_jacobian(
19+
op::TransientFETools.TransientFEOperatorFromWeakForm,
20+
duh::DistributedCellField,
21+
cache)
22+
_matdata_jacobians = TransientFETools.fill_initial_jacobians(op,duh)
23+
matdata = _vcat_distributed_matdata(_matdata_jacobians)
24+
allocate_matrix(op.assem_t,matdata)
25+
end
26+
27+
function ODETools.allocate_jacobian(
28+
op::TransientFETools.TransientFEOperatorFromWeakForm,
29+
duh::DistributedMultiFieldFEFunction,
30+
cache)
31+
_matdata_jacobians = TransientFETools.fill_initial_jacobians(op,duh)
32+
matdata = _vcat_distributed_matdata(_matdata_jacobians)
33+
allocate_matrix(op.assem_t,matdata)
34+
end
35+
36+
function ODETools.jacobians!(
37+
A::AbstractMatrix,
38+
op::TransientFETools.TransientFEOperatorFromWeakForm,
39+
t::Real,
40+
xh::TransientDistributedCellField,
41+
γ::Tuple{Vararg{Real}},
42+
cache)
43+
_matdata_jacobians = TransientFETools.fill_jacobians(op,t,xh,γ)
44+
matdata = _vcat_distributed_matdata(_matdata_jacobians)
45+
assemble_matrix_add!(A,op.assem_t, matdata)
46+
A
47+
end
48+
49+
function _vcat_distributed_matdata(_matdata)
50+
term_to_cellmat = map_parts(a->a[1],local_views(_matdata[1]))
51+
term_to_cellidsrows = map_parts(a->a[2],local_views(_matdata[1]))
52+
term_to_cellidscols = map_parts(a->a[3],local_views(_matdata[1]))
53+
for j in 2:length(_matdata)
54+
term_to_cellmat_j = map_parts(a->a[1],local_views(_matdata[j]))
55+
term_to_cellidsrows_j = map_parts(a->a[2],local_views(_matdata[j]))
56+
term_to_cellidscols_j = map_parts(a->a[3],local_views(_matdata[j]))
57+
term_to_cellmat = map_parts((a,b)->vcat(a,b),local_views(term_to_cellmat),local_views(term_to_cellmat_j))
58+
term_to_cellidsrows = map_parts((a,b)->vcat(a,b),local_views(term_to_cellidsrows),local_views(term_to_cellidsrows_j))
59+
term_to_cellidscols = map_parts((a,b)->vcat(a,b),local_views(term_to_cellidscols),local_views(term_to_cellidscols_j))
60+
end
61+
map_parts( (a,b,c) -> (a,b,c),
62+
local_views(term_to_cellmat),
63+
local_views(term_to_cellidsrows),
64+
local_views(term_to_cellidscols)
65+
)
66+
end
Lines changed: 61 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,61 @@
1+
# Transient MultiField
2+
struct TransientMultiFieldDistributedCellField{A} <: TransientDistributedCellField
3+
cellfield::A
4+
derivatives::Tuple
5+
transient_single_fields::Vector{<:TransientDistributedCellField} # used to iterate
6+
end
7+
8+
# Constructors
9+
function TransientFETools.TransientCellField(multi_field::DistributedMultiFieldFEFunction,derivatives::Tuple)
10+
transient_single_fields = _to_transient_single_distributed_fields(multi_field,derivatives)
11+
TransientMultiFieldDistributedCellField(multi_field,derivatives,transient_single_fields)
12+
end
13+
14+
# Get single index
15+
function Base.getindex(f::TransientMultiFieldDistributedCellField,ifield::Integer)
16+
single_field = f.cellfield[ifield]
17+
single_derivatives = ()
18+
for ifield_derivatives in f.derivatives
19+
single_derivatives = (single_derivatives...,getindex(ifield_derivatives,ifield))
20+
end
21+
TransientSingleFieldDistributedCellField(single_field,single_derivatives)
22+
end
23+
24+
# Get multiple indices
25+
function Base.getindex(f::TransientMultiFieldDistributedCellField,indices::Vector{<:Int})
26+
cellfield = DistributedMultiFieldCellField(f.cellfield[indices],DomainStyle(f.cellfield))
27+
derivatives = ()
28+
for derivative in f.derivatives
29+
derivatives = (derivatives...,DistributedMultiFieldCellField(derivative[indices],DomainStyle(derivative)))
30+
end
31+
transient_single_fields = _to_transient_single_distributed_fields(cellfield,derivatives)
32+
TransientMultiFieldDistributedCellField(cellfield,derivatives,transient_single_fields)
33+
end
34+
35+
function _to_transient_single_distributed_fields(multi_field,derivatives)
36+
transient_single_fields = TransientDistributedCellField[]
37+
for ifield in 1:num_fields(multi_field)
38+
single_field = multi_field[ifield]
39+
single_derivatives = ()
40+
for ifield_derivatives in derivatives
41+
single_derivatives = (single_derivatives...,getindex(ifield_derivatives,ifield))
42+
end
43+
transient_single_field = TransientSingleFieldDistributedCellField(single_field,single_derivatives)
44+
push!(transient_single_fields,transient_single_field)
45+
end
46+
transient_single_fields
47+
end
48+
49+
# Iterate functions
50+
Base.iterate(f::TransientMultiFieldDistributedCellField) = iterate(f.transient_single_fields)
51+
Base.iterate(f::TransientMultiFieldDistributedCellField,state) = iterate(f.transient_single_fields,state)
52+
53+
# Time derivative
54+
function ∂t(f::TransientMultiFieldDistributedCellField)
55+
cellfield, derivatives = first_and_tail(f.derivatives)
56+
transient_single_field_derivatives = TransientDistributedCellField[]
57+
for transient_single_field in f.transient_single_fields
58+
push!(transient_single_field_derivatives,∂t(transient_single_field))
59+
end
60+
TransientMultiFieldDistributedCellField(cellfield,derivatives,transient_single_field_derivatives)
61+
end

test/HeatEquationTests.jl

Lines changed: 77 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,77 @@
1+
module HeatEquationTests
2+
3+
using Gridap
4+
using Gridap.ODEs
5+
using GridapDistributed
6+
using PartitionedArrays
7+
using Test
8+
9+
function main(parts)
10+
θ = 0.2
11+
12+
u(x,t) = (1.0-x[1])*x[1]*(1.0-x[2])*x[2]*t
13+
u(t::Real) = x -> u(x,t)
14+
f(t) = x -> ∂t(u)(x,t)-Δ(u(t))(x)
15+
16+
domain = (0,1,0,1)
17+
partition = (4,4)
18+
model = CartesianDiscreteModel(parts,domain,partition)
19+
20+
order = 2
21+
22+
reffe = ReferenceFE(lagrangian,Float64,order)
23+
V0 = FESpace(
24+
model,
25+
reffe,
26+
conformity=:H1,
27+
dirichlet_tags="boundary"
28+
)
29+
U = TransientTrialFESpace(V0,u)
30+
31+
Ω = Triangulation(model)
32+
degree = 2*order
33+
= Measure(Ω,degree)
34+
35+
#
36+
m(u,v) = (u*v)dΩ
37+
a(u,v) = ((v)(u))dΩ
38+
b(t,v) = (v*f(t))dΩ
39+
40+
res(t,u,v) = a(u,v) + m(∂t(u),v) - b(t,v)
41+
jac(t,u,du,v) = a(du,v)
42+
jac_t(t,u,dut,v) = m(dut,v)
43+
44+
op = TransientFEOperator(res,jac,jac_t,U,V0)
45+
op_constant = TransientConstantMatrixFEOperator(m,a,b,U,V0)
46+
47+
t0 = 0.0
48+
tF = 1.0
49+
dt = 0.1
50+
51+
U0 = U(0.0)
52+
uh0 = interpolate_everywhere(u(0.0),U0)
53+
54+
ls = LUSolver()
55+
ode_solver = ThetaMethod(ls,dt,θ)
56+
57+
sol_t = solve(ode_solver,op,uh0,t0,tF)
58+
sol_t_const = solve(ode_solver,op_constant,uh0,t0,tF)
59+
60+
l2(w) = w*w
61+
62+
tol = 1.0e-6
63+
64+
for (uh_tn, tn) in sol_t
65+
e = u(tn) - uh_tn
66+
el2 = sqrt(sum( (l2(e))dΩ ))
67+
@test el2 < tol
68+
end
69+
70+
for (uh_tn, tn) in sol_t_const
71+
e = u(tn) - uh_tn
72+
el2 = sqrt(sum( (l2(e))dΩ ))
73+
@test el2 < tol
74+
end
75+
end
76+
77+
end #module

0 commit comments

Comments
 (0)