Hi,
I am trying to solve the Stokes flow (having same form as the elastic equilibrium equation) on a square grid with Periodic Boundary Condition (PBC) along both the x and y directions. The code is generating somewhat sensible outputs but I am having large errors at the corners.
Could you please suggest me how I can get rid of these larger errors at the corners of the square domain. I have copied my code below where I have implemented the PBCs on both directions. Is this the correct way of implementing PBC? At the end of this post I have attached the plots of the solution vector components.
Thanks!
### Linear elasticity example FEniCSx tutorial ###
######################################################################
### importing packages and defining parameters starts ###
######################################################################
from __future__ import annotations
from pathlib import Path
from typing import Union
from mpi4py import MPI
from petsc4py import PETSc
import dolfinx.fem as fem
import numpy as np
import scipy.sparse.linalg
from dolfinx import default_scalar_type
from dolfinx.common import Timer, TimingType, list_timings
from dolfinx.io import XDMFFile
from dolfinx.mesh import create_unit_square, locate_entities_boundary
from ufl import (SpatialCoordinate, TestFunction, TrialFunction, as_vector, dx,
exp, grad, inner, pi, sin, cos)
import dolfinx_mpc.utils
from dolfinx_mpc import LinearProblem, MultiPointConstraint
import pandas as pd
import pyvista
from dolfinx import mesh, fem, plot, io
from petsc4py import PETSc
default_scalar_type = PETSc.ScalarType
import ufl
#L = 1
#W = 0.2
mu = 0.001 ### 0.001, 1
#rho = 1
#delta = W / L
#gamma = 0.4 * delta**2
#beta = 1.25
#lambda_ = beta
lambda_ = 1
print("Value of mu: {}, and lambda: {}".format(mu,lambda_))
#g = gamma
######################################################################
### importing packages and defining parameters ends ###
######################################################################
######################################################################
### defining mesh starts ###
######################################################################
# domain = mesh.create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([L, W, W])],
# [20, 6, 6], cell_type=mesh.CellType.hexahedron)
# nx, ny = 50, 50
# domain = mesh.create_rectangle(MPI.COMM_WORLD, [np.array([-2, -2]), np.array([2, 2])],
# [nx, ny], mesh.CellType.triangle)
# Get PETSc int and scalar types
complex_mode = True if np.dtype(default_scalar_type).kind == 'c' else False
# Create mesh and finite element
NX = 50
NY = 50
mesh = create_unit_square(MPI.COMM_WORLD, NX, NY)
#V = fem.functionspace(mesh, ("Lagrange", 1, (mesh.geometry.dim, )))
V = fem.functionspace(mesh, ("Lagrange", 2, (mesh.geometry.dim, )))
tol = 250 * np.finfo(default_scalar_type).resolution
# Listing space points ###
dof_coordinates = V.tabulate_dof_coordinates()
df = pd.DataFrame(dof_coordinates)
df.to_csv('outdata/planenodesPeriodic.csv', index = None, header=False)
######################################################################
### defining mesh ends ###
######################################################################
######################################################################
### boundary conditions starts ###
######################################################################
# def dirichletboundary(x):
# return np.logical_or(np.isclose(x[1], 0, atol=tol), np.isclose(x[1], 1, atol=tol))
# # Create Dirichlet boundary condition
# facets = locate_entities_boundary(mesh, 1, dirichletboundary)
# topological_dofs = fem.locate_dofs_topological(V, 1, facets)
# zero = np.array([0, 0], dtype=default_scalar_type)
# bc = fem.dirichletbc(zero, topological_dofs, V)
# bcs = [bc]
bcs = []
# Sub domain for Periodic boundary condition
def periodic_boundary0(x):
return np.isclose(x[0], 1, atol=tol)
def periodic_boundary1(x):
return np.isclose(x[1], 1, atol=tol)
# def periodic_relation(x):
# out_x = np.zeros_like(x)
# out_x[0] = 1 - x[0]
# out_x[1] = x[1]
# out_x[2] = x[2]
# return out_x
def periodic_relation0(x):
out_x = np.zeros_like(x)
out_x[0] = 1 - x[0]
out_x[1] = x[1]
out_x[2] = x[2]
return out_x
def periodic_relation1(x):
out_x = np.zeros_like(x)
out_x[0] = x[0]
out_x[1] = 1 - x[1]
out_x[2] = x[2]
return out_x
with Timer("~PERIODIC: Initialize MPC"):
mpc = MultiPointConstraint(V)
mpc.create_periodic_constraint_geometrical(V, periodic_boundary0, periodic_relation0, bcs)
mpc.create_periodic_constraint_geometrical(V, periodic_boundary1, periodic_relation1, bcs)
#mpc.create_periodic_constraint_geometrical(V, periodic_boundary, periodic_relation, bcs)
mpc.finalize()
### traction free BC starts on other faces ###
T = fem.Constant(mesh, default_scalar_type((0, 0)))
### traction free BC ends on other faces ###
### defining surface integral over the boundary starts ###
ds = ufl.Measure("ds", domain=mesh)
### defining surface integral over the boundary ends ###
######################################################################
### boundary conditions end ###
######################################################################
######################################################################
### variational formulation starts ###
######################################################################
def epsilon(u):
return ufl.sym(ufl.grad(u)) # Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)
def sigma(u):
return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2 * mu * epsilon(u)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
x = SpatialCoordinate(mesh)
#f = as_vector((sin(2 * pi * x[0]), 0))
f = as_vector((1000 * cos(2 * pi * x[0] + 2 * pi * x[1]), 1000 * cos(2 * pi * x[0] + 2 * pi * x[1])))
#f = fem.Constant(mesh, default_scalar_type((0, 0)))
#f = as_vector((x[0] * sin(5.0 * pi * x[1]) + 1.0 * exp(-(dx_ * dx_ + dy_ * dy_) / 0.02), 0.3 * x[1]))
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
L = ufl.dot(f, v) * ufl.dx + ufl.dot(T, v) * ds
######################################################################
### variational formulation ends ###
######################################################################
######################################################################
### solving variational problem starts ###
######################################################################
# problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
# #problem = LinearProblem(a, L, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
# uh = problem.solve()
# Setup MPC system
with Timer("~PERIODIC: Initialize varitional problem"):
problem = LinearProblem(a, L, mpc, bcs=bcs)
solver = problem.solver
# Give PETSc solver options a unique prefix
solver_prefix = "dolfinx_mpc_solve_{}".format(id(solver))
solver.setOptionsPrefix(solver_prefix)
petsc_options: dict[str, Union[str, int, float]]
if complex_mode or default_scalar_type == np.float32:
petsc_options = {"ksp_type": "preonly", "pc_type": "lu"}
else:
petsc_options = {"ksp_type": "cg", "ksp_rtol": 1e-6, "pc_type": "hypre", "pc_hypre_type": "boomeramg",
"pc_hypre_boomeramg_max_iter": 1, "pc_hypre_boomeramg_cycle_type": "v" # ,
# "pc_hypre_boomeramg_print_statistics": 1
}
# Set PETSc options
opts = PETSc.Options() # type: ignore
opts.prefixPush(solver_prefix)
if petsc_options is not None:
for k, v in petsc_options.items():
opts[k] = v
opts.prefixPop()
solver.setFromOptions()
with Timer("~PERIODIC: Assemble and solve MPC problem"):
uh = problem.solve()
# solver.view()
it = solver.getIterationNumber()
print("Constrained solver iterations {0:d}".format(it))
######################################################################
### solving variational problem ends ###
######################################################################
usol = uh.x.array
numpoints = int(len(usol)/2)
usol = np.reshape(usol, (numpoints,2))
#print(usol.tolist())
df = pd.DataFrame(usol)
df.to_csv('outdata/planedisplacementsPeriodic.csv', index = None, header=False)
print("Number of nodes {}".format(len(dof_coordinates)))
print("Number of sol points {}".format(len(usol)))
