Dear all,
I’m rather new to FeniCS and to FEM so have a hard time currently. During solving of the first problem of the non-linear model. I can not understand how you can add boundary conditions (Neuman and Dirichlet). All information on the web did not help me up to now.
I try to solve the time-dependent non-linear PDE with different boundary conditions at different domain positions (Dirichlet and Neumann).
import meshio
import gmsh
import pygmsh
from dolfinx import log, default_scalar_type
from dolfinx.fem.petsc import NonlinearProblem, LinearProblem, assemble_matrix, create_vector, assemble_vector, apply_lifting, set_bc
from dolfinx.nls.petsc import NewtonSolver
from dolfinx.io import VTXWriter, XDMFFile, gmshio
from dolfinx import nls
import numpy as np
import ufl
import dolfinx
print(f"DOLFINx version: {dolfinx.__version__} based on GIT commit: {dolfinx.git_commit_hash} of https://github.com/FEniCS/dolfinx/")
from mpi4py import MPI
from dolfinx import fem, mesh, plot
from petsc4py.PETSc import ScalarType
from petsc4py import PETSc #solvers
from dolfinx.cpp.mesh import create_meshtags
from dolfinx.fem import (Constant, dirichletbc, Function, FunctionSpace,
form, assemble_scalar, locate_dofs_geometrical, locate_dofs_topological, assemble, DirichletBC)
from dolfinx.mesh import create_unit_square, locate_entities, locate_entities_boundary
from ufl import (SpatialCoordinate, TestFunction, TrialFunction,
dx, grad, inner, Measure, dot, jump, FacetNormal, avg, FiniteElement, MixedElement)
from ufl import (FacetNormal, FiniteElement, Identity, Measure, TestFunction, TestFunctions, TrialFunction, VectorElement, MixedElement, as_vector, div, dot, ds, dx, inner, lhs, grad, nabla_grad, rhs, sym, curl, derivative, split)
from dolfinx.plot import vtk_mesh
## We define our problem specific parameters
dt = Constant(mesh, PETSc.ScalarType(1./1600))
nu = Constant(mesh, PETSc.ScalarType(1/500)) # nu = 1/Re
relax = Constant(mesh, PETSc.ScalarType(1e-6)) # relaxation constant
U_0 = 2 # longitunal/lenghtwise velocity at coordinate (x,y)=(0,H/2)
# Define constants
L = 1
a = 3.0
rho = 1.0
kappa= 0.25
# Define the parameters and functions involved in the equations
K_x, K_y = 1.0, 1.0 # Replace with actual values
K = [K_x, K_y]
a = 3
gamma = 0.7
delta = 0.0025
lambda_ = 0.5
T = 1.0
nu=0.25
num_steps = T / dt
# Boundary condtions ---------------------------------------------------------
u_plate = 25
u0 = 23
u0_1 = 21
u0_2 = 14
# Time ----------------------------------------------------------------------
TotalTime = 1000.0
t0=time.time()
num_steps = 50 # number of time steps
dt = TotalTime / num_steps # time step size
t = 0 # current time
# Mesh
with XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "r") as xdmf:
mesh = xdmf.read_mesh(name="Grid")
ct = xdmf.read_meshtags(mesh, name="Grid")
mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim - 1)
with XDMFFile(MPI.COMM_WORLD, "facet_mesh.xdmf", "r") as xdmf:
ft = xdmf.read_meshtags(mesh, name="Grid")
ds = Measure('ds', domain=mesh, subdomain_data=ft)
# Define the variational problem
V = FunctionSpace(mesh, ("CG", 1))
Q = FunctionSpace(mesh, ("DG", 0))
u_elem = VectorElement('CG', mesh.ufl_cell(), 2) # displacement(2D)
v_elem= FiniteElement('CG', mesh.ufl_cell(), 1) # pore water pressure(scalar)
V = FunctionSpace(mesh, MixedElement([u_elem, v_elem])) #Define vector function space for [u,p]
u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)
f = Constant(mesh, PETSc.ScalarType(0))
n = FacetNormal(mesh)
#Trial function
xt = Function(V)
ut, vt = split(xt)
#Define function for step
x_n = Function(V)
print(x_n.vector.array)
u_n, v_n = split(x_n)
V0, submap = V.sub(0).collapse()
#Define boundary conditions and Neuman conditon in the above domain
# in the circle
Timp_facets = ft.find(FACE_MARK)
Timp_dofs = locate_dofs_topological(V, mesh.topology.dim-1, Timp_facets)
bc = dirichletbc(default_scalar_type(u_plate), Timp_dofs, V)
def boundary_D(x):
return np.logical_or(np.isclose(x[0], 0), np.isclose(x[0], 1))
dofs_D = locate_dofs_geometrical(V, boundary_D)
u_bc = Function(V)
u_bc.interpolate(u_0)
bc = dirichletbc(u_bc, Timp_dofs)
Res = dot(u,v)*dx
Res += - dt*dot(K*grad(u), grad(v))*dx + dot(K*F1(u)*grad(u),v)*dx*dt
Res += - inner(p,q)*dx - a*inner(p, q)*dx*dt + gamma*inner((phiF(u, lambda_) - p), q)*dx*dt
Jac = derivative(Res, xt) # jacobian
Prob = NonlinearProblem(Res, xt, bc, Jac)
solver = nls.petsc.NewtonSolver(MPI.COMM_WORLD, Prob)
# Set nonlinear solver parameters
solver.rtol = 1e-7 #relative_tolerance
solver.atol = 1e-7 #absolute_tolerance
solver.max_it = 100 #maximum_iterations
solver.report = True
# Solve system & output results
# ----------------------------------------------------------------
# time stepping
t_i = 0.0
t_f = 1.0
Nsteps = 100
time = np.linspace(t_i, t_f, Nsteps+1)
# output file
#xdmf = XDMFFile(domain.comm, "terzaghi.xdmf", "w")
#xdmf.write_mesh(domain)
with XDMFFile(MPI.COMM_WORLD, "ux.xdmf", "w") as xdmf1:
xdmf1.write_mesh(mesh)
with XDMFFile(MPI.COMM_WORLD, "uy.xdmf", "w") as xdmf2:
xdmf2.write_mesh(mesh)
with XDMFFile(MPI.COMM_WORLD, "v.xdmf", "w") as xdmf3:
xdmf3.write_mesh(mesh)
t = 0
u, v = xt.split()
u.name = "u"
v.name = "v"
dx = u.sub(0).collapse()
dy = u.sub(1).collapse()
xdmf1.write_function(dx, t)
xdmf2.write_function(dy, t)
xdmf3.write_function(v, t)
for (i, dt) in enumerate(np.diff(time)):
t = time[i+1]
print('-----------------------------------------')
print('>> t =', t, '[sec]')
# solve system
solver.solve(xt)
x_n.vector.array[:] = xt.vector.array[:]
# data for visualization
u, v = xt.split()
dx = u.sub(0).collapse()
dy = u.sub(1).collapse()
xdmf1.write_function(dx, t)
xdmf2.write_function(dy, t)
xdmf3.write_function(v, t)
#xdmf_file.write(poro, t)
xdmf1.close()
xdmf2.close()
xdmf3.close()