Hi, I am a student trying to solve a Stokes system that has a liquid phase bubble within a solid phase, which I represent with two separate viscosities. I am following the example shown here at the bottom example with varying permeabilities.
My code is here:
# Set domain geometries L = 400 H = 200 resolution = 32 domain = Rectangle(Point(0.0,0.0), Point(L,H)) liquid = Circle(Point(100, 0.5*H), 50) solid = domain - liquid domain.set_subdomain(1, liquid) domain.set_subdomain(2, solid) # Export mesh as xml file, read the xml file domain_mesh = generate_mesh(domain, resolution) File("meshes/twomaterials/mesh.xml") << domain_mesh mesh = Mesh("meshes/twomaterials/mesh.xml") # Define subdomains ---------------------------- # Sub domain for no-slip (mark whole boundary, inflow and outflow will overwrite) class Noslip(SubDomain): def inside(self, x, on_boundary): return on_boundary # Sub domain for inflow (left) class Inflow(SubDomain): def inside(self, x, on_boundary): return x < DOLFIN_EPS and on_boundary # Sub domain for outflow (right) class Outflow(SubDomain): def inside(self, x, on_boundary): return x > L - 1e3*DOLFIN_EPS and on_boundary sub_domains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1) sub_domains.set_all(3) noslip = Noslip() noslip.mark(sub_domains, 0) inflow = Inflow() inflow.mark(sub_domains, 1) outflow = Outflow() outflow.mark(sub_domains, 2) file = File("meshes/twomaterials/subdomains.xml") file << sub_domains # defining problem space ----------------- # Taylor-Hood function space VE = VectorElement("DG", mesh.ufl_cell(), 2) # previously was CG, not DG QE = FiniteElement("DG", mesh.ufl_cell(), 1) WE = VE * QE W = FunctionSpace(mesh, WE) V = FunctionSpace(mesh, VE) Q = FunctionSpace(mesh, QE) # boundary conditions----------------------- # No-slip boundary condition for the liquid pocket noslip = Constant((0, 0)) bc0 = DirichletBC(W.sub(0), noslip, sub_domains, 0) # subdomain marked 0 (noslip) # Inflow boundary condition for velocity YMAX = H inflow = Expression(("4*(x*(YMAX-x))/(YMAX*YMAX)", "0."), YMAX=YMAX, element = V.ufl_element()) bc1 = DirichletBC(W.sub(0), inflow, sub_domains, 1) # subdomain marked 1 (inflow) # Boundary condition for pressure at outflow # x0 = 0 zero = Constant(0) bc2 = DirichletBC(W.sub(1), zero, sub_domains, 2) # subdomain marked 2 (outflow) # Collect boundary conditions bcs = [bc0, bc1, bc2] # "Define subdomain markers and integration measure" markers = MeshFunction('size_t', mesh, 2, mesh.domains()) dx = Measure('dx', domain=mesh, subdomain_data=markers) # Define an Expression for viscosity class Viscosity(UserExpression): def __init__(self, markers, **kwargs): super().__init__(**kwargs) self.markers = markers def eval_cell(self, values, x, cell): if self.markers[cell.index]==1: # liquid values==1 # viscosity in liquid is 1 else: values==3e15 # solid viscosity much higher def value_shape(self): # from https://fenicsproject.discourse.group/t/problems-interpolating-a-userexpression-and-plotting-it/1303/2 #return (1,) return () eta = Viscosity(markers) # Define variational problem I = Identity(2) def epsilon(u): # strain rate tensor return 0.5*(grad(u) + grad(u).T) def sigma(u,p): # Cauchy stress tensor return 2*eta*epsilon(u) - p*I # from https://fenicsproject.org/pub/tutorial/html/._ftut1008.html # Define trial and test functions u, p = TrialFunctions(W) v, q = TestFunctions(W) f = Constant((0, 0)) # set f, body force, to 0 a = inner(sigma(u,p), grad(v))*dx + q*div(u)*dx L = inner(f, v)*dx # Compute solution w = Function(W) solve(a == L, w, bcs) (u, p) = w.split(True) plot(u) plt.show() plot(p) plt.show()
The error I keep getting is shown below. Since there are no values being computed, my plots also come out blank. I have spent many hours trying to figure out where it is coming from but I am not even sure what part of the code is causing this malfunction Any help would be greatly appreciated!
/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/matplotlib/quiver.py:715: RuntimeWarning: divide by zero encountered in double_scalars length = a * (widthu_per_lenu / (self.scale * self.width)) /anaconda3/envs/fenicsproject/lib/python3.8/site-packages/matplotlib/quiver.py:715: RuntimeWarning: invalid value encountered in multiply length = a * (widthu_per_lenu / (self.scale * self.width)) /anaconda3/envs/fenicsproject/lib/python3.8/site-packages/matplotlib/quiver.py:767: RuntimeWarning: invalid value encountered in less short = np.repeat(length < minsh, 8, axis=1) /anaconda3/envs/fenicsproject/lib/python3.8/site-packages/matplotlib/quiver.py:780: RuntimeWarning: invalid value encountered in less tooshort = length < self.minlength