Skip to content

Commit

Permalink
Checkpoint, for the records. Code not functional.
Browse files Browse the repository at this point in the history
  • Loading branch information
@amartinhuertas
amartinhuertas committed Jul 16, 2024
1 parent 96452da commit 4c24c21
Show file tree
Hide file tree
Showing 2 changed files with 245 additions and 18 deletions.
110 changes: 110 additions & 0 deletions canvas.jl
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,110 @@
using Gridap
using GridapGeosciences

p_ex(x) = -x[1]x[2]x[3]
u_ex(x) = VectorValue(
x[2]x[3] - 3x[1]x[1]x[2]x[3] / (x[1]x[1] + x[2]x[2] + x[3]x[3]),
x[1]x[3] - 3x[1]x[2]x[2]x[3] / (x[1]x[1] + x[2]x[2] + x[3]x[3]),
x[1]x[2] - 3x[1]x[2]x[3]x[3] / (x[1]x[1] + x[2]x[2] + x[3]x[3])
)
# function uθϕ_ex(θϕ)
# u_ex(θϕ2xyz(θϕ))
# end
# function pθϕ_ex(θϕ)
# p_ex(θϕ2xyz(θϕ))
# end
# function Gridap.Fields.gradient(::typeof(p_ex))
# gradient_unit_sphere(pθϕ_ex)∘xyz2θϕ
# end
# function Gridap.Fields.divergence(::typeof(u_ex))
# divergence_unit_sphere(uθϕ_ex)∘xyz2θϕ
# end
f_ex(x) = -12x[1]x[2]x[3] # divergence(u_ex)(x)
g_ex(x) = u_ex(x)+gradient(p_ex)(x)

# using Distributions
# θϕ=Point(rand(Uniform(0,2*pi)),rand(Uniform(-pi/2,pi/2)))
# θϕ2xyz(θϕ)

function assemble_darcy_problem(model, order)
rt_reffe = ReferenceFE(raviart_thomas, Float64, order)
lg_reffe = ReferenceFE(lagrangian, Float64, order)
V = FESpace(model, rt_reffe, conformity=:Hdiv)
U = TrialFESpace(V)

Q = FESpace(model, lg_reffe; conformity=:L2)
P = TrialFESpace(Q)

# X = MultiFieldFESpace([U, P])
# Y = MultiFieldFESpace([V, Q])

Ω = Triangulation(model)
degree = 10 # 2*order
dΩ = Measure(Ω, degree)
dω = Measure(Ω,degree,ReferenceDomain())

a11(u,v)=∫(v⋅u)dΩ
B11=assemble_matrix(a11,U,V)

a12(p,v)=∫(∇⋅v*p)dΩ
B12=assemble_matrix(a12,P,V)

a21(u,q)=∫(∇⋅u*q)dΩ
B21=assemble_matrix(a12,P,V)

a22(p,q)=∫(p*q)dΩ
B22=assemble_matrix(a22,P,Q)

# a((u, p), (v, q)) = ∫(v ⋅ u + p*q)dΩ + ∫(∇⋅(u)*q -∇⋅(v)*p)dΩ
# l((v, q)) = ∫(q*f_ex + q*p_ex)dΩ
# AffineFEOperator(a, l, X, Y)
B11, B12, B21, B22
end

function solve_darcy_problem(op)
solve(op)
end

function compute_darcy_errors(model, order, xh)
uh,ph=xh
eph = ph-p_ex
euh = uh-u_ex
degree=2*order
Ω = Triangulation(model)
dΩ = Measure(Ω, degree)
err_p = sum(∫(eph*eph)dΩ)
err_u_l2 = sum(∫(euh⋅euh)dΩ)
err_u_div = sum(∫(euh⋅euh + (∇⋅(euh))*(∇⋅(euh)))dΩ)
err_p, err_u_l2, err_u_div
surface=sum(∫(1)dΩ)
surface, surface, surface
end

order=0
n=1
model = CubedSphereDiscreteModel(n ; radius=1.0)
B211,B212,B221,B222=assemble_darcy_problem(model, order);

for geom_order in 1:8
order=0
model = CubedSphereDiscreteModel(n, geom_order; radius=1.0)
B111,B112,B121,B122=assemble_darcy_problem(model, order)
println(norm(B111-B211)/norm(B211), " ",
norm(B112-B212)/norm(B212), " ",
norm(B121-B221)/norm(B221), " ",
norm(B122-B222)/norm(B222))
# println(B1)
end

#for n in [4,8,16,25,40]
for n in [1,]
model = CubedSphereDiscreteModel(n; radius=1.0)
op=assemble_darcy_problem(model, order)
println(Array(op.op.matrix))
# xh=solve_darcy_problem(op)
# err_p, err_u_l2, err_u_div = compute_darcy_errors(model, order, xh)
# println(err_p, " ", err_u_l2, " ", err_u_div)
end



153 changes: 135 additions & 18 deletions test/mpi/DarcyCubedSphereTests.jl
View file
Original file line number Diff line number Diff line change
@@ -1,14 +1,11 @@
module DarcyCubedSphereTestsMPI
# module DarcyCubedSphereTestsMPI
using PartitionedArrays
using Test
using FillArrays
using Gridap
using GridapPETSc
using GridapGeosciences

include("../DarcyCubedSphereTests.jl")
include("../ConvergenceAnalysisTools.jl")

function petsc_gamg_options()
"""
-ksp_type gmres -ksp_rtol 1.0e-06 -ksp_atol 0.0
Expand All @@ -24,23 +21,143 @@ module DarcyCubedSphereTestsMPI
-pc_type lu -pc_factor_mat_solver_type mumps
"""
end

p_ex(x) = -x[1]x[2]x[3]
u_ex(x) = VectorValue(
x[2]x[3] - 3x[1]x[1]x[2]x[3] / (x[1]x[1] + x[2]x[2] + x[3]x[3]),
x[1]x[3] - 3x[1]x[2]x[2]x[3] / (x[1]x[1] + x[2]x[2] + x[3]x[3]),
x[1]x[2] - 3x[1]x[2]x[3]x[3] / (x[1]x[1] + x[2]x[2] + x[3]x[3])
)
f_ex(x) = -12x[1]x[2]x[3]
function uθϕ_ex(θϕ)
u_ex(θϕ2xyz(θϕ))
end
function pθϕ_ex(θϕ)
p_ex(θϕ2xyz(θϕ))
end
function Gridap.Fields.gradient(::typeof(p_ex))
gradient_unit_sphere(pθϕ_ex)∘xyz2θϕ
end
function Gridap.Fields.divergence(::typeof(u_ex))
divergence_unit_sphere(uθϕ_ex)∘xyz2θϕ
end

# using Distributions
# θϕ=Point(rand(Uniform(0,2*pi)),rand(Uniform(-pi/2,pi/2)))
# θϕ2xyz(θϕ)

function assemble_darcy_problem(model, order)
rt_reffe = ReferenceFE(raviart_thomas, Float64, order)
lg_reffe = ReferenceFE(lagrangian, Float64, order)
V = FESpace(model, rt_reffe, conformity=:Hdiv)
U = TrialFESpace(V)

Q = FESpace(model, lg_reffe; conformity=:L2)
P = TrialFESpace(Q)

X = MultiFieldFESpace([U, P])
Y = MultiFieldFESpace([V, Q])

Ω = Triangulation(model)
degree = 10 # 2*order
dΩ = Measure(Ω, degree)

a11(u,v)=∫(v⋅u)dΩ
B11=assemble_matrix(a11,U,V)

a12(p,v)=∫(∇⋅v*p)dΩ
B12=assemble_matrix(a12,P,V)

a21(u,q)=∫(∇⋅u*q)dΩ
B21=assemble_matrix(a12,P,V)

a22(p,q)=∫(p*q)dΩ
B22=assemble_matrix(a22,P,Q)

B11, B12, B21, B22

# a((u, p), (v, q)) = ∫(v ⋅ u - (∇ ⋅ v)p + (∇ ⋅ u)q + p*q)dΩ
# l((v, q)) = ∫(q*f + q*p_ex)dΩ
# AffineFEOperator(a, l, X, Y)
end

function solve_darcy_problem(op)
solve(op)
end

function compute_darcy_errors(model, order, xh, op)
uh,ph=xh
eph = ph-p_ex
euh = uh-u_ex
degree=10 #2*order
Ω = Triangulation(model)
dΩ = Measure(Ω, degree)
err_p = sum(∫(eph*eph)dΩ)
err_u_l2 = sum(∫(euh⋅euh)dΩ)
err_u_div = sum(∫(euh⋅euh + (∇⋅(euh))*(∇⋅(euh)))dΩ)
err_p, err_u_l2, err_u_div
# U,P=op.trial
# uh=interpolate(u_ex,U)
# ph=interpolate(p_ex,P)
# eph = ph-p_ex
# euh = uh-u_ex
# err_p = sum(∫(eph*eph)dΩ)
# err_u_l2 = sum(∫(euh⋅euh)dΩ)
# err_u_div = sum(∫(euh⋅euh + (∇⋅(euh))*(∇⋅(euh)))dΩ)
# err_p, err_u_l2, err_u_div
end

function main(distribute,parts)
ranks = distribute(LinearIndices((prod(parts),)))
GridapPETSc.with(args=split(petsc_gamg_options())) do
num_refs=[2,3,4,5]
hs=[2.0/2^n for n in num_refs]
model_args_series=zip(Fill(ranks,length(num_refs)),num_refs)
t = PartitionedArrays.PTimer(ranks,verbose=true)
PartitionedArrays.tic!(t)
hs1,k1errors,s1=convergence_study(solve_darcy,hs,model_args_series,1,4,PETScLinearSolver())
PartitionedArrays.toc!(t,"DarcyCubedSphereConvergenceStudy")
display(t)
println(hs1)
println(k1errors)
println(s1)
GridapPETSc.with(args=split(petsc_mumps_options())) do
order=0
num_uniform_refinements=0
model=CubedSphereDiscreteModel(
ranks,
num_uniform_refinements;
radius=1.0,
adaptive=false,
order=-1)
B211,B212,B221,B222=assemble_darcy_problem(model, order);

for geom_order in 1:8
model = CubedSphereDiscreteModel(ranks,
num_uniform_refinements;
radius=1.0,
order=geom_order,
adaptive=true)
B111,B112,B121,B122=assemble_darcy_problem(model, order);
println(norm(B111.matrix_partition.item_ref[]-B211.matrix_partition.item_ref[])/norm(B211.matrix_partition.item_ref[]), " ",
norm(B112.matrix_partition.item_ref[]-B212.matrix_partition.item_ref[])/norm(B212.matrix_partition.item_ref[]), " ",
norm(B121.matrix_partition.item_ref[]-B221.matrix_partition.item_ref[])/norm(B221.matrix_partition.item_ref[]), " ",
norm(B122.matrix_partition.item_ref[]-B222.matrix_partition.item_ref[])/norm(B222.matrix_partition.item_ref[]))
# println(B1)
end

# for num_uniform_refinements=0:0
# order=0
# model=CubedSphereDiscreteModel(
# ranks,
# num_uniform_refinements;
# radius=1.0,
# adaptive=true,
# order=2)
# op=assemble_darcy_problem(model, order)
# xh=solve_darcy_problem(op)
# println(Array(op.op.matrix.matrix_partition.item_ref[]))
# println(compute_darcy_errors(model, order, xh, op))
# # println(norm(op.op.matrix.item_ref[]))
# # for num_uniform_refinements=1:6
# # order=0
# # op=assemble_darcy_problem(model, order)
# # xh=solve_darcy_problem(op)
# # println(compute_darcy_errors(model, order, xh))
# # end
# end
end
end
with_mpi() do distribute
main(distribute,4)
end
main(distribute,1)
end;

end #module

0 comments on commit 4c24c21

Please sign in to comment.