Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Transient nsi #23

Open
wants to merge 44 commits into
base: master
Choose a base branch
from
Open
Changes from 1 commit
Commits
Show all changes
44 commits
Select commit Hold shift + click to select a range
88b7744
Added 2D cylinder benchmark GMSH model
May 27, 2020
53e4ef2
First draft of the transient incompressible Navier-Stokes tutorial
May 27, 2020
cdff0f9
- Added cylinder.gif
oriolcg Jun 4, 2020
b7be42a
Added Drag and Lift coefficients plot to transient_ins tutorial
oriolcg Jun 4, 2020
29ea29f
Formating for documentation
oriolcg Jun 4, 2020
911a2e9
Updated Project.toml
oriolcg Jun 5, 2020
5a6a228
Merge branch 'master' of https://github.com/gridap/Tutorials into tra…
oriolcg Jun 5, 2020
8966dd7
First (non working) version of transient FSI problem.
oriolcg Jun 10, 2020
98d4c7c
Added an auxiliar transient_fsi_debug.jl for debugging purposes.
oriolcg Jun 10, 2020
badf074
Merge branch 'correct_normals' into transient_nsi
oriolcg Jun 10, 2020
9eca4b6
Corrected $n * \sigma$ in transient_inc_navier_stokes.jl
oriolcg Jun 10, 2020
936e09f
Merge branch 'correct_normals' into transient_nsi
oriolcg Jun 10, 2020
c5b9d91
Merge branch 'preparing_version_0.11' of https://github.com/gridap/Tu…
oriolcg Jun 11, 2020
b467202
Updated transient_inc_navier_stokes.jl tutorial to the new product co…
oriolcg Jun 11, 2020
6194356
Merge branch 'transient_nsi' of https://github.com/gridap/Tutorials i…
oriolcg Jun 11, 2020
2e639b7
Debugging transient fsi tutorial
oriolcg Jun 25, 2020
f5fcd15
further developments on the transient fsi tutorial
oriolcg Jul 3, 2020
43caad8
Further developments in transient_fsi_debug.jl
oriolcg Jul 10, 2020
08f8377
Merge branch 'master' of https://github.com/gridap/Tutorials into tra…
Aug 3, 2020
264e465
Updated transient NS tutorial
Aug 3, 2020
69d0471
Minor modifications in transient NSI tutorial.
Aug 4, 2020
69f2a28
Merge branch 'master' of https://github.com/gridap/Tutorials into tra…
oriolcg Feb 3, 2021
796941c
Added *pvd files in .gitignore
oriolcg Feb 3, 2021
1d87637
Updated transient_fsi tutorial (not working yet)
oriolcg Feb 4, 2021
7575992
Deleted temporary files
oriolcg Feb 16, 2021
c5b3091
Merge branch 'master' of https://github.com/gridap/Tutorials into tra…
oriolcg May 25, 2021
8bd4a1b
Working on the transient tutorials
oriolcg May 25, 2021
9a27679
minor changes in transient INS
oriolcg May 25, 2021
0c61594
further devs in transient_inc_NS
oriolcg May 27, 2021
49dc9c5
Temporarily adding manifest
fverdugo Jun 3, 2021
49ff1fb
Moving to Gridap 0.16
fverdugo Jun 3, 2021
76daf77
Renaming functions according Gridap 0.16
fverdugo Jun 3, 2021
44ee68c
Misc fixes to poisson_dev_fe in order to let
amartinhuertas Jun 4, 2021
e378083
Updated comment in poisson dev fe tutorial
amartinhuertas Jun 4, 2021
57cb559
Updated transient NSI
oriolcg Jun 5, 2021
e58e20c
Merge branch 'master' of https://github.com/gridap/Tutorials into tra…
oriolcg Jun 5, 2021
03a0183
Merge branch 'gridap_v0.16' of https://github.com/gridap/Tutorials in…
oriolcg Jun 5, 2021
6ef3565
Updates in transient NS
oriolcg Jun 8, 2021
d6998fc
modifications in transient NS
oriolcg Jun 10, 2021
5fb4317
Minor changes in transient NS
oriolcg Jun 10, 2021
ba942f9
Merge branch 'transient_fsi' into transient_nsi
oriolcg Jun 10, 2021
fbfa091
Changes in transient NS
oriolcg Jun 10, 2021
0b6ce5e
Link to GridapODES#master
oriolcg Jun 10, 2021
a302e5d
Fix in build.jl
oriolcg Jun 15, 2021
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
Debugging transient fsi tutorial
  • Loading branch information
oriolcg committed Jun 25, 2020
commit 2e639b76778be9b312c38e415104325286f60eeb
298 changes: 220 additions & 78 deletions src/transient_fsi_debug.jl
Original file line number Diff line number Diff line change
@@ -1,9 +1,13 @@
using Gridap
using GridapODEs
using GridapODEs.ODETools
using GridapODEs.TransientFETools
using LinearAlgebra: tr, inv, det
using LineSearches: BackTracking
using WriteVTK
import GridapODEs.TransientFETools: ∂t

model = DiscreteModelFromFile("../runs/elasticFlag_coarse.json")
#model = DiscreteModelFromFile("../runs/elasticFlag_coarse.json")
model = DiscreteModelFromFile("../models/elasticFlag.json")
model_solid = DiscreteModel(model,"solid")
model_fluid = DiscreteModel(model,"fluid")
writevtk(model,"model")
@@ -32,33 +36,35 @@ function lame_parameters(E,ν)
μ = E/(2*(1+ν))
(λ, μ)
end
const (λ_s,μ_s) = lame_parameters(1.0, 0.33)
const (λ_m,μ_m) = lame_parameters(1.0, 0.33)
const μ_f = Um * ⌀ / Re
const ρ_s = 1.0
const ρ_f = 1.0
const (λ_m,μ_m) = lame_parameters(1.0e-6, 0.33)
const (λ_s,μ_s) = lame_parameters(1.4e6, 0.4)
const ρ_s = 1.0e4
const ρ_f = 1.0e3
const μ_f = ρ_f * Um * ⌀ / Re

# ### Constitutive laws
@law F(∇u) = ∇u + one(∇u)
@law J(∇u) = det(F(∇u))
@law Finv(∇u) = inv(F(∇u))
@law FinvT(∇u) = (Finv(∇u)')
@law E(∇u) = 0.5 * ((F(∇u)')*F(∇u) - one(F(∇u)))
@law E(∇u) = 0.5 * ((F(∇u)')F(∇u) - one(F(∇u)))
@law S(∇u) = 2*μ_s*E(∇u) + λ_s*tr(E(∇u))*one(E(∇u))
@law σ_dev(∇v,Finv) = μ_f*(∇v*Finv + (Finv')*(∇v'))
@law σ_dev(∇v,Finv) = μ_f*(∇vFinv + (Finv')(∇v'))
@law σ_m(ε) = λ_m*tr(ε)*one(ε) + 2*μ_m*ε
@law conv(c,∇v) = (∇v') * c
@law dconv(dc,∇dv,c,∇v) = conv(c,∇dv) + conv(dc,∇v)
@law conv(c,∇v) = (∇v') c
@law Sm(∇u) = 2*μ_m*E(∇u) + λ_m*tr(E(∇u))*one(E(∇u))

# Derivatives:
@law dF(∇du) = ∇du
@law dJ(∇u,∇du) = J(F(∇u))*tr(inv(F(∇u))*dF(∇du))
@law dE(∇u,∇du) = 0.5 * ((dF(∇du)')*F(∇u) + (F(∇u)')*dF(∇du))
@law dJ(∇u,∇du) = J(F(∇u))*tr(inv(F(∇u))dF(∇du))
@law dE(∇u,∇du) = 0.5 * ((dF(∇du)')F(∇u) + (F(∇u)')dF(∇du))
@law dS(∇u,∇du) = 2*μ_s*dE(∇u,∇du) + λ_s*tr(dE(∇u,∇du))*one(E(∇u))
@law dFinv(∇u,∇du) = -Finv(∇u) * dF(∇du) * Finv(∇u)
@law dFinv(∇u,∇du) = -Finv(∇u) dF(∇du) Finv(∇u)
@law dFinvT(∇u,∇du) = (dFinv(∇u,∇du)')
@law dconv(dc,∇dv,c,∇v) = conv(c,∇dv) + conv(dc,∇v)
@law dSm(∇u,∇du) = 2*μ_m*dE(∇u,∇du) + λ_m*tr(dE(∇u,∇du))*one(E(∇u))


# Test FE Spaces
k = 2
Vu = TestFESpace(
model=model,
@@ -76,123 +82,259 @@ Vv = TestFESpace(
dirichlet_tags=["inlet", "noSlip", "cylinder","fixed"])
Vuf = TestFESpace(
model=model_fluid,
#model=model,
valuetype=VectorValue{2,Float64},
reffe=:Lagrangian,
order=k,
conformity =:H1,
dirichlet_tags=["inlet", "noSlip", "cylinder","interface","outlet"])
Vvf = TestFESpace(
model=model_fluid,
#model=model,
valuetype=VectorValue{2,Float64},
reffe=:Lagrangian,
order=k,
conformity =:H1,
dirichlet_tags=["inlet", "noSlip", "cylinder","interface"])
Qf = TestFESpace(
model=model_fluid,
#model=model,
valuetype=Float64,
order=k-1,
reffe=:Lagrangian,
conformity=:C0)

# Trial FE Spaces
Uuf = TrialFESpace(Vuf,[u_noSlip(0.0), u_noSlip(0.0), u_noSlip(0.0), u_noSlip(0.0), u_noSlip(0.0)])
Uvf = TrialFESpace(Vvf,[u_in(0.0), u_noSlip(0.0), u_noSlip(0.0), u_noSlip(0.0)])
Utu = TransientTrialFESpace(Vu,[u_noSlip, u_noSlip, u_noSlip, u_noSlip, u_noSlip])
Utv = TransientTrialFESpace(Vv,[u_in, u_noSlip, u_noSlip, u_noSlip])
Pf = TrialFESpace(Qf)

# Multifield FE Spaces
Y = MultiFieldFESpace([Vuf,Vvf,Qf])
X = MultiFieldFESpace([Uuf,Uvf,Pf])
Yt = MultiFieldFESpace([Vu,Vv,Qf])
Xt = MultiFieldFESpace([Ut_u,Ut_v,Pf])
Xt = MultiFieldFESpace([Utu,Utv,Pf])

function a_vp_f(x, y)
v, p = x
φ, q = y
inner(ε(φ), 2*μ_f*ε(v)) - (∇ * φ) * p + q * (∇ * v)
# Stokes Bilinear forms
function aST_ϕ_f(x, y)
u, v, p = x
ϕ, φ, q = y
(∇(ϕ) ⊙ σ_m(ε(u)))
end

function a_uvp_f(x, y)
function aST_φ_f(x, y)
u, v, p = x
ϕ, φ, q = y
inner(∇(ϕ), σ_m(ε(u))) + inner(ε(φ), 2*μ_f*ε(v)) - (∇ * φ) * p + q * (∇ * v)
visc_term = ε(φ) ⊙ ( 2*μ_f*ε(v) )
pres_term = (∇⋅φ) * p
visc_term + pres_term
end

function at_uvp_f(x, xt, y)
function aST_q_f(x, y)
u, v, p = x
ut, vt, pt = xt
ϕ, φ, q = y
inner(φ,vt)
q * (∇⋅v)
end

function a_uv_s(x, y)
u,v = x
ϕ,φ = y
inner(ε(ϕ), σ_m(ε(u))) + inner(ε(φ), 2*μ_f*ε(v))
# FSI Bilinear forms
function aFSI_ϕ_f(x, y)
u, v, p = x
ϕ, φ, q = y
(∇(ϕ) ⊙ (F(∇(u))⋅Sm(∇(u))))
#(∇(ϕ) ⊙ σ_m(ε(u)))
end

function at_uv_s(x, xt, y)
u,v = x
ut, vt = xt
ϕ,φ = y
inner(ϕ,ut) + inner(φ,vt)
function aFSI_ϕ_Γi(x,y)
u,v,p = x
ϕ,φ,q = y
- (ϕ ⋅ (n_Γi⋅(F(∇(u))⋅Sm(∇(u)))) )
#- (ϕ ⋅ (n_Γi⋅σ_m(ε(u))))
end
function aFSI_φ_f(x, xt, y)
u, v, p = x
ut, vt, pt = xt
ϕ, φ, q = y
temp_term = φ ⋅ ( J(∇(u)) * ρ_f * vt )
conv_term = φ ⋅ ( J(∇(u)) * ρ_f * conv( Finv(∇(u))⋅(v-ut), ∇(v)) )
visc_term = ( ∇(φ) ⊙ ( J(∇(u)) * σ_dev(∇(v),Finv(∇(u))) ⋅ FinvT(∇(u))) )
pres_term = (∇⋅φ) * J(∇(u)) * p * tr(FinvT(∇(u)))
temp_term + conv_term + visc_term + pres_term
end
function aFSI_q_f(x, y)
u, v, p = x
ϕ, φ, q = y
q * (∇⋅v)
#q * (∇⋅( (J(∇(u))*Finv(∇(u)))⋅v ))

res_uvp_f(x,y) = a_uvp_f(x,y)
res_uvp_f(t,x,xt,y) = at_uvp_f(x,xt,y) + a_uvp_f(x,y)
res_uv_s(t,x,xt,y) = at_uvp_s(x,xt,y) + a_uvp_s(x,y)
end
function aFSI_ϕ_s(x, xt, y)
u,v,p = x
ut, vt,pt = xt
ϕ,φ,q = y
(ϕ⋅ut) + 0.0*(u⋅ϕ) - (ϕ⋅v)
end
function aFSI_φ_s(x, xt, y)
u,v,p = x
ut, vt,pt = xt
ϕ,φ,q = y
(φ⋅(ρ_s*vt)) + 0.0*(φ⋅(ρ_s*v)) + (∇(φ) ⊙ (F(∇(u))⋅S(∇(u))))
end

jac_uvp_f(x,dx,y) = a_uvp_f(dx,y)
jac_uvp_f(t,x,xt,dx,y) = a_uvp_f(dx,y)
jac_t_uvp_f(t,x,xt,dxt,y) = at_uvp_f(x,dxt,y)
jac_uv_s(t,x,xt,dx,y) = a_uv_s(dx,y)
jac_t_uv_s(t,x,xt,dxt,y) = at_uv_s(x,dxt,y)
# FSI Jacobians
function daFSI_du_ϕ_f(x,dx, y)
u, v, p = x
du, dv, dp = dx
ϕ, φ, q = y
(∇(ϕ) ⊙ (dF(∇(du))⋅Sm(∇(u)) + F(∇(u))⋅dSm(∇(u),∇(du))))
#(∇(ϕ) ⊙ σ_m(ε(du)))
end
function daFSI_du_ϕ_Γi(x,dx,y)
u,v,p = x
du,dv,dp = dx
ϕ,φ,q = y
- (ϕ ⋅ (n_Γi⋅(dF(∇(du))⋅Sm(∇(u)) + F(∇(u))⋅dSm(∇(u),∇(du)))) )
#- (ϕ ⋅ (n_Γi⋅σ_m(ε(du))) )
end
function daFSI_du_φ_f(x, xt, dx, y)
u, v, p = x
ut, vt, pt = xt
du, dv, dp = dx
ϕ, φ, q = y
temp_term = φ ⋅ ( dJ(∇(u),∇(du)) * ρ_f * vt )
conv_term = φ ⋅ ( ( dJ(∇(u),∇(du)) * ρ_f * conv( Finv(∇(u))⋅(v-ut), ∇(v)) ) +
( J(∇(u)) * ρ_f * conv( dFinv(∇(u),∇(du))⋅(v-ut), ∇(v)) ) )
visc_term = ∇(φ) ⊙ ( dJ(∇(u),∇(du)) * σ_dev(∇(v),Finv(∇(u))) ⋅ FinvT(∇(u)) +
J(∇(u)) * σ_dev(∇(v),dFinv(∇(u),∇(du))) ⋅ FinvT(∇(u)) +
J(∇(u)) * σ_dev(∇(v),Finv(∇(u))) ⋅ dFinvT(∇(u),∇(du)) )
pres_term = (∇⋅φ) * p * ( dJ(∇(u),∇(du)) * tr(FinvT(∇(u))) +
J(∇(u)) * tr(dFinvT(∇(u),∇(du))) )
temp_term + conv_term + visc_term + pres_term
end
function daFSI_dv_φ_f(x, xt, dx, y)
u, v, p = x
ut, vt, pt = xt
du, dv, dp = dx
ϕ, φ, q = y
conv_term = φ ⋅ ( J(∇(u)) * ρ_f * dconv( Finv(∇(u))⋅dv, ∇(dv), Finv(∇(u))⋅(v-ut) , ∇(v)) )
visc_term = ( ∇(φ) ⊙ ( J(∇(u)) * σ_dev(∇(dv),Finv(∇(u))) ⋅ FinvT(∇(u))) )
conv_term + visc_term
end
function daFSI_dp_φ_f(x, dx, y)
u, v, p = x
du, dv, dp = dx
ϕ, φ, q = y
pres_term = (∇⋅φ) * J(∇(u)) * dp * tr(FinvT(∇(u)))
end
function daFSI_dut_φ_f(x, dxt, y)
u, v, p = x
dut, dvt, dpt = dxt
ϕ, φ, q = y
conv_term = - φ ⋅ ( J(∇(u)) * ρ_f * conv( Finv(∇(u))⋅dut, ∇(v)) )
end
function daFSI_dvt_φ_f(x, dxt, y)
u, v, p = x
dut, dvt, dpt = dxt
ϕ, φ, q = y
temp_term = φ ⋅ ( J(∇(u)) * ρ_f * dvt )
end
function daFSI_ϕ_s(x, dx, y)
u,v,p = x
du,dv,dp = dx
ϕ,φ,q = y
0.0*(du⋅ϕ) - (ϕ⋅dv)
end
function daFSI_φ_s(x, dx, y)
u,v,p = x
du,dv,dp = dx
ϕ,φ,q = y
0.0*(φ⋅(ρ_s*dv)) + (∇(φ) ⊙ ( dF(∇(du))⋅S(∇(u)) + (F(∇(u))⋅dS(∇(u),∇(du))) ) )
end
function daFSI_dt_s(x, dxt, y)
u, v, p = x
dut, dvt, dpt = dxt
ϕ, φ, q = y
ϕ⋅dut + (φ⋅(ρ_s*dvt))
end

# Triangulations and quadratures
degree = 2*k
bdegree = 2*k
trian = Triangulation(model)
trian_solid = Triangulation(model_solid)
trian_fluid = Triangulation(model_fluid)
degree = 2*k
trian_Γi = BoundaryTriangulation(model_fluid,"interface")
n_Γi = get_normal_vector(trian_Γi)
quad_solid = CellQuadrature(trian_solid,degree)
quad_fluid = CellQuadrature(trian_fluid,degree)
quad_Γi = CellQuadrature(trian_Γi,bdegree)

t_Ωf = FETerm(res_uvp_f, jac_uvp_f, trian_fluid, quad_fluid)
op = FEOperator(X,Y,t_Ωf)
xh = solve(op)
# Stokes FE Operator
res_ST_f(x,y) = aST_ϕ_f(x,y) + aST_φ_f(x,y) + aST_q_f(x,y)
jac_ST_f(x,dx,y) = aST_ϕ_f(dx,y) + aST_φ_f(dx,y) + aST_q_f(dx,y)
tST_Ωf = FETerm(res_ST_f, jac_ST_f, trian_fluid, quad_fluid)
opST = FEOperator(X,Y,tST_Ωf)

# vh = restrict(xh.blocks[1],trian_fluid)
# ph = restrict(xh.blocks[2],trian_fluid)
# writevtk(trian_fluid,"test_vp_f_0.0.vtu",cellfields = ["vh" => vh, "ph" => ph])
# FSI FE Operator
res_FSI_f(t,x,xt,y) = aFSI_ϕ_f(x,y) + aFSI_φ_f(x,xt,y) + aFSI_q_f(x,y)
jac_FSI_f(t,x,xt,dx,y) = daFSI_du_ϕ_f(x,dx,y) + daFSI_du_φ_f(x,xt,dx,y) + daFSI_dv_φ_f(x,xt,dx,y) + daFSI_dp_φ_f(x,dx,y) + aFSI_q_f(dx,y)
jac_t_FSI_f(t,x,xt,dxt,y) = daFSI_dut_φ_f(x,dxt,y) + daFSI_dvt_φ_f(x,dxt,y)
res_FSI_s(t,x,xt,y) = aFSI_ϕ_s(x,xt,y) + aFSI_φ_s(x,xt,y)
jac_FSI_s(t,x,xt,dx,y) = daFSI_ϕ_s(x,dx,y) + daFSI_φ_s(x,dx,y)
jac_t_FSI_s(t,x,xt,dxt,y) = daFSI_dt_s(x,dxt,y)
tFSI_Ωf = FETerm(res_FSI_f, jac_FSI_f, jac_t_FSI_f, trian_fluid, quad_fluid)
tFSI_Ωs = FETerm(res_FSI_s, jac_FSI_s, jac_t_FSI_s, trian_solid, quad_solid)
tFSI_Γi = FETerm(aFSI_ϕ_Γi,daFSI_du_ϕ_Γi,trian_Γi,quad_Γi)
opFSI = TransientFEOperator(Xt,Yt,tFSI_Ωf,tFSI_Ωs,tFSI_Γi)

# U0u = Utu(0.0)
# uh0 = interpolate_everywhere(U0u, u_noSlip(0.0))
# vh0 = xh[1]
# ph0 = xh[2]
# xh0 = Gridap.MultiField.MultiFieldFEFunction(Xt(0.0),[uh0,vh0,ph0])
# Output function
function writePVD(filePath::String, trian::Triangulation, sol; append=false)
outfiles = paraview_collection(filePath, append=append) do pvd
for (i, (xh, t)) in enumerate(sol)
uh = xh.blocks[1]
vh = xh.blocks[2]
ph = xh.blocks[3]
pvd[t] = createvtk(
trian,
filePath * "_$t.vtu",
cellfields = ["uh" => uh, "vh" => vh, "ph" => ph],
)
end
end
end
folderName = "fsi-results"
fileName = "fields"
if !isdir(folderName)
mkdir(folderName)
end
filePath = join([folderName, fileName], "/")

# Solve Stokes problem
xh = solve(opST)
writePVD(filePath, trian, [(xh, 0.0)])

# Solve FSI problem
xh0 = interpolate(Xt(0.0),xh)
uh0f = restrict(xh0.blocks[1],trian_fluid)
vh0f = restrict(xh0.blocks[2],trian_fluid)
ph0f = restrict(xh0.blocks[3],trian_fluid)
writevtk(trian_fluid,"test_uvp_f_0.0.vtu",cellfields = ["uh" => uh0f, "vh" => vh0f, "ph" => ph0f])
writevtk(trian,"test_uvp_0.0.vtu",cellfields = ["uh" => xh0[1], "vh" => xh0[2], "ph" => xh0[3]])

t_t_Ωf = FETerm(res_uvp_f, jac_uvp_f, jac_t_uvp_f, trian_fluid, quad_fluid)
t_t_Ωs = FETerm(res_uv_s, jac_uv_s, jac_t_uv_s, trian_solid, quad_solid)
opt = TransientFEOperator(Xt,Yt,t_t_Ωf,t_t_Ωs)
ls = LUSolver()
odes = ThetaMethod(ls, 0.001, 0.5)
nls = NLSolver(
show_trace = true,
method = :newton,
linesearch = BackTracking(),
)
odes = ThetaMethod(nls, 0.05, 0.5)
solver = TransientFESolver(odes)
sol_t = solve(solver, opt, xh0, 0.0, 0.002)
println(length(xh0))
for (xh,t) in sol_t
println(typeof(t), length(xh[1]))
# uh = restrict(xh.blocks[1],trian_fluid)
# vh = restrict(xh.blocks[2],trian_fluid)
# ph = restrict(xh.blocks[3],trian_fluid)
# writevtk(trian_fluid,"test_uvp_f_$t.vtu",cellfields = ["uh" => uh, "vh" => vh, "ph" => ph])
writevtk(trian,"test_uvp_$t.vtu",cellfields = ["uh" => xh0[1], "vh" => xh0[2], "ph" => xh0[3]])
end
sol_t = solve(solver, opFSI, xh0, 0.0, 10.0)
writePVD(filePath, trian, sol_t, append=true)


# println(length(xh0.blocks))
# for (xh,t) in sol_t
# println(typeof(t), length(xh[1]))
# # uh = restrict(xh.blocks[1],trian_fluid)
# # vh = restrict(xh.blocks[2],trian_fluid)
# # ph = restrict(xh.blocks[3],trian_fluid)
# # writevtk(trian_fluid,"test_uvp_f_$t.vtu",cellfields = ["uh" => uh, "vh" => vh, "ph" => ph])
# writevtk(trian,"test_uvp_$t.vtu",cellfields = ["uh" => xh0[1], "vh" => xh0[2], "ph" => xh0[3]])
# end

# vh = restrict(xh.blocks[1],trian_fluid)
# ph = restrict(xh.blocks[2],trian_fluid)
# writevtk(trian_fluid,"test.vtu",cellfields = ["vh" => vh, "ph" => ph])