From fd26e5088f5f81ecdf2d361540900ac4e7df9101 Mon Sep 17 00:00:00 2001 From: Oriol Colomes Date: Thu, 18 Mar 2021 17:30:30 +0100 Subject: [PATCH 01/15] Added TransientCellField and TransientFEBasis --- src/TransientFETools/TransientCellField.jl | 34 ++++++++++++++++++++++ src/TransientFETools/TransientFETools.jl | 6 ++++ 2 files changed, 40 insertions(+) create mode 100644 src/TransientFETools/TransientCellField.jl diff --git a/src/TransientFETools/TransientCellField.jl b/src/TransientFETools/TransientCellField.jl new file mode 100644 index 0000000..7c167dc --- /dev/null +++ b/src/TransientFETools/TransientCellField.jl @@ -0,0 +1,34 @@ +# Transient CellField +struct TransientCellField{A} <: CellField + cellfield::A + derivatives::Tuple{Vararg{A}} +end + +# CellField methods +get_data(f::TransientCellField) = get_data(f.cellfield) +get_triangulation(f::TransientCellField) = get_triangulation(f.cellfield) +DomainStyle(::Type{<:TransientCellField{A}}) where A = DomainStyle(A) +gradient(f::TransientCellField) = gradient(f.cellfield) +∇∇(f::TransientCellField) = ∇∇(f.cellfield) +change_domain(f::TransientCellField,trian::Triangulation,target_domain::DomainStyle) = change_domain(f.cellfield,trian,target_domain) + +# Transient FEBasis +struct TransientFEBasis{A} <: FEBasis + febasis::A + derivatives::Tuple{Vararg{A}} +end + +# FEBasis methods +get_data(f::TransientFEBasis) = get_data(f.febasis) +get_triangulation(f::TransientFEBasis) = get_triangulation(f.febasis) +DomainStyle(::Type{<:TransientFEBasis{A}}) where A = DomainStyle(A) +BasisStyle(::Type{<:TransientFEBasis{A}}) where A = BasisStyle(A) +gradient(f::TransientFEBasis) = gradient(f.febasis) +∇∇(f::TransientFEBasis) = ∇∇(f.febasis) +change_domain(f::TransientFEBasis,trian::Triangulation,target_domain::DomainStyle) = change_domain(f.febasis,trian,target_domain) + +# Time derivative +function ∂t(f::Union{TransientCellField,TransientFEBasis}) + cellfield, derivatives = first_and_tail(f.derivatives) + TransientCellField(cellfield,derivatives) +end \ No newline at end of file diff --git a/src/TransientFETools/TransientFETools.jl b/src/TransientFETools/TransientFETools.jl index a82137b..e9a48ad 100644 --- a/src/TransientFETools/TransientFETools.jl +++ b/src/TransientFETools/TransientFETools.jl @@ -85,6 +85,10 @@ export test_transient_fe_solution export Transient2ndOrderFEOperator +using Gridap.CellData: CellField +using Gridap.CellData: DomainStyle +using Gridap.FESpaces: FEBasis + include("TransientFESpaces.jl") include("TransientFEOperators.jl") @@ -95,6 +99,8 @@ include("TransientFESolvers.jl") include("TransientFESolutions.jl") +include("TransientCellField.jl") + export FETerm function FETerm(args...) Helpers.@unreachable """\n From d8793fc9621b5de9b80a43209f5d93d826822e0c Mon Sep 17 00:00:00 2001 From: Oriol Colomes Date: Thu, 18 Mar 2021 21:44:31 +0100 Subject: [PATCH 02/15] Implemented new API --- src/TransientFETools/ODEOperatorInterfaces.jl | 21 ++++--- src/TransientFETools/TransientCellField.jl | 4 +- src/TransientFETools/TransientFEOperators.jl | 59 ++++++++++--------- src/TransientFETools/TransientFETools.jl | 15 +++-- 4 files changed, 56 insertions(+), 43 deletions(-) diff --git a/src/TransientFETools/ODEOperatorInterfaces.jl b/src/TransientFETools/ODEOperatorInterfaces.jl index 47b932e..619bd01 100644 --- a/src/TransientFETools/ODEOperatorInterfaces.jl +++ b/src/TransientFETools/ODEOperatorInterfaces.jl @@ -61,10 +61,11 @@ function residual!( xhF::Tuple{Vararg{AbstractVector}}, ode_cache) Xh, = ode_cache - xh = () - for i in 1:get_order(op)+1 - xh = (xh...,EvaluationFunction(Xh[i],xhF[i])) + dxh = () + for i in 2:get_order(op)+1 + dxh = (dxh...,EvaluationFunction(Xh[i],xhF[i])) end + xh=TransientCellField(EvaluationFunction(Xh[1],xhF[1]),dxh) residual!(b,op.feop,t,xh,ode_cache) end @@ -86,10 +87,11 @@ function jacobian!( γᵢ::Real, ode_cache) Xh, = ode_cache - xh = () - for i in 1:get_order(op)+1 - xh = (xh...,EvaluationFunction(Xh[i],xhF[i])) + dxh = () + for i in 2:get_order(op)+1 + dxh = (dxh...,EvaluationFunction(Xh[i],xhF[i])) end + xh=TransientCellField(EvaluationFunction(Xh[1],xhF[1]),dxh) jacobian!(A,op.feop,t,xh,i,γᵢ,ode_cache) end @@ -104,9 +106,10 @@ function jacobians!( γ::Tuple{Vararg{Real}}, ode_cache) Xh, = ode_cache - xh = () - for i in 1:get_order(op)+1 - xh = (xh...,EvaluationFunction(Xh[i],xhF[i])) + dxh = () + for i in 2:get_order(op)+1 + dxh = (dxh...,EvaluationFunction(Xh[i],xhF[i])) end + xh=TransientCellField(EvaluationFunction(Xh[1],xhF[1]),dxh) jacobians!(J,op.feop,t,xh,γ,ode_cache) end diff --git a/src/TransientFETools/TransientCellField.jl b/src/TransientFETools/TransientCellField.jl index 7c167dc..6d1dc79 100644 --- a/src/TransientFETools/TransientCellField.jl +++ b/src/TransientFETools/TransientCellField.jl @@ -31,4 +31,6 @@ change_domain(f::TransientFEBasis,trian::Triangulation,target_domain::DomainStyl function ∂t(f::Union{TransientCellField,TransientFEBasis}) cellfield, derivatives = first_and_tail(f.derivatives) TransientCellField(cellfield,derivatives) -end \ No newline at end of file +end + +∂tt(f::Union{TransientCellField,TransientFEBasis}) = ∂t(∂t(f::Union{TransientCellField,TransientFEBasis})) \ No newline at end of file diff --git a/src/TransientFETools/TransientFEOperators.jl b/src/TransientFETools/TransientFEOperators.jl index 6f7c007..38dd440 100644 --- a/src/TransientFETools/TransientFEOperators.jl +++ b/src/TransientFETools/TransientFEOperators.jl @@ -112,18 +112,18 @@ end function TransientConstantFEOperator(m::Function,a::Function,b::Function, trial,test) - res(t,(u,ut),v) = m(ut,v) + a(u,v) - b(v) - jac(t,(u,ut),du,v) = a(du,v) - jac_t(t,(u,ut),dut,v) = m(dut,v) + res(t,u,v) = m(∂t(u),v) + a(u,v) - b(v) + jac(t,u,du,v) = a(du,v) + jac_t(t,u,dut,v) = m(dut,v) assem_t = SparseMatrixAssembler(trial,test) TransientFEOperatorFromWeakForm{Constant}(res,(jac,jac_t),assem_t,(trial,∂t(trial)),test,1) end function TransientAffineFEOperator(m::Function,a::Function,b::Function, trial,test) - res(t,(u,ut),v) = m(t,ut,v) + a(t,u,v) - b(t,v) - jac(t,(u,ut),du,v) = a(t,du,v) - jac_t(t,(u,ut),dut,v) = m(t,dut,v) + res(t,u,v) = m(t,∂t(u),v) + a(t,u,v) - b(t,v) + jac(t,u,du,v) = a(t,du,v) + jac_t(t,u,dut,v) = m(t,dut,v) assem_t = SparseMatrixAssembler(trial,test) TransientFEOperatorFromWeakForm{Affine}(res,(jac,jac_t),assem_t,(trial,∂t(trial)),test,1) end @@ -137,10 +137,10 @@ end function TransientConstantFEOperator(m::Function,c::Function,a::Function,b::Function, trial,test) - res(t,(u,ut,utt),v) = m(utt,v) + c(ut,v) + a(u,v) - b(v) - jac(t,(u,ut,utt),du,v) = a(du,v) - jac_t(t,(u,ut,utt),dut,v) = c(dut,v) - jac_tt(t,(u,ut,utt),dutt,v) = m(dutt,v) + res(t,u,v) = m(∂tt(u),v) + c(∂t(u),v) + a(u,v) - b(v) + jac(t,u,du,v) = a(du,v) + jac_t(t,u,dut,v) = c(dut,v) + jac_tt(t,u,dutt,v) = m(dutt,v) assem_t = SparseMatrixAssembler(trial,test) trial_t = ∂t(trial) trial_tt = ∂t(trial_t) @@ -150,10 +150,10 @@ end function TransientAffineFEOperator(m::Function,c::Function,a::Function,b::Function, trial,test) - res(t,(u,ut,utt),v) = m(t,utt,v) + c(t,ut,v) + a(t,u,v) - b(t,v) - jac(t,(u,ut,utt),du,v) = a(t,du,v) - jac_t(t,(u,ut,utt),dut,v) = c(t,dut,v) - jac_tt(t,(u,ut,utt),dutt,v) = m(t,dutt,v) + res(t,u,v) = m(t,∂tt(u),v) + c(t,∂t(u),v) + a(t,u,v) - b(t,v) + jac(t,u,du,v) = a(t,du,v) + jac_t(t,u,dut,v) = c(t,dut,v) + jac_tt(t,u,dutt,v) = m(t,dutt,v) assem_t = SparseMatrixAssembler(trial,test) trial_t = ∂t(trial) trial_tt = ∂t(trial_t) @@ -204,10 +204,11 @@ get_order(op::TransientFEOperatorFromWeakForm) = op.order function allocate_residual(op::TransientFEOperatorFromWeakForm,uh::FEFunction,cache) V = get_test(op) v = get_cell_shapefuns(V) - xh = () - for i in 1:get_order(op)+1 - xh = (xh...,uh) + dxh = () + for i in 1:get_order(op) + dxh = (dxh...,uh) end + xh = TransientCellField(uh,dxh) vecdata = collect_cell_vector(V,op.res(0.0,xh,v)) allocate_vector(op.assem_t,vecdata) end @@ -216,7 +217,7 @@ function residual!( b::AbstractVector, op::TransientFEOperatorFromWeakForm, t::Real, - xh::Tuple{Vararg{FEFunction}}, + xh::TransientCellField,#Tuple{Vararg{FEFunction}}, cache) V = get_test(op) v = get_cell_shapefuns(V) @@ -226,10 +227,11 @@ function residual!( end function allocate_jacobian(op::TransientFEOperatorFromWeakForm,uh::FEFunction,cache) - xh = () - for i in 1:get_order(op)+1 - xh = (xh...,uh) + dxh = () + for i in 1:get_order(op) + dxh = (dxh...,uh) end + xh = TransientCellField(uh,dxh) _matdata = () for i in 1:get_order(op)+1 _matdata = (_matdata...,matdata_jacobian(op,0.0,xh,i,0.0)) @@ -242,7 +244,7 @@ function jacobian!( A::AbstractMatrix, op::TransientFEOperatorFromWeakForm, t::Real, - xh::Tuple{Vararg{FEFunction}}, + xh::TransientCellField,#Tuple{Vararg{FEFunction}}, i::Integer, γᵢ::Real, cache) @@ -255,7 +257,7 @@ function jacobians!( A::AbstractMatrix, op::TransientFEOperatorFromWeakForm, t::Real, - xh::Tuple{Vararg{FEFunction}}, + xh::TransientCellField,#Tuple{Vararg{FEFunction}}, γ::Tuple{Vararg{Real}}, cache) _matdata = () @@ -288,7 +290,7 @@ end function matdata_jacobian( op::TransientFEOperatorFromWeakForm, t::Real, - xh::Tuple{Vararg{FEFunction}}, + xh::TransientCellField,#Tuple{Vararg{FEFunction}}, i::Integer, γᵢ::Real) Uh = evaluate(get_trial(op),nothing) @@ -311,15 +313,16 @@ function test_transient_fe_operator(op::TransientFEOperator,uh) @test isa(U0,FESpace) r = allocate_residual(op,uh,cache) @test isa(r,AbstractVector) - residual!(r,op,0.0,(uh,uh),cache) + xh = TransientCellField(uh,(uh,)) + residual!(r,op,0.0,xh,cache) @test isa(r,AbstractVector) J = allocate_jacobian(op,uh,cache) @test isa(J,AbstractMatrix) - jacobian!(J,op,0.0,(uh,uh),1,1.0,cache) + jacobian!(J,op,0.0,xh,1,1.0,cache) @test isa(J,AbstractMatrix) - jacobian!(J,op,0.0,(uh,uh),2,1.0,cache) + jacobian!(J,op,0.0,xh,2,1.0,cache) @test isa(J,AbstractMatrix) - jacobians!(J,op,0.0,(uh,uh),(1.0,1.0),cache) + jacobians!(J,op,0.0,xh,(1.0,1.0),cache) @test isa(J,AbstractMatrix) cache = update_cache!(cache,op,0.0) true diff --git a/src/TransientFETools/TransientFETools.jl b/src/TransientFETools/TransientFETools.jl index e9a48ad..3823da9 100644 --- a/src/TransientFETools/TransientFETools.jl +++ b/src/TransientFETools/TransientFETools.jl @@ -83,14 +83,21 @@ import GridapODEs.ODETools: GenericODESolution import Base: iterate export test_transient_fe_solution -export Transient2ndOrderFEOperator - +export TransientCellField using Gridap.CellData: CellField -using Gridap.CellData: DomainStyle using Gridap.FESpaces: FEBasis +import Gridap.CellData: get_data +import Gridap.CellData: get_triangulation +import Gridap.CellData: DomainStyle +import Gridap.CellData: gradient +import Gridap.CellData: ∇∇ +import Gridap.CellData: change_domain +import Gridap.FESpaces: BasisStyle include("TransientFESpaces.jl") +include("TransientCellField.jl") + include("TransientFEOperators.jl") include("ODEOperatorInterfaces.jl") @@ -99,8 +106,6 @@ include("TransientFESolvers.jl") include("TransientFESolutions.jl") -include("TransientCellField.jl") - export FETerm function FETerm(args...) Helpers.@unreachable """\n From df0ad71314fa80881655dd62362719c2a4913bf4 Mon Sep 17 00:00:00 2001 From: Oriol Colomes Date: Thu, 18 Mar 2021 21:45:19 +0100 Subject: [PATCH 03/15] TransientFETests passing --- .../Transient2ndOrderFEOperatorsTests.jl | 8 ++++---- .../TransientFEOperatorsTests.jl | 6 +++--- test/TransientFEsTests/TransientFETests.jl | 20 ++++++++++--------- 3 files changed, 18 insertions(+), 16 deletions(-) diff --git a/test/TransientFEsTests/Transient2ndOrderFEOperatorsTests.jl b/test/TransientFEsTests/Transient2ndOrderFEOperatorsTests.jl index a2966ec..684a3be 100644 --- a/test/TransientFEsTests/Transient2ndOrderFEOperatorsTests.jl +++ b/test/TransientFEsTests/Transient2ndOrderFEOperatorsTests.jl @@ -34,10 +34,10 @@ c(ut,v) = ∫(v*ut)dΩ a(u,v) = ∫(∇(v)⊙∇(u))dΩ b(v,t) = ∫(v*f(t))dΩ -res(t,(u,ut,utt),v) = m(utt,v) + c(ut,v) + a(u,v) - b(v,t) -jac(t,(u,ut,utt),du,v) = a(du,v) -jac_t(t,(u,ut,utt),dut,v) = c(dut,v) -jac_tt(t,(u,ut,utt),dutt,v) = m(dutt,v) +res(t,u,v) = m(∂tt(u),v) + c(∂t(u),v) + a(u,v) - b(v,t) +jac(t,u,du,v) = a(du,v) +jac_t(t,u,dut,v) = c(dut,v) +jac_tt(t,u,dutt,v) = m(dutt,v) op = TransientFEOperator(res,jac,jac_t,jac_tt,U,V0) diff --git a/test/TransientFEsTests/TransientFEOperatorsTests.jl b/test/TransientFEsTests/TransientFEOperatorsTests.jl index ed12332..28133c6 100644 --- a/test/TransientFEsTests/TransientFEOperatorsTests.jl +++ b/test/TransientFEsTests/TransientFEOperatorsTests.jl @@ -45,9 +45,9 @@ a(u,v) = ∫(∇(v)⊙∇(u))dΩ m(u,v) = ∫(v*u)dΩ b(v,t) = ∫(v*f(t))dΩ -res(t,(u,ut),v) = a(u,v) + m(ut,v) - b(v,t) -jac(t,(u,ut),du,v) = a(du,v) -jac_t(t,(u,ut),dut,v) = m(dut,v) +res(t,u,v) = a(u,v) + m(∂t(u),v) - b(v,t) +jac(t,u,du,v) = a(du,v) +jac_t(t,u,dut,v) = m(dut,v) op = TransientFEOperator(res,jac,jac_t,U,V0) diff --git a/test/TransientFEsTests/TransientFETests.jl b/test/TransientFEsTests/TransientFETests.jl index bbdefe0..5877faa 100644 --- a/test/TransientFEsTests/TransientFETests.jl +++ b/test/TransientFEsTests/TransientFETests.jl @@ -67,9 +67,9 @@ dΩ = Measure(Ω,degree) a(u,v) = ∫(∇(v)⋅∇(u))dΩ b(v,t) = ∫(v*f(t))dΩ -res(t,(u,ut),v) = a(u,v) + ∫(ut*v)dΩ - b(v,t) -jac(t,(u,ut),du,v) = a(du,v) -jac_t(t,(u,ut),dut,v) = ∫(dut*v)dΩ +res(t,u,v) = a(u,v) + ∫(∂t(u)*v)dΩ - b(v,t) +jac(t,u,du,v) = a(du,v) +jac_t(t,u,dut,v) = ∫(dut*v)dΩ U0 = U(0.0) _res(u,v) = a(u,v) + 10.0*∫(u*v)dΩ - b(v,0.0) @@ -91,9 +91,10 @@ cache = allocate_cache(odeop) r = allocate_residual(op,uh,cache) J = allocate_jacobian(op,uh,cache) uh10 = interpolate_everywhere(0.0,U0)#10.0) -residual!(r,op,0.0,(uh,uh10),cache) -jacobian!(J,op,1.0,(uh,uh10),1,1.0,cache) -jacobian!(J,op,1.0,(uh,uh10),2,10.0,cache) +xh = TransientCellField(uh,(uh10,)) +residual!(r,op,0.0,xh,cache) +jacobian!(J,op,1.0,xh,1,1.0,cache) +jacobian!(J,op,1.0,xh,2,10.0,cache) @test all(r.≈_r) @test all(J.≈_J) @@ -116,9 +117,10 @@ odes = ThetaMethod(nls,dt,1.0) solver = TransientFESolver(odes) # Return a specialization of TransientFESolver @test test_transient_fe_solver(solver,op,uh0,t0,tF) -residual!(r,op,0.1,(uh,uh),cache) -jacobian!(J,op,1.0,(uh,uh10),1,1.0,cache) -jacobian!(J,op,1.0,(uh,uh10),2,10.0,cache) +xh = TransientCellField(uh,(uh,)) +residual!(r,op,0.1,xh,cache) +jacobian!(J,op,1.0,xh,1,1.0,cache) +jacobian!(J,op,1.0,xh,2,10.0,cache) u0 = get_free_dof_values(uh0) odes From bdbba6696c66ce685ce4fe1551e19ee5fd405d51 Mon Sep 17 00:00:00 2001 From: Oriol Colomes Date: Thu, 25 Mar 2021 21:53:23 +0100 Subject: [PATCH 04/15] Added functions for transient MultiFieldCellField --- src/TransientFETools/TransientCellField.jl | 29 ++++++++++++++++++++++ src/TransientFETools/TransientFETools.jl | 1 + 2 files changed, 30 insertions(+) diff --git a/src/TransientFETools/TransientCellField.jl b/src/TransientFETools/TransientCellField.jl index 6d1dc79..aca7497 100644 --- a/src/TransientFETools/TransientCellField.jl +++ b/src/TransientFETools/TransientCellField.jl @@ -12,6 +12,35 @@ gradient(f::TransientCellField) = gradient(f.cellfield) ∇∇(f::TransientCellField) = ∇∇(f.cellfield) change_domain(f::TransientCellField,trian::Triangulation,target_domain::DomainStyle) = change_domain(f.cellfield,trian,target_domain) +# MultiFieldCellField methods +num_fields(f::TransientCellField{A}) where A = num_fields(f.cellfield) +function Base.getindex(f::TransientCellField{A},i::Integer) where A + singleCellfield = getindex(f.single_cellfields,i) + singleDerivatives = (getindex(i_derivatives,i) for i_derivatives in f.derivatives) + TransientCellField(singleCellfield,singleDerivatives) +end +function Base.iterate(f::TransientCellField{A}) where A + cellField, state1 = iterate(f.cellfield) + derivatives = () + state2 = [] + for i_derivatives in f.derivatives + single_derivative, single_state = iterate(i_derivatives) + derivatives = (derivatives...,single_derivative) + push!(state2,single_state) + end + TransientCellField(cellField,derivatives),(state1,state2) +end +function Base.iterate(f::TransientCellField{A},state) where A + state1, state2 = state + cellField, state1 = iterate(f.cellfield,state[1]) + derivatives = () + for (i,i_derivatives) in enumerate(f.derivatives) + single_derivative, state2[i] = iterate(i_derivatives) + derivatives = (derivatives...,single_derivative) + end + TransientCellField(cellField,derivatives),(state1, state2) +end + # Transient FEBasis struct TransientFEBasis{A} <: FEBasis febasis::A diff --git a/src/TransientFETools/TransientFETools.jl b/src/TransientFETools/TransientFETools.jl index 3823da9..2d7f097 100644 --- a/src/TransientFETools/TransientFETools.jl +++ b/src/TransientFETools/TransientFETools.jl @@ -85,6 +85,7 @@ export test_transient_fe_solution export TransientCellField using Gridap.CellData: CellField +using Gridap.MultiField: MultiFieldCellField using Gridap.FESpaces: FEBasis import Gridap.CellData: get_data import Gridap.CellData: get_triangulation From b04e25f0675d1aec10a74b677fa326c89733ac41 Mon Sep 17 00:00:00 2001 From: Oriol Colomes Date: Thu, 25 Mar 2021 21:53:49 +0100 Subject: [PATCH 05/15] Modified tests --- test/TransientFEsTests/BoundaryHeatEquationTests.jl | 6 +++--- test/TransientFEsTests/DGHeatEquationTests.jl | 6 +++--- test/TransientFEsTests/ForwardEulerHeatEquationTests.jl | 6 +++--- test/TransientFEsTests/HeatEquationAutoDiffTests.jl | 6 +++--- test/TransientFEsTests/HeatEquationTests.jl | 6 +++--- test/TransientFEsTests/HeatVectorEquationTests.jl | 6 +++--- test/TransientFEsTests/StokesEquationAutoDiffTests.jl | 6 +++--- test/TransientFEsTests/StokesEquationTests.jl | 6 +++--- test/TransientFEsTests/VectorHeatEquationTests.jl | 8 ++++---- 9 files changed, 28 insertions(+), 28 deletions(-) diff --git a/test/TransientFEsTests/BoundaryHeatEquationTests.jl b/test/TransientFEsTests/BoundaryHeatEquationTests.jl index f07f76f..9f1afce 100644 --- a/test/TransientFEsTests/BoundaryHeatEquationTests.jl +++ b/test/TransientFEsTests/BoundaryHeatEquationTests.jl @@ -54,9 +54,9 @@ b(v,t) = ∫(v*f(t))dΩ m(ut,v) = ∫(ut*v)dΩ b_Γ(v,t) = ∫(v*(∇(u(t))⋅nb))dΓ -res(t,(u,ut),v) = a(u,v) + m(ut,v) - b(v,t) - b_Γ(v,t) -jac(t,(u,ut),du,v) = a(du,v) -jac_t(t,(u,ut),dut,v) = m(dut,v) +res(t,u,v) = a(u,v) + m(∂t(u),v) - b(v,t) - b_Γ(v,t) +jac(t,u,du,v) = a(du,v) +jac_t(t,u,dut,v) = m(dut,v) op = TransientFEOperator(res,jac,jac_t,U,V0) diff --git a/test/TransientFEsTests/DGHeatEquationTests.jl b/test/TransientFEsTests/DGHeatEquationTests.jl index 3cfaf1a..31c7329 100644 --- a/test/TransientFEsTests/DGHeatEquationTests.jl +++ b/test/TransientFEsTests/DGHeatEquationTests.jl @@ -54,9 +54,9 @@ b_Γ(v,t) = ∫( (γ/h)*v*u(t) - (∇(v)⋅nb)*u(t) )dΓ a_Λ(u,v) = ∫( (γ/h)*jump(v*ns)⊙jump(u*ns) - jump(v*ns)⊙mean(∇(u)) - mean(∇(v))⊙jump(u*ns) )dΛ -res(t,(u,ut),v) = a(u,v) + m(ut,v) + a_Γ(u,v) + a_Λ(u,v) - b(v,t) - b_Γ(v,t) -jac(t,(u,ut),du,v) = a(du,v) + a_Γ(du,v) + a_Λ(du,v) -jac_t(t,(u,ut),dut,v) = m(dut,v) +res(t,u,v) = a(u,v) + m(∂t(u),v) + a_Γ(u,v) + a_Λ(u,v) - b(v,t) - b_Γ(v,t) +jac(t,u,du,v) = a(du,v) + a_Γ(du,v) + a_Λ(du,v) +jac_t(t,u,dut,v) = m(dut,v) op = TransientFEOperator(res,jac,jac_t,U,V0) diff --git a/test/TransientFEsTests/ForwardEulerHeatEquationTests.jl b/test/TransientFEsTests/ForwardEulerHeatEquationTests.jl index 63427bc..5062111 100644 --- a/test/TransientFEsTests/ForwardEulerHeatEquationTests.jl +++ b/test/TransientFEsTests/ForwardEulerHeatEquationTests.jl @@ -43,9 +43,9 @@ quad = CellQuadrature(trian,degree) a(u,v) = ∇(v)⋅∇(u) b(v,t) = v*f(t) -res(t,(u,ut),v) = a(u,v) + ut*v - b(v,t) -jac(t,(u,ut),du,v) = a(du,v) -jac_t(t,(u,ut),dut,v) = dut*v +res(t,u,v) = a(u,v) + ∂t(u)*v - b(v,t) +jac(t,u,du,v) = a(du,v) +jac_t(t,u,dut,v) = dut*v t_Ω = FETerm(res,jac,jac_t,trian,quad) op = TransientFEOperator(U,V0,t_Ω) diff --git a/test/TransientFEsTests/HeatEquationAutoDiffTests.jl b/test/TransientFEsTests/HeatEquationAutoDiffTests.jl index db15846..5f3e975 100644 --- a/test/TransientFEsTests/HeatEquationAutoDiffTests.jl +++ b/test/TransientFEsTests/HeatEquationAutoDiffTests.jl @@ -41,9 +41,9 @@ dΩ = Measure(Ω,degree) a(u,v) = ∫(∇(v)⋅∇(u))dΩ b(v,t) = ∫(v*f(t))dΩ -res(t,(u,ut),v) = a(u,v) + ∫(ut*v)dΩ - b(v,t) -jac(t,(u,ut),du,v) = a(du,v) -jac_t(t,(u,ut),dut,v) = ∫(dut*v)dΩ +res(t,u,v) = a(u,v) + ∫(∂t(u)*v)dΩ - b(v,t) +jac(t,u,du,v) = a(du,v) +jac_t(t,u,dut,v) = ∫(dut*v)dΩ U₀ = evaluate(U,nothing) dv = get_cell_shapefuns(V0) diff --git a/test/TransientFEsTests/HeatEquationTests.jl b/test/TransientFEsTests/HeatEquationTests.jl index 70545a6..8dad651 100644 --- a/test/TransientFEsTests/HeatEquationTests.jl +++ b/test/TransientFEsTests/HeatEquationTests.jl @@ -40,9 +40,9 @@ dΩ = Measure(Ω,degree) a(u,v) = ∫(∇(v)⋅∇(u))dΩ b(v,t) = ∫(v*f(t))dΩ -res(t,(u,ut),v) = a(u,v) + ∫(ut*v)dΩ - b(v,t) -jac(t,(u,ut),du,v) = a(du,v) -jac_t(t,(u,ut),dut,v) = ∫(dut*v)dΩ +res(t,u,v) = a(u,v) + ∫(∂t(u)*v)dΩ - b(v,t) +jac(t,u,du,v) = a(du,v) +jac_t(t,u,dut,v) = ∫(dut*v)dΩ op = TransientFEOperator(res,jac,jac_t,U,V0) diff --git a/test/TransientFEsTests/HeatVectorEquationTests.jl b/test/TransientFEsTests/HeatVectorEquationTests.jl index 96475e6..e563f4a 100644 --- a/test/TransientFEsTests/HeatVectorEquationTests.jl +++ b/test/TransientFEsTests/HeatVectorEquationTests.jl @@ -42,9 +42,9 @@ a(u,v) = ∫(∇(v)⊙∇(u))dΩ m(u,v) = ∫(u⋅v)dΩ b(v,t) = ∫(v⋅f(t))dΩ -res(t,(u,ut),v) = a(u,v) + m(ut,v) - b(v,t) -jac(t,(u,ut),du,v) = a(du,v) -jac_t(t,(u,ut),dut,v) = m(dut,v) +res(t,u,v) = a(u,v) + m(∂t(u),v) - b(v,t) +jac(t,u,du,v) = a(du,v) +jac_t(t,u,dut,v) = m(dut,v) op = TransientFEOperator(res,jac,jac_t,U,V0) diff --git a/test/TransientFEsTests/StokesEquationAutoDiffTests.jl b/test/TransientFEsTests/StokesEquationAutoDiffTests.jl index cbcd8f6..750cc70 100644 --- a/test/TransientFEsTests/StokesEquationAutoDiffTests.jl +++ b/test/TransientFEsTests/StokesEquationAutoDiffTests.jl @@ -63,9 +63,9 @@ m(ut,v) = ∫(ut⋅v)dΩ X = TransientMultiFieldFESpace([U,P]) Y = MultiFieldFESpace([V0,Q]) -res(t,((u,p),(ut,pt)),(v,q)) = a(u,v) + m(ut,v) - ∫((∇⋅v)*p)dΩ + ∫(q*(∇⋅u))dΩ - b((v,q),t) -jac(t,((u,p),(ut,pt)),(du,dp),(v,q)) = a(du,v) - ∫((∇⋅v)*dp)dΩ + ∫(q*(∇⋅du))dΩ -jac_t(t,((u,p),(ut,pt)),(dut,dpt),(v,q)) = m(dut,v) +res(t,(u,p),(v,q)) = a(u,v) + m(∂t(u),v) - ∫((∇⋅v)*p)dΩ + ∫(q*(∇⋅u))dΩ - b((v,q),t) +jac(t,(u,p),(du,dp),(v,q)) = a(du,v) - ∫((∇⋅v)*dp)dΩ + ∫(q*(∇⋅du))dΩ +jac_t(t,(u,p),(dut,dpt),(v,q)) = m(dut,v) b((v,q)) = b((v,q),0.0) diff --git a/test/TransientFEsTests/StokesEquationTests.jl b/test/TransientFEsTests/StokesEquationTests.jl index d198cac..2156ab5 100644 --- a/test/TransientFEsTests/StokesEquationTests.jl +++ b/test/TransientFEsTests/StokesEquationTests.jl @@ -62,9 +62,9 @@ m(ut,v) = ∫(ut⋅v)dΩ X = TransientMultiFieldFESpace([U,P]) Y = MultiFieldFESpace([V0,Q]) -res(t,((u,p),(ut,pt)),(v,q)) = a(u,v) + m(ut,v) - ∫((∇⋅v)*p)dΩ + ∫(q*(∇⋅u))dΩ - b((v,q),t) -jac(t,((u,p),(ut,pt)),(du,dp),(v,q)) = a(du,v) - ∫((∇⋅v)*dp)dΩ + ∫(q*(∇⋅du))dΩ -jac_t(t,((u,p),(ut,pt)),(dut,dpt),(v,q)) = m(dut,v) +res(t,(u,p),(v,q)) = a(u,v) + m(∂t(u),v) - ∫((∇⋅v)*p)dΩ + ∫(q*(∇⋅u))dΩ - b((v,q),t) +jac(t,(u,p),(du,dp),(v,q)) = a(du,v) - ∫((∇⋅v)*dp)dΩ + ∫(q*(∇⋅du))dΩ +jac_t(t,(u,p),(dut,dpt),(v,q)) = m(dut,v) b((v,q)) = b((v,q),0.0) diff --git a/test/TransientFEsTests/VectorHeatEquationTests.jl b/test/TransientFEsTests/VectorHeatEquationTests.jl index b85be08..b4cee2a 100644 --- a/test/TransientFEsTests/VectorHeatEquationTests.jl +++ b/test/TransientFEsTests/VectorHeatEquationTests.jl @@ -47,11 +47,11 @@ m(ut,v) = ∫(ut⋅v)dΩ X = TransientMultiFieldFESpace([U,U]) Y = MultiFieldFESpace([V0,V0]) -_res(t,u,ut,v) = a(u,v) + m(ut,v) - b(v,t) +_res(t,u,v) = a(u,v) + m(∂t(u),v) - b(v,t) -res(t,((u1,u2),(u1t,u2t)),(v1,v2)) = _res(t,u1,u1t,v1) + _res(t,u2,u2t,v2) -jac(t,(x,xt),(du1,du2),(v1,v2)) = a(du1,v1) + a(du2,v2) -jac_t(t,(x,xt),(du1t,du2t),(v1,v2)) = m(du1t,v1) + m(du2t,v2) +res(t,(u1,u2),(v1,v2)) = _res(t,u1,v1) + _res(t,u2,v2) +jac(t,x,(du1,du2),(v1,v2)) = a(du1,v1) + a(du2,v2) +jac_t(t,x,(du1t,du2t),(v1,v2)) = m(du1t,v1) + m(du2t,v2) op = TransientFEOperator(res,jac,jac_t,X,Y) From f5600039c8955b5459137518e86cc9bcc037bfc1 Mon Sep 17 00:00:00 2001 From: Oriol Colomes Date: Thu, 1 Apr 2021 12:21:30 +0200 Subject: [PATCH 06/15] Fix in TransientCellField --- src/TransientFETools/TransientCellField.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/TransientFETools/TransientCellField.jl b/src/TransientFETools/TransientCellField.jl index aca7497..c4412ce 100644 --- a/src/TransientFETools/TransientCellField.jl +++ b/src/TransientFETools/TransientCellField.jl @@ -14,8 +14,8 @@ change_domain(f::TransientCellField,trian::Triangulation,target_domain::DomainSt # MultiFieldCellField methods num_fields(f::TransientCellField{A}) where A = num_fields(f.cellfield) -function Base.getindex(f::TransientCellField{A},i::Integer) where A - singleCellfield = getindex(f.single_cellfields,i) +function Base.getindex(f::TransientCellField{A},i::Integer) where A + singleCellfield = getindex(f.cellfield,i) singleDerivatives = (getindex(i_derivatives,i) for i_derivatives in f.derivatives) TransientCellField(singleCellfield,singleDerivatives) end @@ -62,4 +62,4 @@ function ∂t(f::Union{TransientCellField,TransientFEBasis}) TransientCellField(cellfield,derivatives) end -∂tt(f::Union{TransientCellField,TransientFEBasis}) = ∂t(∂t(f::Union{TransientCellField,TransientFEBasis})) \ No newline at end of file +∂tt(f::Union{TransientCellField,TransientFEBasis}) = ∂t(∂t(f::Union{TransientCellField,TransientFEBasis})) From 305f70d5aa2af9e6a3c9a131684ecee750233ddb Mon Sep 17 00:00:00 2001 From: Oriol Colomes Date: Sat, 3 Apr 2021 21:05:30 +0200 Subject: [PATCH 07/15] Fixes for autodiff in new API --- src/TransientFETools/TransientCellField.jl | 10 +++++++--- src/TransientFETools/TransientFEOperators.jl | 17 ++++++++++++----- 2 files changed, 19 insertions(+), 8 deletions(-) diff --git a/src/TransientFETools/TransientCellField.jl b/src/TransientFETools/TransientCellField.jl index aca7497..3b31885 100644 --- a/src/TransientFETools/TransientCellField.jl +++ b/src/TransientFETools/TransientCellField.jl @@ -1,7 +1,7 @@ # Transient CellField struct TransientCellField{A} <: CellField cellfield::A - derivatives::Tuple{Vararg{A}} + derivatives::Tuple#{Vararg{A,B} where B} end # CellField methods @@ -15,8 +15,12 @@ change_domain(f::TransientCellField,trian::Triangulation,target_domain::DomainSt # MultiFieldCellField methods num_fields(f::TransientCellField{A}) where A = num_fields(f.cellfield) function Base.getindex(f::TransientCellField{A},i::Integer) where A - singleCellfield = getindex(f.single_cellfields,i) - singleDerivatives = (getindex(i_derivatives,i) for i_derivatives in f.derivatives) + singleCellfield = getindex(f.cellfield,i) + singleDerivatives = () + for i_derivatives in f.derivatives + singleDerivatives = (singleDerivatives...,getindex(i_derivatives,i)) + end + # singleDerivatives = (getindex(i_derivatives,i) for i_derivatives in f.derivatives) TransientCellField(singleCellfield,singleDerivatives) end function Base.iterate(f::TransientCellField{A}) where A diff --git a/src/TransientFETools/TransientFEOperators.jl b/src/TransientFETools/TransientFEOperators.jl index 38dd440..02045bb 100644 --- a/src/TransientFETools/TransientFEOperators.jl +++ b/src/TransientFETools/TransientFEOperators.jl @@ -171,20 +171,27 @@ function TransientFEOperator(res::Function,jac::Function,jac_t::Function, end function TransientFEOperator(res::Function,trial,test;order::Integer=1) - jacs = () - for i in 1:order+1 + function jac_0(t,x,dx0,dv) + function res_0(y) + x0 = TransientCellField(y,x.derivatives) + res(t,x0,dv) + end + jacobian(res_0,x.cellfield) + end + jacs = (jac_0,) + for i in 1:order function jac_i(t,x,dxi,dv) function res_i(y) - x[i] = y + derivatives = (x.derivatives[1:i-1],y,x.derivatives[i+1:end]) res(t,x,dv) end - jacobian(res_i,x[i]) + jacobian(res_i,x.derivatives[i]) end jacs = (jacs...,jac_i) end # jac(t,u,ut,du,dv) = jacobian(x->res(t,x,ut,dv),u) # jac_t(t,u,ut,dut,dv) = jacobian(xt->res(t,u,xt,dv),ut) - TransientFEOperator(res,jacs,trial,test) + TransientFEOperator(res,jacs...,trial,test) end function SparseMatrixAssembler( From 9e82d54a5c9837e32cc039ee9d7779e510ae2fe3 Mon Sep 17 00:00:00 2001 From: Oriol Colomes Date: Sat, 3 Apr 2021 21:05:54 +0200 Subject: [PATCH 08/15] Added test for mixed-dims (not fully working yet) --- .../FreeSurfacePotentialFlowTest.jl | 112 ++++++++++++++++++ 1 file changed, 112 insertions(+) create mode 100644 test/TransientFEsTests/FreeSurfacePotentialFlowTest.jl diff --git a/test/TransientFEsTests/FreeSurfacePotentialFlowTest.jl b/test/TransientFEsTests/FreeSurfacePotentialFlowTest.jl new file mode 100644 index 0000000..04f150b --- /dev/null +++ b/test/TransientFEsTests/FreeSurfacePotentialFlowTest.jl @@ -0,0 +1,112 @@ +module FreeSurfacePotentialFlowTests + +using Gridap +using Gridap.Geometry +using GridapODEs.TransientFETools +using GridapODEs.ODETools + +# Parameters +L = 2*π +H = 1.0 +n = 4 +order = 2 +g = 9.81 +ξ = 0.1 +λ = L/2 +k = 2*π/L +h = L/n +ω = √(g*k*tanh(k*H)) +t₀ = 0.0 +tf = 8*π +Δt = h/(2*λ*ω) +θ = 0.5 + +# Exact solution +ϕₑ(x,t) = ω/k * ξ * (cosh(k*(x[2]))) / sinh(k*H) * sin(k*x[1] - ω*t) +ηₑ(x,t) = ξ * cos(k*x[1] - ω*t) +ϕₑ(t::Real) = x -> ϕₑ(x,t) +ηₑ(t::Real) = x -> ηₑ(x,t) + +# Domain +domain = (0,L,0,H) +partition = (n,n) +model = CartesianDiscreteModel(domain,partition;isperiodic=(true,false)) + +# Boundaries +labels = get_face_labeling(model) +add_tag_from_tags!(labels,"bottom",[1,2,5]) +add_tag_from_tags!(labels,"free_surface",[3,4,6]) +bgface_to_mask = get_face_mask(labels,[3,4,6],1) +model_Γ =BoundaryDiscreteModel(Polytope{1},model,bgface_to_mask) + +# Triangulation +Ω = Triangulation(model) +Γ = Triangulation(model_Γ) +dΩ = Measure(Ω,2*order) +dΓ = Measure(Γ,2*order) + +# FE spaces +reffe = ReferenceFE(lagrangian,Float64,order) +V = TestFESpace(model,reffe,conformity=:H1) +V_Γ = TestFESpace(model_Γ,reffe,conformity=:H1) +U = TransientTrialFESpace(V) +U_Γ = TransientTrialFESpace(V_Γ) +X = TransientMultiFieldFESpace([U,U_Γ]) +Y = MultiFieldFESpace([V,V_Γ]) + +# Weak form +α = 2/Δt + +# Optimal transient FE Operator: +m((ϕt,ηt),(w,v)) = ∫( 0.5*(α/g*(w*ϕt) + v*ϕt) - (w*ηt) )dΓ +a((ϕ,η),(w,v)) = ∫( ∇(ϕ)⋅∇(w) )dΩ + ∫( 0.5*(α*(w*η) + g*v*η) )dΓ +b((w,v)) = ∫( 0.0*w )dΓ +# op = TransientConstantFEOperator(m,a,b,X,Y) + +# TransientFEOperator exploiting automatic differentiation (testing purposes) +res(t,x,y) = m(∂t(x),y) + a(x,y) - b(y) +jac(t,x,dx,y) = a(dx,y) +jac_t(t,x,dxt,y) = m(dxt,y) +# res(t,(ϕ,η),(w,v)) = m(∂t((ϕ,η)),(w,v)) + a((ϕ,η),(w,v)) - b((w,v))#∫( ∇(ϕ)⋅∇(w) )dΩ + ∫( -∂t(η)*w + (∂t(ϕ)+g*η) * 0.5*(v + α/g*w) )dΓ +# jac(t,(ϕ,η),(dϕ,dη),(w,v)) = a((dϕ,dη),(w,v))#∫( ∇(dϕ)⋅∇(w) )dΩ + ∫( g*dη * 0.5*(v + α/g*w) )dΓ +# jac_t(t,(ϕ,η),(dϕt,dηt),(w,v)) = m((dϕt,dηt),(w,v))#∫( -dηt*w + dϕt * 0.5*(v + α/g*w) )dΓ +#op = TransientFEOperator(res,jac,jac_t,X,Y) +op = TransientFEOperator(res,X,Y) + +# Solver +ls = LUSolver() +odes = ThetaMethod(ls,Δt,θ) +solver = TransientFESolver(odes) + +# Initial solution +x₀ = interpolate_everywhere([ϕₑ(0.0),ηₑ(0.0)],X(0.0)) + +# Solution +sol_t = solve(solver,op,x₀,t₀,tf) + +# Post-process +l2_Ω(w) = √(∑(∫(w*w)dΩ)) +l2_Γ(v) = √(∑(∫(v*v)dΓ)) +E_kin(w) = 0.5*∑( ∫(∇(w)⋅∇(w))dΩ ) +E_pot(v) = g*0.5*∑( ∫(v*v)dΓ ) +Eₑ = 0.5*g*ξ^2*L + +# folderName = "ϕFlow-results" +# if !isdir(folderName) +# mkdir(folderName) +# end +# filePath_Ω = folderName * "/fields_Ω" +# filePath_Γ = folderName * "/fields_Γ" +# pvd_Ω = paraview_collection(filePath_Ω, append=false) +# pvd_Γ = paraview_collection(filePath_Γ, append=false) +for ((ϕn,ηn),tn) in sol_t + E = E_kin(ϕn) + E_pot(ηn) + error_ϕ = l2_Ω(ϕn-ϕₑ(tn)) + error_η = l2_Γ(ηn-ηₑ(tn)) + println(E/Eₑ," ", error_ϕ," ",error_η) + + # pvd_Ω[tn] = createvtk(Ω,filePath_Ω * "_$tn.vtu",cellfields = ["phi" => ϕn]) + # pvd_Γ[tn] = createvtk(Γ,filePath_Γ * "_$tn.vtu",cellfields = ["eta" => ηn]) +end + +end \ No newline at end of file From f0b533fe587f5c2e185212ef69ea235d3de8b8cb Mon Sep 17 00:00:00 2001 From: Oriol Colomes Date: Mon, 12 Apr 2021 13:01:08 +0200 Subject: [PATCH 09/15] Added FreeSurfacePotentialFlowTests.jl (test transient problem for mixed-dimensional PDEs) --- ...st.jl => FreeSurfacePotentialFlowTests.jl} | 29 ++++++------------- test/TransientFEsTests/runtests.jl | 2 ++ 2 files changed, 11 insertions(+), 20 deletions(-) rename test/TransientFEsTests/{FreeSurfacePotentialFlowTest.jl => FreeSurfacePotentialFlowTests.jl} (69%) diff --git a/test/TransientFEsTests/FreeSurfacePotentialFlowTest.jl b/test/TransientFEsTests/FreeSurfacePotentialFlowTests.jl similarity index 69% rename from test/TransientFEsTests/FreeSurfacePotentialFlowTest.jl rename to test/TransientFEsTests/FreeSurfacePotentialFlowTests.jl index 04f150b..649be12 100644 --- a/test/TransientFEsTests/FreeSurfacePotentialFlowTest.jl +++ b/test/TransientFEsTests/FreeSurfacePotentialFlowTests.jl @@ -4,11 +4,12 @@ using Gridap using Gridap.Geometry using GridapODEs.TransientFETools using GridapODEs.ODETools +using Test # Parameters L = 2*π H = 1.0 -n = 4 +n = 8 order = 2 g = 9.81 ξ = 0.1 @@ -17,7 +18,7 @@ k = 2*π/L h = L/n ω = √(g*k*tanh(k*H)) t₀ = 0.0 -tf = 8*π +tf = 2*π Δt = h/(2*λ*ω) θ = 0.5 @@ -67,11 +68,7 @@ b((w,v)) = ∫( 0.0*w )dΓ res(t,x,y) = m(∂t(x),y) + a(x,y) - b(y) jac(t,x,dx,y) = a(dx,y) jac_t(t,x,dxt,y) = m(dxt,y) -# res(t,(ϕ,η),(w,v)) = m(∂t((ϕ,η)),(w,v)) + a((ϕ,η),(w,v)) - b((w,v))#∫( ∇(ϕ)⋅∇(w) )dΩ + ∫( -∂t(η)*w + (∂t(ϕ)+g*η) * 0.5*(v + α/g*w) )dΓ -# jac(t,(ϕ,η),(dϕ,dη),(w,v)) = a((dϕ,dη),(w,v))#∫( ∇(dϕ)⋅∇(w) )dΩ + ∫( g*dη * 0.5*(v + α/g*w) )dΓ -# jac_t(t,(ϕ,η),(dϕt,dηt),(w,v)) = m((dϕt,dηt),(w,v))#∫( -dηt*w + dϕt * 0.5*(v + α/g*w) )dΓ -#op = TransientFEOperator(res,jac,jac_t,X,Y) -op = TransientFEOperator(res,X,Y) +op = TransientFEOperator(res,jac,jac_t,X,Y) # Solver ls = LUSolver() @@ -91,22 +88,14 @@ E_kin(w) = 0.5*∑( ∫(∇(w)⋅∇(w))dΩ ) E_pot(v) = g*0.5*∑( ∫(v*v)dΓ ) Eₑ = 0.5*g*ξ^2*L -# folderName = "ϕFlow-results" -# if !isdir(folderName) -# mkdir(folderName) -# end -# filePath_Ω = folderName * "/fields_Ω" -# filePath_Γ = folderName * "/fields_Γ" -# pvd_Ω = paraview_collection(filePath_Ω, append=false) -# pvd_Γ = paraview_collection(filePath_Γ, append=false) +tol = 1.0e-2 for ((ϕn,ηn),tn) in sol_t E = E_kin(ϕn) + E_pot(ηn) error_ϕ = l2_Ω(ϕn-ϕₑ(tn)) error_η = l2_Γ(ηn-ηₑ(tn)) - println(E/Eₑ," ", error_ϕ," ",error_η) - - # pvd_Ω[tn] = createvtk(Ω,filePath_Ω * "_$tn.vtu",cellfields = ["phi" => ϕn]) - # pvd_Γ[tn] = createvtk(Γ,filePath_Γ * "_$tn.vtu",cellfields = ["eta" => ηn]) + @test abs(E/Eₑ-1.0) <= tol + @test error_ϕ <= tol + @test error_η <= tol end -end \ No newline at end of file +end diff --git a/test/TransientFEsTests/runtests.jl b/test/TransientFEsTests/runtests.jl index 95a0ca2..c8a55be 100644 --- a/test/TransientFEsTests/runtests.jl +++ b/test/TransientFEsTests/runtests.jl @@ -24,4 +24,6 @@ using Test @testset "DGHeatEquationTests" begin include("DGHeatEquationTests.jl") end +@testset "FreeSurfacePotentialFlowTests" begin include("FreeSurfacePotentialFlowTests.jl") end + end # module From 624cad4ca5f14b3e3097ceee180d3b5ba4d30124 Mon Sep 17 00:00:00 2001 From: Oriol Colomes Date: Mon, 12 Apr 2021 13:07:10 +0200 Subject: [PATCH 10/15] Fixed DiffEqsWrappersTests --- test/DiffEqsWrappersTests/DiffEqsTests.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/DiffEqsWrappersTests/DiffEqsTests.jl b/test/DiffEqsWrappersTests/DiffEqsTests.jl index 59ba5a1..149d58a 100644 --- a/test/DiffEqsWrappersTests/DiffEqsTests.jl +++ b/test/DiffEqsWrappersTests/DiffEqsTests.jl @@ -40,9 +40,9 @@ function fe_problem(u, n) b(v, t) = ∫( v * f(t) )dΩ m(u, v) = ∫( v * u )dΩ - res(t, (u, ut), v) = a(u, v) + m(ut, v) - b(v, t) - jac(t, (u, ut), du, v) = a(du, v) - jac_t(t, (u, ut), dut, v) = m(dut, v) + res(t, u, v) = a(u, v) + m(∂t(u), v) - b(v, t) + jac(t, u, du, v) = a(du, v) + jac_t(t, u, dut, v) = m(dut, v) op = TransientFEOperator(res, jac, jac_t, U, V0) From 47568999b138e7dda95e2b2897a7fafc99c304ae Mon Sep 17 00:00:00 2001 From: Oriol Colomes Date: Tue, 20 Apr 2021 17:30:47 +0200 Subject: [PATCH 11/15] debugging tmp file --- tmp.jl | 305 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 305 insertions(+) create mode 100644 tmp.jl diff --git a/tmp.jl b/tmp.jl new file mode 100644 index 0000000..988da7d --- /dev/null +++ b/tmp.jl @@ -0,0 +1,305 @@ +using Gridap +using Gridap.Arrays +using Gridap.Geometry +using GridapODEs.TransientFETools +using GridapODEs.ODETools +using Test +using LineSearches: BackTracking +using ForwardDiff +import GridapODEs.TransientFETools: ∂t + +# Problem setting +println("Setting analytical fluid problem parameters") +# Solid properties +E_s = 1.0 +ν_s = 0.4 +ρ_s = 1.0 +# Fluid properties +ρ_f = 1.0 +μ_f = 1.0 +γ_f = 1.0 +# Mesh properties +n_m = 3 +E_m = 1.0 +ν_m = 0.1 +α_m = 1.0e-5 +# Time stepping +t0 = 0.0 +tf = 0.5 +dt = 0.1 +θ = 0.5 + +# Define BC functions +println("Defining Boundary conditions") +#u(x, t) = 1.0e-2 * VectorValue( x[1]^2*x[2] , -x[1]*x[2]^2) * t +#v(x, t) = 1.0e-2 * VectorValue( x[1]^2*x[2] , -x[1]*x[2]^2) +u(x, t) = 1.0e-2 * VectorValue( 1.0 , -1.0) * t +v(x, t) = 1.0e-2 * VectorValue( 1.0 , -1.0) +u(t::Real) = x -> u(x,t) +v(t::Real) = x -> v(x,t) +∂tu(t) = x -> VectorValue(ForwardDiff.derivative(t -> get_array(u(x,t)),t)) +∂tv(t) = x -> VectorValue(ForwardDiff.derivative(t -> get_array(v(x,t)),t)) +∂tu(x,t) = ∂tu(t)(x) +∂tv(x,t) = ∂tv(t)(x) +T_u = typeof(u) +T_v = typeof(v) +@eval ∂t(::$T_u) = $∂tu +@eval ∂t(::$T_v) = $∂tv +p(x,t) = 0.0 +p(t::Real) = x -> p(x,t) +bconds = ( + # Tags + FSI_Vu_tags = ["boundary"], + FSI_Vv_tags = ["boundary"], + ST_Vu_tags = ["boundary","interface"], + ST_Vv_tags = ["boundary","interface"], + # Values, + FSI_Vu_values = [u], + FSI_Vv_values = [v], + ST_Vu_values = [u(0.0),u(0.0)], + ST_Vv_values = [v(0.0),v(0.0)], +) + +function lame_parameters(E,ν) + λ = (E*ν)/((1+ν)*(1-2*ν)) + μ = E/(2*(1+ν)) + (λ, μ) +end + +# Define Forcing terms +I = TensorValue( 1.0, 0.0, 0.0, 1.0 ) +F(t) = x -> ∇(u(t))(x) + I +J(t) = x -> det(F(t))(x) +E(t) = x -> 0.5 * ((F(t)(x)')⋅F(t)(x) - I) +(λ_s, μ_s) = lame_parameters(E_s,ν_s) +(λ_m, μ_m) = lame_parameters(E_m,ν_m) +S_SV(t) = x -> 2*μ_s*E(t)(x) + λ_s*tr(E(t)(x))*I +fv_ST_Ωf(t) = x -> - μ_f*Δ(v(t))(x) + ∇(p(t))(x) +fu_Ωf(t) = x -> - α_m * Δ(u(t))(x) +fv_Ωf(t) = x -> ρ_f * ∂t(v)(t)(x) - μ_f * Δ(v(t))(x) + ∇(p(t))(x) + ρ_f*( (∇(v(t))(x)')⋅(v(t)(x) - ∂t(u)(t)(x)) ) +fp_Ωf(t) = x -> (∇⋅v(t))(x) +fu_Ωs(t) = x -> ∂t(u)(t)(x) - v(t)(x) +fv_Ωs(t) = x -> ρ_s * ∂t(v)(t)(x) #- (∇⋅(F(t)⋅S_SV(t)))(x) # Divergence of a a doted function not working yet... + +# Discrete model +println("Defining discrete model") +domain = (-1,1,-1,1) +partition = (n_m,n_m) +model = CartesianDiscreteModel(domain,partition) +trian = Triangulation(model) +R = 0.5 +xc = 0.0 +yc = 0.0 +function is_in(coords) + n = length(coords) + x = (1/n)*sum(coords) + d = (x[1]-xc)^2 + (x[2]-yc)^2 - R^2 + d < 1.0e-8 +end +oldcell_to_coods = get_cell_coordinates(trian) +oldcell_to_is_in = collect1d(lazy_map(is_in,oldcell_to_coods)) +incell_to_cell = findall(oldcell_to_is_in) +outcell_to_cell = findall(collect(Bool, .! oldcell_to_is_in)) +model_solid = DiscreteModel(model,incell_to_cell) +model_fluid = DiscreteModel(model,outcell_to_cell) + +# Build fluid-solid interface labelling +println("Defining Fluid-Solid interface") +labeling = get_face_labeling(model_fluid) +new_entity = num_entities(labeling) + 1 +topo = get_grid_topology(model_fluid) +D = num_cell_dims(model_fluid) +for d in 0:D-1 + fluid_boundary_mask = collect(Bool,get_isboundary_face(topo,d)) + fluid_outer_boundary_mask = get_face_mask(labeling,"boundary",d) + fluid_interface_mask = collect(Bool,fluid_boundary_mask .!= fluid_outer_boundary_mask) + dface_list = findall(fluid_interface_mask) + for face in dface_list + labeling.d_to_dface_to_entity[d+1][face] = new_entity + end +end +add_tag!(labeling,"interface",[new_entity]) + +# Triangulations +println("Defining triangulations") +Ω = Triangulation(model) +Ω_s = Triangulation(model_solid) +Ω_f = Triangulation(model_fluid) +Γ_i = BoundaryTriangulation(model_fluid,tags="interface") + +# Quadratures +println("Defining quadratures") +order = 2 +degree = 2*order +bdegree = 2*order +dΩ = Measure(Ω,degree) +dΩs = Measure(Ω_s,degree) +dΩf = Measure(Ω_f,degree) +dΓi = Measure(Γ_i,bdegree) + +# Test FE Spaces +println("Defining FE spaces") +# Reference FE +reffeᵤ = ReferenceFE(lagrangian,VectorValue{2,Float64},order) +reffeₚ = ReferenceFE(lagrangian,Float64,order-1) +# Test FE Spaces +Vu_FSI = TestFESpace(model, reffeᵤ, conformity=:H1, dirichlet_tags=bconds[:FSI_Vu_tags]) +Vv_FSI = TestFESpace(model, reffeᵤ, conformity=:H1, dirichlet_tags=bconds[:FSI_Vv_tags]) +Vu_ST = TestFESpace(model_fluid, reffeᵤ, conformity=:H1, dirichlet_tags=bconds[:ST_Vu_tags]) +Vv_ST = TestFESpace(model_fluid, reffeᵤ, conformity=:H1, dirichlet_tags=bconds[:ST_Vv_tags]) +Q = TestFESpace(model_fluid, reffeₚ, conformity=:C0, constraint=:zeromean) + +# Trial FE Spaces +Uu_ST = TrialFESpace(Vu_ST,bconds[:ST_Vu_values]) +Uv_ST = TrialFESpace(Vv_ST,bconds[:ST_Vv_values]) +Uu_FSI = TransientTrialFESpace(Vu_FSI,bconds[:FSI_Vu_values]) +Uv_FSI = TransientTrialFESpace(Vv_FSI,bconds[:FSI_Vv_values]) +P = TrialFESpace(Q) + +# Multifield FE Spaces +Y_ST = MultiFieldFESpace([Vu_ST,Vv_ST,Q]) +X_ST = MultiFieldFESpace([Uu_ST,Uv_ST,P]) +Y_FSI = MultiFieldFESpace([Vu_FSI,Vv_FSI,Q]) +X_FSI = TransientMultiFieldFESpace([Uu_FSI,Uv_FSI,P]) + +# Stokes problem for initial solution +println("Defining Stokes operator") +a((u,v,p),(ϕ,φ,q)) = ∫( ϕ⋅u + ε(φ) ⊙ (2.0*μ_f*ε(v)) - (∇⋅φ) * p + q * (∇⋅v) )dΩf +l((ϕ,φ,q)) = ∫( φ⋅fv_ST_Ωf(0.0) )dΩf +res(x,y) = a(x,y) - l(y) +jac(x,dx,y) = a(dx,y) +op_ST = FEOperator(res,jac,X_ST,Y_ST) + +# Map operations +F(∇u) = ∇u + I +J(∇u) = det(F(∇u)) +Finv(∇u) = inv(F(∇u)) +FinvT(∇u) = (Finv(∇u)') +E(∇u) = 0.5 * ((F(∇u)')⋅F(∇u) - I) + +# Map operations derivatives +dF(∇du) = ∇du +dJ(∇u,∇du) = J(∇u)*tr(Finv(∇u)⋅dF(∇du)) +dFinv(∇u,∇du) = -Finv(∇u) ⋅ dF(∇du) ⋅ Finv(∇u) +dFinvT(∇u,∇du) = (dFinv(∇u,∇du)') +dE(∇u,∇du) = 0.5 * ((dF(∇du)')⋅F(∇u) + (F(∇u)')⋅dF(∇du)) + +# Fluid constitutive laws +# Cauchy stress +σᵥ_Ωf(μ,u) = 2.0*μ*ε(u) +σᵥ_Ωf(μ) = (∇v,Finv) -> μ*(∇v⋅Finv + (Finv')⋅(∇v')) +# First Piola-Kirchhoff stress +Pᵥ_Ωf(μ,u,v) = (J∘∇(u)) * (σᵥ_Ωf(μ)∘(∇(v),Finv∘∇(u)))' ⋅ (FinvT∘∇(u)) +function Pₚ_Ωf(u,p) + typeof(- (J∘∇(u)) * p * tr((FinvT∘∇(u)))) + - (J∘∇(u)) * p * tr((FinvT∘∇(u))) +end +# First Piola-Kirchhoff stress Jacobian +dPᵥ_Ωf_du(μ,u,du,v) = (dJ∘(∇(u),∇(du))) * (σᵥ_Ωf(μ)∘(∇(v),(Finv∘∇(u))))' ⋅ (FinvT∘∇(u)) + + (J∘∇(u)) * (σᵥ_Ωf(μ)∘(∇(v),(dFinv∘(∇(u),∇(du)))))' ⋅ (FinvT∘∇(u)) + + (J∘∇(u)) * (σᵥ_Ωf(μ)∘(∇(v),(Finv∘∇(u))))' ⋅ (dFinvT∘(∇(u),∇(du))) +dPₚ_Ωf_du(u,du,p) = - p * ( (dJ∘(∇(u),∇(du))) * tr((FinvT∘∇(u))) + + (J∘∇(u)) * tr((dFinvT∘(∇(u),∇(du)))) ) +dPᵥ_Ωf_dv(μ,u,dv) = Pᵥ_Ωf(μ,u,dv) +dPₚ_Ωf_dp(u,dp) = Pₚ_Ωf(u,dp) +# Convective term +conv(c,∇v) = (∇v') ⋅ c +dconv(dc,∇dv,c,∇v) = conv(c,∇dv) + conv(dc,∇v) + +# Solid constitutive laws +# Second Piola-Kirchhoff stress (Saint-Venant) +Sₛᵥ_Ωs(λ,μ,u) = 2*μ*(E∘∇(u)) + λ*tr((E∘∇(u)))*(I∘∇(u)) +dSₛᵥ_Ωs_du(λ,μ,u,du) = 2*μ*(dE∘(∇(u),∇(du))) + λ*tr((dE∘(∇(u),∇(du))))*(I∘∇(u)) +# First Piola-Kirchhoff stress (Saint-Venant) +Pₛᵥ_Ωs(λ,μ,u) = (F∘∇(u)) ⋅ Sₛᵥ_Ωs(λ,μ,u) +dPₛᵥ_Ωs_du(λ,μ,u,du) = (dF∘∇(du)) ⋅ Sₛᵥ_Ωs(λ,μ,u) + (F∘∇(u)) ⋅ dSₛᵥ_Ωs_du(λ,μ,u,du) + +# # FSI problem +println("Defining FSI operator") +n_Γi = get_normal_vector(Γ_i) +a_f((u,v,p),(ϕ,φ,q)) = ∫( φ ⋅ ( (J∘∇(u)) * ρ_f * ∂t(v) ) )dΩ +da_f((u,v,p),(du,dv,dp),(ϕ,φ,q)) = ∫( φ ⋅ ( (dJ∘(∇(u),∇(du))) * ρ_f * ∂t(v) ) )dΩ +da_t_f((u,v,p),(dut,dvt,dpt),(ϕ,φ,q)) = ∫( φ ⋅ ( (J∘∇(u)) * ρ_f * dvt ) )dΩ +# da_t_f((u,v,p),(dut,dvt,dpt),(ϕ,φ,q)) = ∫( φ ⋅ ( ρ_f * dvt) )dΩf +a_fsi((u,v,p),(ϕ,φ,q)) = ∫( ϕ⋅u + φ⋅v + q*p )dΩ +l_f(t,(ϕ,φ,q)) = ∫( ϕ⋅fu_Ωf(t) )dΩ +l_fsi(t,(ϕ,φ,q)) = ∫( φ⋅fv_Ωf(t) )dΩ +res(t,x,y) = a_f(x,y) + a_fsi(x,y) + a_fsi(∂t(x),y) - l_f(t,y) - l_fsi(t,y) +jac(t,x,dx,y) = da_f(x,dx,y) + a_fsi(dx,y) +jac_t(t,x,dxt,y) = da_t_f(x,dxt,y) + a_fsi(dxt,y) +op_FSI = TransientFEOperator(res,jac,jac_t,X_FSI,Y_FSI) +#op_FSI = TransientFEOperator(res,X_FSI,Y_FSI) +# a_f((u,v,p),(ut,vt,pt),(ϕ,φ,q)) = +# ∫( α_m * (∇(ϕ) ⊙ ∇(u)) )dΩf + +# ∫( φ ⋅ ( (J∘∇(u)) * ρ_f * vt ) )dΩf +# da_f((u,v,p),(ut,vt,pt),(du,dv,dp),(ϕ,φ,q)) = +# ∫( α_m * (∇(ϕ) ⊙ ∇(du)) )dΩf + +# ∫( φ ⋅ ( (dJ∘(∇(u),∇(du))) * ρ_f * vt ) )dΩf +# da_t_f((u,v,p),(dut,dvt,dpt),(ϕ,φ,q)) = ∫( φ ⋅ ( (J∘∇(u)) * ρ_f * dvt ) )dΩf +# a_fsi((u,v,p),(ϕ,φ,q)) = ∫( ϕ⋅u + φ⋅v + q*p )dΩf + ∫( ϕ⋅u + φ⋅v + q*p )dΩs +# l_f(t,(ϕ,φ,q)) = ∫( ϕ⋅fu_Ωf(t) )dΩf +# l_fsi(t,(ϕ,φ,q)) = ∫( φ⋅fv_Ωf(t) )dΩf +# res(t,(x,xt),y) = a_f(x,xt,y) + a_fsi(x,y) + a_fsi(xt,y) - l_f(t,y) - l_fsi(t,y) +# jac(t,(x,xt),dx,y) = da_f(x,xt,dx,y) + a_fsi(dx,y) +# jac_t(t,(x,xt),dxt,y) = da_t_f(x,dxt,y) + a_fsi(dxt,y) +# op_FSI = TransientFEOperator(res,jac,jac_t,X_FSI,Y_FSI) + +# Setup output files +folderName = "fsi-results" +fileName = "fields" +if !isdir(folderName) + mkdir(folderName) +end +filePath = join([folderName, fileName], "/") + +# Solve Stokes problem +println("Defining Stokes solver") +xh = solve(op_ST) +writevtk(Ω_f, filePath, cellfields = ["uh"=>xh[1], "vh"=>xh[2], "ph"=>xh[3]]) + +# Compute Stokes solution L2-norm +l2(w) = w⋅w +eu_ST = u(0.0) - xh[1] +ev_ST = v(0.0) - xh[2] +ep_ST = p(0.0) - xh[3] +eul2_ST = sqrt(∑( ∫(l2(eu_ST))dΩf )) +evl2_ST = sqrt(∑( ∫(l2(ev_ST))dΩf )) +epl2_ST = sqrt(∑( ∫(l2(ep_ST))dΩf )) +println("Stokes L2-norm u: ", eul2_ST) +println("Stokes L2-norm v: ", evl2_ST) +println("Stokes L2-norm p: ", epl2_ST) +@test eul2_ST < 1.0e-10 +@test evl2_ST < 1.0e-10 +@test epl2_ST < 1.0e-10 + +# Solve FSI problem +println("Defining FSI solver") +xh0 = interpolate_everywhere([u(0.0),v(0.0),p(0.0)],X_FSI(0.0)) +nls = NLSolver( +show_trace = true, +method = :newton, +linesearch = BackTracking(), +ftol = 1.0e-10, +iterations = 50 +) +odes = ThetaMethod(nls, dt, 0.5) +solver = TransientFESolver(odes) +xht = solve(solver, op_FSI, xh0, t0, tf) + +for (i,(xh, t)) in enumerate(xht) + println("STEP: $i, TIME: $t") + println("============================") + + # Compute errors + eu = u(t) - xh[1] + ev = v(t) - xh[2] + ep = p(t) - xh[3] + eul2 = sqrt(∑( ∫(l2(eu))dΩ )) + evl2 = sqrt(∑( ∫(l2(ev))dΩ )) + epl2 = sqrt(∑( ∫(l2(ep))dΩ )) + + # Test checks + println("FSI L2-norm u: ", eul2) + println("FSI L2-norm v: ", evl2) + println("FSI L2-norm p: ", epl2) +end \ No newline at end of file From 03b65fb8ea9f06933ec408fa73d8e019ae18e446 Mon Sep 17 00:00:00 2001 From: Oriol Colomes Date: Thu, 9 Dec 2021 21:48:34 +0100 Subject: [PATCH 12/15] Adding Free-surface with potential flow test (mixed-dimensional PDE) --- .../FreeSurfacePotentialFlowTests.jl | 61 ++++++++++--------- 1 file changed, 33 insertions(+), 28 deletions(-) diff --git a/test/TransientFEsTests/FreeSurfacePotentialFlowTests.jl b/test/TransientFEsTests/FreeSurfacePotentialFlowTests.jl index 649be12..b222544 100644 --- a/test/TransientFEsTests/FreeSurfacePotentialFlowTests.jl +++ b/test/TransientFEsTests/FreeSurfacePotentialFlowTests.jl @@ -37,19 +37,17 @@ model = CartesianDiscreteModel(domain,partition;isperiodic=(true,false)) labels = get_face_labeling(model) add_tag_from_tags!(labels,"bottom",[1,2,5]) add_tag_from_tags!(labels,"free_surface",[3,4,6]) -bgface_to_mask = get_face_mask(labels,[3,4,6],1) -model_Γ =BoundaryDiscreteModel(Polytope{1},model,bgface_to_mask) # Triangulation -Ω = Triangulation(model) -Γ = Triangulation(model_Γ) +Ω = Interior(model) +Γ = Boundary(model,tags="free_surface") dΩ = Measure(Ω,2*order) dΓ = Measure(Γ,2*order) # FE spaces reffe = ReferenceFE(lagrangian,Float64,order) -V = TestFESpace(model,reffe,conformity=:H1) -V_Γ = TestFESpace(model_Γ,reffe,conformity=:H1) +V = TestFESpace(Ω,reffe,conformity=:H1) +V_Γ = TestFESpace(Γ,reffe,conformity=:H1) U = TransientTrialFESpace(V) U_Γ = TransientTrialFESpace(V_Γ) X = TransientMultiFieldFESpace([U,U_Γ]) @@ -62,40 +60,47 @@ Y = MultiFieldFESpace([V,V_Γ]) m((ϕt,ηt),(w,v)) = ∫( 0.5*(α/g*(w*ϕt) + v*ϕt) - (w*ηt) )dΓ a((ϕ,η),(w,v)) = ∫( ∇(ϕ)⋅∇(w) )dΩ + ∫( 0.5*(α*(w*η) + g*v*η) )dΓ b((w,v)) = ∫( 0.0*w )dΓ -# op = TransientConstantFEOperator(m,a,b,X,Y) +op_const = TransientConstantFEOperator(m,a,b,X,Y) # TransientFEOperator exploiting automatic differentiation (testing purposes) res(t,x,y) = m(∂t(x),y) + a(x,y) - b(y) jac(t,x,dx,y) = a(dx,y) jac_t(t,x,dxt,y) = m(dxt,y) -op = TransientFEOperator(res,jac,jac_t,X,Y) +op_trans = TransientFEOperator(res,jac,jac_t,X,Y) +op_ad = TransientFEOperator(res,X,Y) # Solver ls = LUSolver() -odes = ThetaMethod(ls,Δt,θ) -solver = TransientFESolver(odes) +ode_solver = ThetaMethod(ls,Δt,θ) # Initial solution x₀ = interpolate_everywhere([ϕₑ(0.0),ηₑ(0.0)],X(0.0)) -# Solution -sol_t = solve(solver,op,x₀,t₀,tf) - -# Post-process -l2_Ω(w) = √(∑(∫(w*w)dΩ)) -l2_Γ(v) = √(∑(∫(v*v)dΓ)) -E_kin(w) = 0.5*∑( ∫(∇(w)⋅∇(w))dΩ ) -E_pot(v) = g*0.5*∑( ∫(v*v)dΓ ) -Eₑ = 0.5*g*ξ^2*L - -tol = 1.0e-2 -for ((ϕn,ηn),tn) in sol_t - E = E_kin(ϕn) + E_pot(ηn) - error_ϕ = l2_Ω(ϕn-ϕₑ(tn)) - error_η = l2_Γ(ηn-ηₑ(tn)) - @test abs(E/Eₑ-1.0) <= tol - @test error_ϕ <= tol - @test error_η <= tol +function test_operator(op) + # Solution + sol_t = solve(ode_solver,op,x₀,t₀,tf) + + # Post-process + l2_Ω(w) = √(∑(∫(w*w)dΩ)) + l2_Γ(v) = √(∑(∫(v*v)dΓ)) + E_kin(w) = 0.5*∑( ∫(∇(w)⋅∇(w))dΩ ) + E_pot(v) = g*0.5*∑( ∫(v*v)dΓ ) + Eₑ = 0.5*g*ξ^2*L + + tol = 1.0e-2 + for ((ϕn,ηn),tn) in sol_t + E = E_kin(ϕn) + E_pot(ηn) + error_ϕ = l2_Ω(ϕn-ϕₑ(tn)) + error_η = l2_Γ(ηn-ηₑ(tn)) + @test abs(E/Eₑ-1.0) <= tol + @test error_ϕ <= tol + @test error_η <= tol + end + println("xxx") end +test_operator(op_const) +test_operator(op_trans) +# test_operator(op_ad) # Not working yet + end From 3044f0fff1633b19a473a8b92913333b49fcdb50 Mon Sep 17 00:00:00 2001 From: Oriol Colomes Date: Thu, 9 Dec 2021 22:39:09 +0100 Subject: [PATCH 13/15] Fixing autodiff with new API --- src/TransientFETools/TransientFEOperators.jl | 7 +++---- .../FreeSurfacePotentialFlowTests.jl | 1 - test/TransientFEsTests/HeatEquationAutoDiffTests.jl | 12 +++++------- .../TransientFEsTests/StokesEquationAutoDiffTests.jl | 9 +++++---- 4 files changed, 13 insertions(+), 16 deletions(-) diff --git a/src/TransientFETools/TransientFEOperators.jl b/src/TransientFETools/TransientFEOperators.jl index ea712b6..9346c63 100644 --- a/src/TransientFETools/TransientFEOperators.jl +++ b/src/TransientFETools/TransientFEOperators.jl @@ -204,15 +204,14 @@ function TransientFEOperator(res::Function,trial,test;order::Integer=1) for i in 1:order function jac_i(t,x,dxi,dv) function res_i(y) - derivatives = (x.derivatives[1:i-1],y,x.derivatives[i+1:end]) - res(t,x,dv) + derivatives = (x.derivatives[1:i-1]...,y,x.derivatives[i+1:end]...) + xi = TransientCellField(x.cellfield,derivatives) + res(t,xi,dv) end jacobian(res_i,x.derivatives[i]) end jacs = (jacs...,jac_i) end - # jac(t,u,ut,du,dv) = jacobian(x->res(t,x,ut,dv),u) - # jac_t(t,u,ut,dut,dv) = jacobian(xt->res(t,u,xt,dv),ut) TransientFEOperator(res,jacs...,trial,test) end diff --git a/test/TransientFEsTests/FreeSurfacePotentialFlowTests.jl b/test/TransientFEsTests/FreeSurfacePotentialFlowTests.jl index b222544..c1bf7d5 100644 --- a/test/TransientFEsTests/FreeSurfacePotentialFlowTests.jl +++ b/test/TransientFEsTests/FreeSurfacePotentialFlowTests.jl @@ -96,7 +96,6 @@ function test_operator(op) @test error_ϕ <= tol @test error_η <= tol end - println("xxx") end test_operator(op_const) diff --git a/test/TransientFEsTests/HeatEquationAutoDiffTests.jl b/test/TransientFEsTests/HeatEquationAutoDiffTests.jl index bf9a17b..0d549ad 100644 --- a/test/TransientFEsTests/HeatEquationAutoDiffTests.jl +++ b/test/TransientFEsTests/HeatEquationAutoDiffTests.jl @@ -9,9 +9,6 @@ using GridapODEs.TransientFETools using Gridap.FESpaces: get_algebraic_operator using Gridap.Arrays: test_array -import Gridap: ∇ -import GridapODEs.TransientFETools: ∂t - θ = 0.2 u(x,t) = (1.0-x[1])*x[1]*(1.0-x[2])*x[2]*t @@ -49,12 +46,13 @@ U₀ = evaluate(U,nothing) dv = get_fe_basis(V0) du = get_trial_fe_basis(U₀) uh = FEFunction(U₀,rand(num_free_dofs(U₀))) +uh_t = TransientCellField(uh,(uh,)) -cell_j = get_array(jac(0.5,(uh,uh),du,dv)) -cell_j_t = get_array(jac_t(0.5,(uh,uh),du,dv)) +cell_j = get_array(jac(0.5,uh_t,du,dv)) +cell_j_t = get_array(jac_t(0.5,uh_t,du,dv)) -cell_j_auto = get_array(jacobian(x->res(0.5,(x,uh),dv),uh)) -cell_j_t_auto = get_array(jacobian(x->res(0.5,(uh,x),dv),uh)) +cell_j_auto = get_array(jacobian(x->res(0.5,TransientCellField(x,(uh,)),dv),uh)) +cell_j_t_auto = get_array(jacobian(x->res(0.5,TransientCellField(uh,(x,)),dv),uh)) test_array(cell_j_auto,cell_j,≈) test_array(cell_j_t_auto,cell_j_t,≈) diff --git a/test/TransientFEsTests/StokesEquationAutoDiffTests.jl b/test/TransientFEsTests/StokesEquationAutoDiffTests.jl index a6da49f..55dbcae 100644 --- a/test/TransientFEsTests/StokesEquationAutoDiffTests.jl +++ b/test/TransientFEsTests/StokesEquationAutoDiffTests.jl @@ -75,12 +75,13 @@ X₀ = evaluate(X,nothing) dy = get_fe_basis(Y) dx = get_trial_fe_basis(X₀) xh = FEFunction(X₀,rand(num_free_dofs(X₀))) +xh_t = TransientCellField(xh,(xh,)) -cell_j = get_array(jac(0.5,(xh,xh),dx,dy)) -cell_j_t = get_array(jac_t(0.5,(xh,xh),dx,dy)) +cell_j = get_array(jac(0.5,xh_t,dx,dy)) +cell_j_t = get_array(jac_t(0.5,xh_t,dx,dy)) -cell_j_auto = get_array(jacobian(x->res(0.5,(x,xh),dy),xh)) -cell_j_t_auto = get_array(jacobian(x->res(0.5,(xh,x),dy),xh)) +cell_j_auto = get_array(jacobian(x->res(0.5,TransientCellField(x,(xh,)),dy),xh)) +cell_j_t_auto = get_array(jacobian(x->res(0.5,TransientCellField(xh,(x,)),dy),xh)) for i in 1:length(cell_j) test_array(cell_j[i].array[1,1],cell_j_auto[i].array[1,1],≈) From cdbec51a4cca0b76916f112399108241bd6d6af6 Mon Sep 17 00:00:00 2001 From: Oriol Colomes Date: Thu, 9 Dec 2021 22:54:47 +0100 Subject: [PATCH 14/15] Adapted README.md snipped --- README.md | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/README.md b/README.md index 6a2df6b..f871fae 100644 --- a/README.md +++ b/README.md @@ -14,9 +14,6 @@ using ForwardDiff using GridapODEs.ODETools using GridapODEs.TransientFETools -import Gridap: ∇ -import GridapODEs.TransientFETools: ∂t - θ = 0.5 u(x,t) = (1.0-x[1])*x[1]*(1.0-x[2])*x[2]*t @@ -48,9 +45,9 @@ a(u,v) = ∫( ∇(v)⋅∇(u) )dΩ b(v,t) = ∫( v*f(t) )dΩ m(u,v) = ∫( v*u )dΩ -res(t,(u,ut),v) = a(u,v) + m(ut,v) - b(v,t) -jac(t,(u,ut),du,v) = a(du,v) -jac_t(t,(u,ut),dut,v) = m(dut,v) +res(t,u,v) = a(u,v) + m(∂t(u),v) - b(v,t) +jac(t,u,du,v) = a(du,v) +jac_t(t,u,dut,v) = m(dut,v) op = TransientFEOperator(res,jac,jac_t,U,V0) From 6cebe513edbac3e258a6b6a7533b84390dfbfb64 Mon Sep 17 00:00:00 2001 From: Oriol Colomes Date: Thu, 9 Dec 2021 22:56:31 +0100 Subject: [PATCH 15/15] removing tmp file --- tmp.jl | 305 --------------------------------------------------------- 1 file changed, 305 deletions(-) delete mode 100644 tmp.jl diff --git a/tmp.jl b/tmp.jl deleted file mode 100644 index 988da7d..0000000 --- a/tmp.jl +++ /dev/null @@ -1,305 +0,0 @@ -using Gridap -using Gridap.Arrays -using Gridap.Geometry -using GridapODEs.TransientFETools -using GridapODEs.ODETools -using Test -using LineSearches: BackTracking -using ForwardDiff -import GridapODEs.TransientFETools: ∂t - -# Problem setting -println("Setting analytical fluid problem parameters") -# Solid properties -E_s = 1.0 -ν_s = 0.4 -ρ_s = 1.0 -# Fluid properties -ρ_f = 1.0 -μ_f = 1.0 -γ_f = 1.0 -# Mesh properties -n_m = 3 -E_m = 1.0 -ν_m = 0.1 -α_m = 1.0e-5 -# Time stepping -t0 = 0.0 -tf = 0.5 -dt = 0.1 -θ = 0.5 - -# Define BC functions -println("Defining Boundary conditions") -#u(x, t) = 1.0e-2 * VectorValue( x[1]^2*x[2] , -x[1]*x[2]^2) * t -#v(x, t) = 1.0e-2 * VectorValue( x[1]^2*x[2] , -x[1]*x[2]^2) -u(x, t) = 1.0e-2 * VectorValue( 1.0 , -1.0) * t -v(x, t) = 1.0e-2 * VectorValue( 1.0 , -1.0) -u(t::Real) = x -> u(x,t) -v(t::Real) = x -> v(x,t) -∂tu(t) = x -> VectorValue(ForwardDiff.derivative(t -> get_array(u(x,t)),t)) -∂tv(t) = x -> VectorValue(ForwardDiff.derivative(t -> get_array(v(x,t)),t)) -∂tu(x,t) = ∂tu(t)(x) -∂tv(x,t) = ∂tv(t)(x) -T_u = typeof(u) -T_v = typeof(v) -@eval ∂t(::$T_u) = $∂tu -@eval ∂t(::$T_v) = $∂tv -p(x,t) = 0.0 -p(t::Real) = x -> p(x,t) -bconds = ( - # Tags - FSI_Vu_tags = ["boundary"], - FSI_Vv_tags = ["boundary"], - ST_Vu_tags = ["boundary","interface"], - ST_Vv_tags = ["boundary","interface"], - # Values, - FSI_Vu_values = [u], - FSI_Vv_values = [v], - ST_Vu_values = [u(0.0),u(0.0)], - ST_Vv_values = [v(0.0),v(0.0)], -) - -function lame_parameters(E,ν) - λ = (E*ν)/((1+ν)*(1-2*ν)) - μ = E/(2*(1+ν)) - (λ, μ) -end - -# Define Forcing terms -I = TensorValue( 1.0, 0.0, 0.0, 1.0 ) -F(t) = x -> ∇(u(t))(x) + I -J(t) = x -> det(F(t))(x) -E(t) = x -> 0.5 * ((F(t)(x)')⋅F(t)(x) - I) -(λ_s, μ_s) = lame_parameters(E_s,ν_s) -(λ_m, μ_m) = lame_parameters(E_m,ν_m) -S_SV(t) = x -> 2*μ_s*E(t)(x) + λ_s*tr(E(t)(x))*I -fv_ST_Ωf(t) = x -> - μ_f*Δ(v(t))(x) + ∇(p(t))(x) -fu_Ωf(t) = x -> - α_m * Δ(u(t))(x) -fv_Ωf(t) = x -> ρ_f * ∂t(v)(t)(x) - μ_f * Δ(v(t))(x) + ∇(p(t))(x) + ρ_f*( (∇(v(t))(x)')⋅(v(t)(x) - ∂t(u)(t)(x)) ) -fp_Ωf(t) = x -> (∇⋅v(t))(x) -fu_Ωs(t) = x -> ∂t(u)(t)(x) - v(t)(x) -fv_Ωs(t) = x -> ρ_s * ∂t(v)(t)(x) #- (∇⋅(F(t)⋅S_SV(t)))(x) # Divergence of a a doted function not working yet... - -# Discrete model -println("Defining discrete model") -domain = (-1,1,-1,1) -partition = (n_m,n_m) -model = CartesianDiscreteModel(domain,partition) -trian = Triangulation(model) -R = 0.5 -xc = 0.0 -yc = 0.0 -function is_in(coords) - n = length(coords) - x = (1/n)*sum(coords) - d = (x[1]-xc)^2 + (x[2]-yc)^2 - R^2 - d < 1.0e-8 -end -oldcell_to_coods = get_cell_coordinates(trian) -oldcell_to_is_in = collect1d(lazy_map(is_in,oldcell_to_coods)) -incell_to_cell = findall(oldcell_to_is_in) -outcell_to_cell = findall(collect(Bool, .! oldcell_to_is_in)) -model_solid = DiscreteModel(model,incell_to_cell) -model_fluid = DiscreteModel(model,outcell_to_cell) - -# Build fluid-solid interface labelling -println("Defining Fluid-Solid interface") -labeling = get_face_labeling(model_fluid) -new_entity = num_entities(labeling) + 1 -topo = get_grid_topology(model_fluid) -D = num_cell_dims(model_fluid) -for d in 0:D-1 - fluid_boundary_mask = collect(Bool,get_isboundary_face(topo,d)) - fluid_outer_boundary_mask = get_face_mask(labeling,"boundary",d) - fluid_interface_mask = collect(Bool,fluid_boundary_mask .!= fluid_outer_boundary_mask) - dface_list = findall(fluid_interface_mask) - for face in dface_list - labeling.d_to_dface_to_entity[d+1][face] = new_entity - end -end -add_tag!(labeling,"interface",[new_entity]) - -# Triangulations -println("Defining triangulations") -Ω = Triangulation(model) -Ω_s = Triangulation(model_solid) -Ω_f = Triangulation(model_fluid) -Γ_i = BoundaryTriangulation(model_fluid,tags="interface") - -# Quadratures -println("Defining quadratures") -order = 2 -degree = 2*order -bdegree = 2*order -dΩ = Measure(Ω,degree) -dΩs = Measure(Ω_s,degree) -dΩf = Measure(Ω_f,degree) -dΓi = Measure(Γ_i,bdegree) - -# Test FE Spaces -println("Defining FE spaces") -# Reference FE -reffeᵤ = ReferenceFE(lagrangian,VectorValue{2,Float64},order) -reffeₚ = ReferenceFE(lagrangian,Float64,order-1) -# Test FE Spaces -Vu_FSI = TestFESpace(model, reffeᵤ, conformity=:H1, dirichlet_tags=bconds[:FSI_Vu_tags]) -Vv_FSI = TestFESpace(model, reffeᵤ, conformity=:H1, dirichlet_tags=bconds[:FSI_Vv_tags]) -Vu_ST = TestFESpace(model_fluid, reffeᵤ, conformity=:H1, dirichlet_tags=bconds[:ST_Vu_tags]) -Vv_ST = TestFESpace(model_fluid, reffeᵤ, conformity=:H1, dirichlet_tags=bconds[:ST_Vv_tags]) -Q = TestFESpace(model_fluid, reffeₚ, conformity=:C0, constraint=:zeromean) - -# Trial FE Spaces -Uu_ST = TrialFESpace(Vu_ST,bconds[:ST_Vu_values]) -Uv_ST = TrialFESpace(Vv_ST,bconds[:ST_Vv_values]) -Uu_FSI = TransientTrialFESpace(Vu_FSI,bconds[:FSI_Vu_values]) -Uv_FSI = TransientTrialFESpace(Vv_FSI,bconds[:FSI_Vv_values]) -P = TrialFESpace(Q) - -# Multifield FE Spaces -Y_ST = MultiFieldFESpace([Vu_ST,Vv_ST,Q]) -X_ST = MultiFieldFESpace([Uu_ST,Uv_ST,P]) -Y_FSI = MultiFieldFESpace([Vu_FSI,Vv_FSI,Q]) -X_FSI = TransientMultiFieldFESpace([Uu_FSI,Uv_FSI,P]) - -# Stokes problem for initial solution -println("Defining Stokes operator") -a((u,v,p),(ϕ,φ,q)) = ∫( ϕ⋅u + ε(φ) ⊙ (2.0*μ_f*ε(v)) - (∇⋅φ) * p + q * (∇⋅v) )dΩf -l((ϕ,φ,q)) = ∫( φ⋅fv_ST_Ωf(0.0) )dΩf -res(x,y) = a(x,y) - l(y) -jac(x,dx,y) = a(dx,y) -op_ST = FEOperator(res,jac,X_ST,Y_ST) - -# Map operations -F(∇u) = ∇u + I -J(∇u) = det(F(∇u)) -Finv(∇u) = inv(F(∇u)) -FinvT(∇u) = (Finv(∇u)') -E(∇u) = 0.5 * ((F(∇u)')⋅F(∇u) - I) - -# Map operations derivatives -dF(∇du) = ∇du -dJ(∇u,∇du) = J(∇u)*tr(Finv(∇u)⋅dF(∇du)) -dFinv(∇u,∇du) = -Finv(∇u) ⋅ dF(∇du) ⋅ Finv(∇u) -dFinvT(∇u,∇du) = (dFinv(∇u,∇du)') -dE(∇u,∇du) = 0.5 * ((dF(∇du)')⋅F(∇u) + (F(∇u)')⋅dF(∇du)) - -# Fluid constitutive laws -# Cauchy stress -σᵥ_Ωf(μ,u) = 2.0*μ*ε(u) -σᵥ_Ωf(μ) = (∇v,Finv) -> μ*(∇v⋅Finv + (Finv')⋅(∇v')) -# First Piola-Kirchhoff stress -Pᵥ_Ωf(μ,u,v) = (J∘∇(u)) * (σᵥ_Ωf(μ)∘(∇(v),Finv∘∇(u)))' ⋅ (FinvT∘∇(u)) -function Pₚ_Ωf(u,p) - typeof(- (J∘∇(u)) * p * tr((FinvT∘∇(u)))) - - (J∘∇(u)) * p * tr((FinvT∘∇(u))) -end -# First Piola-Kirchhoff stress Jacobian -dPᵥ_Ωf_du(μ,u,du,v) = (dJ∘(∇(u),∇(du))) * (σᵥ_Ωf(μ)∘(∇(v),(Finv∘∇(u))))' ⋅ (FinvT∘∇(u)) + - (J∘∇(u)) * (σᵥ_Ωf(μ)∘(∇(v),(dFinv∘(∇(u),∇(du)))))' ⋅ (FinvT∘∇(u)) + - (J∘∇(u)) * (σᵥ_Ωf(μ)∘(∇(v),(Finv∘∇(u))))' ⋅ (dFinvT∘(∇(u),∇(du))) -dPₚ_Ωf_du(u,du,p) = - p * ( (dJ∘(∇(u),∇(du))) * tr((FinvT∘∇(u))) + - (J∘∇(u)) * tr((dFinvT∘(∇(u),∇(du)))) ) -dPᵥ_Ωf_dv(μ,u,dv) = Pᵥ_Ωf(μ,u,dv) -dPₚ_Ωf_dp(u,dp) = Pₚ_Ωf(u,dp) -# Convective term -conv(c,∇v) = (∇v') ⋅ c -dconv(dc,∇dv,c,∇v) = conv(c,∇dv) + conv(dc,∇v) - -# Solid constitutive laws -# Second Piola-Kirchhoff stress (Saint-Venant) -Sₛᵥ_Ωs(λ,μ,u) = 2*μ*(E∘∇(u)) + λ*tr((E∘∇(u)))*(I∘∇(u)) -dSₛᵥ_Ωs_du(λ,μ,u,du) = 2*μ*(dE∘(∇(u),∇(du))) + λ*tr((dE∘(∇(u),∇(du))))*(I∘∇(u)) -# First Piola-Kirchhoff stress (Saint-Venant) -Pₛᵥ_Ωs(λ,μ,u) = (F∘∇(u)) ⋅ Sₛᵥ_Ωs(λ,μ,u) -dPₛᵥ_Ωs_du(λ,μ,u,du) = (dF∘∇(du)) ⋅ Sₛᵥ_Ωs(λ,μ,u) + (F∘∇(u)) ⋅ dSₛᵥ_Ωs_du(λ,μ,u,du) - -# # FSI problem -println("Defining FSI operator") -n_Γi = get_normal_vector(Γ_i) -a_f((u,v,p),(ϕ,φ,q)) = ∫( φ ⋅ ( (J∘∇(u)) * ρ_f * ∂t(v) ) )dΩ -da_f((u,v,p),(du,dv,dp),(ϕ,φ,q)) = ∫( φ ⋅ ( (dJ∘(∇(u),∇(du))) * ρ_f * ∂t(v) ) )dΩ -da_t_f((u,v,p),(dut,dvt,dpt),(ϕ,φ,q)) = ∫( φ ⋅ ( (J∘∇(u)) * ρ_f * dvt ) )dΩ -# da_t_f((u,v,p),(dut,dvt,dpt),(ϕ,φ,q)) = ∫( φ ⋅ ( ρ_f * dvt) )dΩf -a_fsi((u,v,p),(ϕ,φ,q)) = ∫( ϕ⋅u + φ⋅v + q*p )dΩ -l_f(t,(ϕ,φ,q)) = ∫( ϕ⋅fu_Ωf(t) )dΩ -l_fsi(t,(ϕ,φ,q)) = ∫( φ⋅fv_Ωf(t) )dΩ -res(t,x,y) = a_f(x,y) + a_fsi(x,y) + a_fsi(∂t(x),y) - l_f(t,y) - l_fsi(t,y) -jac(t,x,dx,y) = da_f(x,dx,y) + a_fsi(dx,y) -jac_t(t,x,dxt,y) = da_t_f(x,dxt,y) + a_fsi(dxt,y) -op_FSI = TransientFEOperator(res,jac,jac_t,X_FSI,Y_FSI) -#op_FSI = TransientFEOperator(res,X_FSI,Y_FSI) -# a_f((u,v,p),(ut,vt,pt),(ϕ,φ,q)) = -# ∫( α_m * (∇(ϕ) ⊙ ∇(u)) )dΩf + -# ∫( φ ⋅ ( (J∘∇(u)) * ρ_f * vt ) )dΩf -# da_f((u,v,p),(ut,vt,pt),(du,dv,dp),(ϕ,φ,q)) = -# ∫( α_m * (∇(ϕ) ⊙ ∇(du)) )dΩf + -# ∫( φ ⋅ ( (dJ∘(∇(u),∇(du))) * ρ_f * vt ) )dΩf -# da_t_f((u,v,p),(dut,dvt,dpt),(ϕ,φ,q)) = ∫( φ ⋅ ( (J∘∇(u)) * ρ_f * dvt ) )dΩf -# a_fsi((u,v,p),(ϕ,φ,q)) = ∫( ϕ⋅u + φ⋅v + q*p )dΩf + ∫( ϕ⋅u + φ⋅v + q*p )dΩs -# l_f(t,(ϕ,φ,q)) = ∫( ϕ⋅fu_Ωf(t) )dΩf -# l_fsi(t,(ϕ,φ,q)) = ∫( φ⋅fv_Ωf(t) )dΩf -# res(t,(x,xt),y) = a_f(x,xt,y) + a_fsi(x,y) + a_fsi(xt,y) - l_f(t,y) - l_fsi(t,y) -# jac(t,(x,xt),dx,y) = da_f(x,xt,dx,y) + a_fsi(dx,y) -# jac_t(t,(x,xt),dxt,y) = da_t_f(x,dxt,y) + a_fsi(dxt,y) -# op_FSI = TransientFEOperator(res,jac,jac_t,X_FSI,Y_FSI) - -# Setup output files -folderName = "fsi-results" -fileName = "fields" -if !isdir(folderName) - mkdir(folderName) -end -filePath = join([folderName, fileName], "/") - -# Solve Stokes problem -println("Defining Stokes solver") -xh = solve(op_ST) -writevtk(Ω_f, filePath, cellfields = ["uh"=>xh[1], "vh"=>xh[2], "ph"=>xh[3]]) - -# Compute Stokes solution L2-norm -l2(w) = w⋅w -eu_ST = u(0.0) - xh[1] -ev_ST = v(0.0) - xh[2] -ep_ST = p(0.0) - xh[3] -eul2_ST = sqrt(∑( ∫(l2(eu_ST))dΩf )) -evl2_ST = sqrt(∑( ∫(l2(ev_ST))dΩf )) -epl2_ST = sqrt(∑( ∫(l2(ep_ST))dΩf )) -println("Stokes L2-norm u: ", eul2_ST) -println("Stokes L2-norm v: ", evl2_ST) -println("Stokes L2-norm p: ", epl2_ST) -@test eul2_ST < 1.0e-10 -@test evl2_ST < 1.0e-10 -@test epl2_ST < 1.0e-10 - -# Solve FSI problem -println("Defining FSI solver") -xh0 = interpolate_everywhere([u(0.0),v(0.0),p(0.0)],X_FSI(0.0)) -nls = NLSolver( -show_trace = true, -method = :newton, -linesearch = BackTracking(), -ftol = 1.0e-10, -iterations = 50 -) -odes = ThetaMethod(nls, dt, 0.5) -solver = TransientFESolver(odes) -xht = solve(solver, op_FSI, xh0, t0, tf) - -for (i,(xh, t)) in enumerate(xht) - println("STEP: $i, TIME: $t") - println("============================") - - # Compute errors - eu = u(t) - xh[1] - ev = v(t) - xh[2] - ep = p(t) - xh[3] - eul2 = sqrt(∑( ∫(l2(eu))dΩ )) - evl2 = sqrt(∑( ∫(l2(ev))dΩ )) - epl2 = sqrt(∑( ∫(l2(ep))dΩ )) - - # Test checks - println("FSI L2-norm u: ", eul2) - println("FSI L2-norm v: ", evl2) - println("FSI L2-norm p: ", epl2) -end \ No newline at end of file