Skip to content

Commit

Permalink
Activated tests
Browse files Browse the repository at this point in the history
  • Loading branch information
JordiManyer committed Dec 11, 2024
1 parent 6062c1d commit 4563876
Show file tree
Hide file tree
Showing 5 changed files with 91 additions and 65 deletions.
3 changes: 1 addition & 2 deletions src/BlockSolvers/StaggeredFEOperators.jl
Original file line number Diff line number Diff line change
Expand Up @@ -210,8 +210,7 @@ function get_operator!(op_k::AffineFEOperator, op::StaggeredAffineFEOperator{NB}
A, b = get_matrix(op_k), get_vector(op_k)
a(uk,vk) = op.biforms[k](xhs,uk,vk)
l(vk) = op.liforms[k](xhs,vk)
uhd = zero(op.trials[k])
assemble_matrix_and_vector!(a,l,A,b,op.assems[k],op.trials[k],op.tests[k],uhd)
assemble_matrix_and_vector!(a,l,A,b,op.assems[k],op.trials[k],op.tests[k],zero(op.trials[k]))
return op_k
end

Expand Down
136 changes: 73 additions & 63 deletions test/BlockSolvers/StaggeredFEOperatorsTests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -7,13 +7,13 @@ using BlockArrays
using GridapSolvers
using GridapSolvers.BlockSolvers, GridapSolvers.LinearSolvers, GridapSolvers.NonlinearSolvers

function test_solution(xh,sol,X,dΩ)
function test_solution(xh,sol,X,dΩ,verbose)
N = length(sol)
xh_exact = interpolate(sol,X)
for k in 1:N
eh_k = xh[k] - xh_exact[k]
e_k = sqrt(sum((eh_k * eh_k)dΩ))
println("Error in field $k: $e_k")
verbose && println("Error in field $k: $e_k")
@test e_k < 1e-10
end
end
Expand Down Expand Up @@ -43,7 +43,7 @@ function driver_affine(model,verbose)
l3((u1,(u2,u3)),v4) = (sol[4] * (u1 + u2) * v4)dΩ

# Define solver: Each block will be solved with a CG solver
lsolver = CGSolver(JacobiLinearSolver();rtol=1.e-10,verbose=verbose)
lsolver = CGSolver(JacobiLinearSolver();rtol=1.e-12,verbose=verbose)
solver = StaggeredFESolver(fill(lsolver,3))

# Create operator from full spaces
Expand All @@ -52,7 +52,7 @@ function driver_affine(model,verbose)
Y = MultiFieldFESpace([V,V,V,V];style=mfs)
op = StaggeredAffineFEOperator([a1,a2,a3],[l1,l2,l3],X,Y)
xh = solve(solver,op)
test_solution(xh,sol,X,dΩ)
test_solution(xh,sol,X,dΩ,verbose)

# Create operator from components
UB1, VB1 = U1, V
Expand All @@ -64,72 +64,82 @@ function driver_affine(model,verbose)
xh = zero(X)
xh, cache = solve!(xh,solver,op);
xh, cache = solve!(xh,solver,op,cache);
test_solution(xh,sol,X,dΩ)
test_solution(xh,sol,X,dΩ,verbose)

return true
end

function test_nonlinear(model,verbose)
order = 1
reffe = ReferenceFE(lagrangian,Float64,order)
V = FESpace(model,reffe;dirichlet_tags="boundary")

sol = [x -> x[1], x -> x[2], x -> x[1] + x[2], x -> 2.0*x[1]]
U1 = TrialFESpace(V,sol[1])
U2 = TrialFESpace(V,sol[2])
U3 = TrialFESpace(V,sol[3])
U4 = TrialFESpace(V,sol[4])

# Define weakforms
Ω = Triangulation(model)
= Measure(Ω,4*order)

F(u::Function) = x -> (u(x) + 1) * u(x)
F(u) = (u + 1) * u
dF(u,du) = 2.0 * u * du + du

j1((),u1,du1,dv1) = (dF(u1,du1) * dv1)dΩ
r1((),u1,v1) = ((F(u1) - F(sol[1])) * v1)dΩ

j2((u1,),(u2,u3),(du2,du3),(dv2,dv3)) = (u1 * dF(u2,du2) * dv2)dΩ + (dF(u3,du3) * dv3)dΩ
r2((u1,),(u2,u3),(v2,v3)) = (u1 * (F(u2) - F(sol[2])) * v2)dΩ + ((F(u3) - F(sol[3])) * v3)dΩ

j3((u1,(u2,u3)),u4,du4,dv4) = (u3 * dF(u4,du4) * dv4)dΩ
r3((u1,(u2,u3)),u4,v4) = (u3 * (F(u4) - F(sol[4])) * v4)dΩ

# Define solver: Each block will be solved with a LU solver
lsolver = LUSolver()
nlsolver = NewtonSolver(lsolver;rtol=1.e-10,verbose=verbose)
solver = StaggeredFESolver(fill(nlsolver,3))

# Create operator from full spaces
mfs = BlockMultiFieldStyle(3,(1,2,1))
X = MultiFieldFESpace([U1,U2,U3,U4];style=mfs)
Y = MultiFieldFESpace([V,V,V,V];style=mfs)
op = StaggeredNonlinearFEOperator([r1,r2,r3],[j1,j2,j3],X,Y)
xh = zero(X)
xh, cache = solve!(xh,solver,op);
test_solution(xh,sol,X,dΩ,verbose)

# Create operator from components
UB1, VB1 = U1, V
UB2, VB2 = MultiFieldFESpace([U2,U3]), MultiFieldFESpace([V,V])
UB3, VB3 = U4, V
op = StaggeredNonlinearFEOperator([r1,r2,r3],[j1,j2,j3],[UB1,UB2,UB3],[VB1,VB2,VB3])

xh = zero(X)
xh, cache = solve!(xh,solver,op,cache);
test_solution(xh,sol,X,dΩ,verbose)

return true
end

############################################################################################

np = (2,2)
ranks = DebugArray(LinearIndices((prod(np),)))
verbose = i_am_main(ranks)
model = CartesianDiscreteModel(ranks,np,(0,1,0,1),(4,4))

@testset "StaggeredAffineFEOperators" driver_affine(model,verbose)

model = CartesianDiscreteModel((0,1,0,1),(4,4))

order = 1
reffe = ReferenceFE(lagrangian,Float64,order)
V = FESpace(model,reffe;dirichlet_tags="boundary")

sol = [x -> x[1], x -> x[2], x -> x[1] + x[2], x -> x[1] - x[2]]
U1 = TrialFESpace(V,sol[1])
U2 = TrialFESpace(V,sol[2])
U3 = TrialFESpace(V,sol[3])
U4 = TrialFESpace(V,sol[4])

# Define weakforms
Ω = Triangulation(model)
= Measure(Ω,4*order)

F(u::Function) = x -> (u(x) + 1) * u(x)
F(u) = (u + 1) * u
dF(u,du) = 2.0 * u * du + du

j1((),u1,du1,dv1) = (dF(u1,du1) * dv1)dΩ
r1((),u1,v1) = ((F(u1) - F(sol[1])) * v1)dΩ

j2((u1,),(u2,u3),(du2,du3),(dv2,dv3)) = (u1 * dF(u2,du2) * dv2)dΩ + (dF(u3,du3) * dv3)dΩ
r2((u1,),(u2,u3),(v2,v3)) = (u1 * (F(u2) - F(sol[2])) * v2)dΩ + ((F(u3) - F(sol[3])) * v3)dΩ

j3((u1,(u2,u3)),u4,du4,dv4) = (u3 * dF(u4,du4) * dv4)dΩ
r3((u1,(u2,u3)),u4,v4) = (u3 * (F(u4) - F(sol[4])) * v4)dΩ

# Define solver: Each block will be solved with a LU solver
lsolver = LUSolver()
nlsolver = NewtonSolver(lsolver;rtol=1.e-10,verbose=verbose)
solver = StaggeredFESolver(fill(nlsolver,3))

# Create operator from full spaces
mfs = BlockMultiFieldStyle(3,(1,2,1))
X = MultiFieldFESpace([U1,U2,U3,U4];style=mfs)
Y = MultiFieldFESpace([V,V,V,V];style=mfs)
op = StaggeredNonlinearFEOperator([r1,r2,r3],[j1,j2,j3],X,Y)
xh = solve(solver,op)
test_solution(xh,sol,X,dΩ)
function driver(model,verbose)
@testset "StaggeredAffineFEOperators" driver_affine(model,verbose)
@testset "StaggeredNonlinearFEOperators" driver_affine(model,verbose)
end

# Create operator from components
UB1, VB1 = U1, V
UB2, VB2 = MultiFieldFESpace([U2,U3]), MultiFieldFESpace([V,V])
UB3, VB3 = U4, V
op = StaggeredNonlinearFEOperator([r1,r2,r3],[j1,j2,j3],[UB1,UB2,UB3],[VB1,VB2,VB3])
# Distributed
function main(distribute,parts)
ranks = distribute(LinearIndices((prod(parts),)))
model = CartesianDiscreteModel(ranks,parts,(0,1,0,1),(8,8))
driver(model,false)
end

xh = zero(X)
xh, cache = solve!(xh,solver,op);
test_solution(xh,sol,X,dΩ)
# Serial
function main()
model = CartesianDiscreteModel((0,1,0,1),(8,8))
driver(model,false)
end

end # module
9 changes: 9 additions & 0 deletions test/BlockSolvers/mpi/StaggeredFEOperatorsTests.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,9 @@
module StaggeredFEOperatorsMPITests
using MPI, PartitionedArrays
include("../StaggeredFEOperatorsTests.jl")

with_mpi() do distribute
StaggeredFEOperatorsTests.main(distribute,(2,2))
end

end
7 changes: 7 additions & 0 deletions test/BlockSolvers/seq/StaggeredFEOperatorsTests.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,7 @@
module StaggeredFEOperatorsSequentialTests
using PartitionedArrays
include("../StaggeredFEOperatorsTests.jl")

StaggeredFEOperatorsTests.main()

end
1 change: 1 addition & 0 deletions test/BlockSolvers/seq/runtests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2,3 +2,4 @@ using Test

@testset "BlockDiagonalSolvers" begin include("BlockDiagonalSolversTests.jl") end
@testset "BlockTriangularSolvers" begin include("BlockTriangularSolversTests.jl") end
@testset "StaggeredFEOperators" begin include("StaggeredFEOperatorsTests.jl") end

0 comments on commit 4563876

Please sign in to comment.