[LeoPart,GeoPart] Variable coefficient in HDG Stokes not works in 2D

1

I’m using geopart and leopart to solve a Stokes problem using HDG with a variable coefficient. I found the same code only works for 3D but not 2D.

Is there any problem with my implementation below?

Mateiral(mu)

Domain.PNG
Domain.PNG825×500 160 KB

2D Velocity solutions:
Sol2D

3D Velocity solutions:

Sol3D
Sol3D704×450 290 KB

Package:

import numpy as np
import dolfin
from dolfin import *

from mpi4py import MPI as pyMPI

from leopart import StokesStaticCondensation, FormsStokes
import geopart.stokes.incompressible
import dolfin_dg as dg

comm = pyMPI.COMM_WORLD
mpi_comm = MPI.comm_world


#load mesh,boundaries and coefficients from file
mark = {"Internal":0, "wall": 1,"inlet": 2,"outlet": 3 }

#read mesh and boundaries from file
mesh = Mesh()
hdf = HDF5File(mesh.mpi_comm(), "mesh_boundaries.h5", "r")
hdf.read(mesh, "/mesh", False)
boundaries = MeshFunction('size_t', mesh, mesh.topology().dim() - 1)
hdf.read(boundaries, "/boundaries")
hdf.close()

#read viscosity coefficient from file
mu_ele = FunctionSpace(mesh, "DG", 0)
mu = Function(mu_ele)
hdf = HDF5File(mesh.mpi_comm(), "mesh_coeffs.h5", "r")
hdf.read(mu, "/mu")
hdf.close()

#output viscosity to paraview
XDMFFile(mpi_comm, "coeff_preview.xdmf").write_checkpoint(mu, "coeffs", 0)

#Define HDG element and function space
element_cls = geopart.stokes.incompressible.HDG2()
W = element_cls.function_space(mesh)

ds = dolfin.Measure('ds',domain=mesh,subdomain_data=boundaries)
n = dolfin.FacetNormal(mesh)

#Define boundary condition
U = element_cls.create_solution_variable(W)

p_in = dolfin.Constant(1.0)         # pressure inlet
p_out = dolfin.Constant(0.0)        # pressure outlet
noslip = dolfin.Constant([0.0]*mesh.geometry().dim())  # no-slip wall

#Boundary conditions
gN1 = (- p_out*dolfin.Identity(mesh.geometry().dim())) * n
Neumann_outlet=dg.DGNeumannBC(ds(mark["outlet"]), gN1)
gN2 = (- p_in*dolfin.Identity(mesh.geometry().dim())) * n
Neumann_inlet=dg.DGNeumannBC(ds(mark["inlet"]), gN2)
Dirichlet_wall=dg.DGDirichletBC(ds(mark["wall"]), noslip)

weak_bcs = [Dirichlet_wall,Neumann_inlet,Neumann_outlet]

#Body force term
f = dolfin.Constant([0.0]*mesh.geometry().dim())
model=geopart.stokes.StokesModel(eta=mu,f=f)

#Form and solve Stokes
A, b = dolfin.PETScMatrix(), dolfin.PETScVector()
element_cls.solve_stokes(W, U, (A, b), weak_bcs, model)
uh, ph = element_cls.get_velocity(U), element_cls.get_pressure(U)

#Output solution p,u to paraview
dolfin.XDMFFile("pressure.xdmf").write_checkpoint(ph, "p")
dolfin.XDMFFile("velocity.xdmf").write_checkpoint(uh, "u")

flux = [dolfin.assemble(dolfin.dot(uh, n)*ds(i)) for i in range(len(mark))]

if comm.Get_rank() == 0:
    for key, value in mark.items():
        print("Flux_%s= %.15lf"%(key,flux[value]))

You can find example code and mesh file

2

I’m a little busy this week but I will take a look as soon as I can.

2 Likes
3

I just did one more testing that using taylor-hood in fenics instead of HDG. I still have the same issue.

Here is the MWE for the taylor-hood @ * hdg_test/2d/cg_test.py

So I’m guessing I may didn’t set up the material properly. But I did following test which looks good to me, is this a bug in fenics?

coefficient
coefficient757×357 217 KB 

import numpy as np
import dolfin
from dolfin import *

#read mesh and boundaries from file
mesh = Mesh()
hdf = HDF5File(mesh.mpi_comm(), "mesh_boundaries.h5", "r")
hdf.read(mesh, "/mesh", False)
boundaries = MeshFunction('size_t', mesh, mesh.topology().dim() - 1)
hdf.read(boundaries, "/boundaries")
hdf.close()

#read viscosity coefficient from file
mu_ele = FunctionSpace(mesh, "DG", 0)
mu = Function(mu_ele)
hdf = HDF5File(mesh.mpi_comm(), "mesh_coeffs.h5", "r")
hdf.read(mu, "/mu")
hdf.close()

#Visulization !!!only for 2D
import matplotlib.tri as mtri
import matplotlib.pyplot as plt
plt.figure(figsize=(10,4))

nodes=mesh.coordinates()
cells=mesh.cells()

#plot mesh with coefficient
triang = mtri.Triangulation(nodes[:,0], nodes[:,1], triangles=cells)
plt.tripcolor(triang, facecolors=mu.vector()[:],linewidth=0.5,edgecolors='k',cmap='rainbow')
plt.colorbar()

plt.axis('equal')
plt.savefig('coefficient.png', bbox_inches='tight',dpi=100)

You can find updated example code and mesh file here:

4

The effect you’re seeing looks like an inherent component of the physical model.

Switching the BCs of the 3D problem to be equivalent with the 2D problem yields a consistent result.

I made the following changes to the BCs in your Taylor-Hood implementation to yield the images below

bc0 = DirichletBC(W.sub(0), noslip, "(near(x[1], 10.0) or near(x[1], 0.0)) and on_boundary")
bc1 = DirichletBC(W.sub(0).sub(2), Constant(0.0), "(near(x[2], 0.0) or near(x[2], 4.0)) and on_boundary")
bcs = [bc0, bc1]

image
image707×497 47.6 KB

image

My recommendation would be to carefully consider what the boundary conditions should be in the case of the 2D model.

With respect to your initial experiments: In the 2D model, the incoming flow introduces a much greater volume flux than the 3D model with no-slip BCs everywhere.

As to why this result is physically valid, my investigation continues… Even with the higher volume flux, it is counter intuitive for the left side (very) high viscosity region to affect the flow as seen here.

1 Like
5

Thanks for checking this. It seems really an anti-intuitive physical effect.

I used COMSOL to run a simulation like this:

image

It gives us similar results we got in fenics.

BTW, the implementation of slip boundary conditions looks quite interesting.

How can I prescribe slip boundary conditions in GeoPart/Leopart?