Merge pull request #70 from gridap/darcy_driver

Darcy driver

test/DarcyTests.jl

@@ -0,0 +1,61 @@
+using Gridap
+using GridapDistributed
+using GridapPETSc
+using GridapPETSc: PetscScalar, PetscInt, PETSC
+using SparseMatricesCSR
+using PartitionedArrays
+
+function main(parts)
+
+  if PETSC.MatMumpsSetIcntl_handle[] == C_NULL
+    @info "Skipping DarcyTests since petsc is not configured with mumps."
+    return nothing
+  end
+
+  domain = (0,1,0,1)
+  partition = (100,100)
+  model = CartesianDiscreteModel(parts,domain,partition)
+
+  k = 1
+  reffe_u = ReferenceFE(raviart_thomas,Float64,k)
+  reffe_p = ReferenceFE(lagrangian,Float64,k)
+
+  V = FESpace(model,reffe_u,dirichlet_tags=[5,6])
+  Q = FESpace(model,reffe_p,conformity=:L2)
+
+  uD = VectorValue(0.0,0.0)
+  U = TrialFESpace(V,uD)
+  P = TrialFESpace(Q)
+
+  Y = MultiFieldFESpace([V, Q])
+  X = MultiFieldFESpace([U, P])
+
+  Ω = Interior(model)
+  degree = 2
+  dΩ = Measure(Ω,degree)
+
+  neumanntags = [8,]
+  Γ = Boundary(model,tags=neumanntags)
+  dΓ = Measure(Γ,degree)
+
+  kinv1 = TensorValue(1.0,0.0,0.0,1.0)
+
+  σ(u) =kinv1⋅u
+
+  a((u,p), (v,q)) = ∫(v⋅(σ∘u) - (∇⋅v)*p + q*(∇⋅u) + 0*p*q)dΩ
+  n_Γ = get_normal_vector(Γ)
+  h = -1.0
+
+  b((v,q)) = ∫((v⋅n_Γ)*h)dΓ
+
+  op = AffineFEOperator(a,b,X,Y)
+
+  options = "-ksp_error_if_not_converged true -ksp_converged_reason -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps"
+  xh = GridapPETSc.with(args=split(options)) do
+    ls = PETScLinearSolver()
+    xh = solve(ls,op)
+  end
+  uh, ph = xh
+
+end
+

test/PLaplacianTests.jl

@@ -13,10 +13,10 @@ function mysnessetup(snes)
   mumpsmat = Ref{GridapPETSc.PETSC.Mat}()
   @check_error_code GridapPETSc.PETSC.SNESSetFromOptions(snes[])
   @check_error_code GridapPETSc.PETSC.SNESGetKSP(snes[],ksp)
-  @check_error_code GridapPETSc.PETSC.KSPView(ksp[],C_NULL)
+  #@check_error_code GridapPETSc.PETSC.KSPView(ksp[],C_NULL)
   @check_error_code GridapPETSc.PETSC.KSPSetType(ksp[],GridapPETSc.PETSC.KSPPREONLY)
   @check_error_code GridapPETSc.PETSC.KSPGetPC(ksp[],pc)
-  @check_error_code GridapPETSc.PETSC.PCView(pc[],C_NULL)
+  #@check_error_code GridapPETSc.PETSC.PCView(pc[],C_NULL)
   @check_error_code GridapPETSc.PETSC.PCSetType(pc[],GridapPETSc.PETSC.PCLU)
   @check_error_code GridapPETSc.PETSC.PCFactorSetMatSolverType(pc[],GridapPETSc.PETSC.MATSOLVERMUMPS)
   @check_error_code GridapPETSc.PETSC.PCFactorSetUpMatSolverType(pc[])
@@ -36,9 +36,9 @@ end
 
 function main(parts,solver)
   if solver == :mumps
-    options = "-snes_type newtonls -snes_linesearch_type basic  -snes_linesearch_damping 1.0 -snes_rtol 1.0e-14 -snes_atol 0.0 -snes_monitor -snes_converged_reason"
+    options = "-snes_type newtonls -snes_linesearch_type basic  -snes_linesearch_damping 1.0 -snes_rtol 1.0e-14 -snes_atol 0.0 -snes_monitor -snes_converged_reason -ksp_converged_reason -ksp_error_if_not_converged true"
   elseif solver == :gmres
-    options = "-snes_type newtonls -snes_linesearch_type basic  -snes_linesearch_damping 1.0 -snes_rtol 1.0e-14 -snes_atol 0.0 -snes_monitor -pc_type jacobi -ksp_type gmres -ksp_monitor -snes_converged_reason"
+    options = "-snes_type newtonls -snes_linesearch_type basic  -snes_linesearch_damping 1.0 -snes_rtol 1.0e-14 -snes_atol 0.0 -snes_monitor -pc_type jacobi -ksp_type gmres -snes_converged_reason -ksp_converged_reason -ksp_error_if_not_converged true"
   else
     error()
   end

test/PartitionedArraysTests.jl

@@ -52,7 +52,7 @@ function partitioned_tests(parts)
 
   end
 
-  GridapPETSc.Init(args=split("-ksp_type gmres -ksp_monitor -pc_type jacobi"))
+  GridapPETSc.Init(args=split("-ksp_type gmres -ksp_converged_reason -ksp_error_if_not_converged true -pc_type jacobi"))
 
   function test_get_local_vector(v::PVector,x::PETScVector)
     if (get_backend(v.values)==mpi)

test/PoissonTests.jl

@@ -34,9 +34,9 @@ end
 
 function main(parts,solver)
   if solver == :mumps
-    options = "-info  -ksp_error_if_not_converged true"
+    options = "-ksp_error_if_not_converged true -ksp_converged_reason"
   elseif solver == :cg
-    options = "-info -pc_type jacobi -ksp_type cg -ksp_monitor -ksp_rtol 1.0e-12"
+    options = "-pc_type jacobi -ksp_type cg -ksp_error_if_not_converged true -ksp_converged_reason -ksp_rtol 1.0e-12"
   else
     error()
   end

test/mpi/DarcyTests.jl

@@ -0,0 +1,3 @@
+include("../DarcyTests.jl")
+nparts = (2,2)
+prun(main,mpi,nparts)

test/mpi/DarcyTestsRun.jl

@@ -0,0 +1,4 @@
+module DarcyTestsRun
+ include("mpiexec.jl")
+ run_mpi_driver(procs=4,file="DarcyTests.jl")
+end # module

test/mpi/GCTests.jl

@@ -23,7 +23,7 @@ function report_memory_and_random_gc(parts)
 end
 
 function main_bis(parts)
-  main(parts,SubAssembledRows())
+  main(parts,:gmres,SubAssembledRows())
   report_memory_and_random_gc(parts)
 end
 

test/mpi/runtests.jl

@@ -4,4 +4,5 @@ using Test
 @time @testset "PLaplacianTests" begin include("PLaplacianTestsRun.jl") end
 @time @testset "GCTests" begin include("GCTestsRun.jl") end
 @time @testset "PoissonTests" begin include("PoissonTestsRun.jl") end
+@time @testset "DarcyTests" begin include("DarcyTestsRun.jl") end
 end # module

test/sequential/DarcyDriver.jl

@@ -0,0 +1,68 @@
+module DarcyDriver
+
+using Gridap
+using GridapPETSc
+using GridapPETSc: PetscScalar, PetscInt, PETSC
+using SparseMatricesCSR
+
+options = "-ksp_converged_reason -ksp_monitor -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps"
+
+out = GridapPETSc.with(args=split(options)) do
+
+  if PETSC.MatMumpsSetIcntl_handle[] == C_NULL
+    @info "Skipping DarcyDriver since petsc is not configured with mumps."
+    return nothing
+  end
+
+  domain = (0,1,0,1)
+  partition = (100,100)
+  model = CartesianDiscreteModel(domain,partition)
+
+  k = 1
+  reffe_u = ReferenceFE(raviart_thomas,Float64,k)
+  reffe_p = ReferenceFE(lagrangian,Float64,k)
+
+  V = FESpace(model,reffe_u,dirichlet_tags=[5,6])
+  Q = FESpace(model,reffe_p,conformity=:L2)
+
+  uD = VectorValue(0.0,0.0)
+  U = TrialFESpace(V,uD)
+  P = TrialFESpace(Q)
+
+  Y = MultiFieldFESpace([V, Q])
+  X = MultiFieldFESpace([U, P])
+
+  Ω = Interior(model)
+  degree = 2
+  dΩ = Measure(Ω,degree)
+
+  neumanntags = [8,]
+  Γ = Boundary(model,tags=neumanntags)
+  dΓ = Measure(Γ,degree)
+
+  kinv1 = TensorValue(1.0,0.0,0.0,1.0)
+  kinv2 = TensorValue(100.0,90.0,90.0,100.0)
+
+  function σ(x,u)
+     if ((abs(x[1]-0.5) <= 0.1) && (abs(x[2]-0.5) <= 0.1))
+        return kinv2⋅u
+     else
+        return kinv1⋅u
+     end
+  end
+
+  px = get_physical_coordinate(Ω)
+  a((u,p), (v,q)) = ∫(v⋅(σ∘(px,u)) - (∇⋅v)*p + q*(∇⋅u) + 0*p*q)dΩ
+  n_Γ = get_normal_vector(Γ)
+  h = -1.0
+
+  b((v,q)) = ∫((v⋅n_Γ)*h)dΓ
+
+  op = AffineFEOperator(a,b,X,Y)
+  ls = PETScLinearSolver()
+  xh = solve(ls,op)
+  uh, ph = xh
+
+end
+
+end # module

test/sequential/runtests.jl

@@ -18,6 +18,8 @@ using Test
 
 @time @testset "ElasticityDriver" begin include("ElasticityDriver.jl") end
 
+@time @testset "DarcyDriver" begin include("DarcyDriver.jl") end
+
 @time @testset "PLaplacianDriver" begin include("PLaplacianDriver.jl") end
 
 end # module