### Segmentation fault when applying BC

78
views
-2
7 weeks ago by

Dear all,

I am solving a phase change problem with special elements. I must detail the assembly of the system because of a special quadrature rule and assemble at each operation in order to update the system correctly (I tested it with classic finite element). However, I have an issue because I have a segmentation error at the second iteration when I apply the boundary conditions.

Does someone know why I have such an issue ? Thank you in advance. Below are my code and the error message.

from __future__ import print_function
from fenics import *
from dolfin import *
from cutfem import *

# Have the compiler generate code for evaluating derivatives
parameters["form_compiler"]["no-evaluate_basis_derivatives"] = False

# MESH AND RELATED PROPERTIES

N = 101
# Larger mesh than the domain we want to study
mesh = RectangleMesh(Point(0.0,0.0),Point(1.0,1.1),N,N)
# Mesh Properties
norm = FacetNormal(mesh)
h = CellSize(mesh)

# Level set expression for the interface
level_set = Expression("fabs(x[0]-0.5)>=fabs(x[1]-0.5) ? fabs(x[0]-0.5)-0.5 : fabs(x[1]-0.5)-0.5 ", degree = 2)

# CONSTANTS AND FUNCTIONS DEFINING THE PROBLEM

# Dimensionless numbers
Pr = Constant(1.0)  # Prandtl
Ste = Constant(1.0) # Stefan
Ra = Constant(10000.0)  # Rayleigh
DeltaT = Constant(0.001)
k = Constant(1.0)

# Stabilization parameters
gamma = Constant(1E-7)
gamma_ghost_pressure = Constant(0.005)
gamma_ghost_temperature = Constant(0.1)
gamma_ghost_velocity = Constant(0.1)

# Time discretization
t0 = 0.0
tf = 0.5
dt = 0.01
num_steps = int((tf-t0)/dt) # number of time steps
t = t0

# Phase change function
def phi(temp):
Tr = 0.0
r = 0.1
return(0.5+0.5*tanh((Tr-temp)/r))
# Viscosity depending on the temperature
def mu(temp):
muL = 1
muS = 100000
return(muL + (muS - muL)*phi(temp))
def epsilon(u):

# Local application of the laser beam on the top wall
width = Constant(0.30)
flux = Expression(("0.","(x[0]>(0.5-e/2) && x[0]<(0.5+e/2)) ? 1000*(x[0]-(0.5-e/2))*(x[0]-(0.5+e/2)) : 0."), degree=2, e=width)

# NORMAL OF THE LEVEL SET FUNCTION

Vn = VectorFunctionSpace(mesh, "CG", 1)
Fn = FunctionSpace(mesh, "CG", 2)
n_K = TrialFunction(Vn)
vn = TestFunction(Vn)

ls = Function(Fn)
ls.interpolate(level_set)

an = normal*dot(n_K,vn)*dx

An = assemble(an)
bn = assemble(Ln)
n_K = Function(Vn)
solve(An,n_K.vector(),bn)

file = File("./LaserBeamCut/Mesh/Level_Set_Normal.pvd")
file << n_K

# MESH CELLS MARKERS

# Compute Mesh - Levelset intersection and corresponding marker
mesh_cutter = MeshCutter(mesh,level_set)
interior_facet_marker = CutFEMTools_compute_interface_ghost_penalty_facets(mesh_cutter.domain_marker())

# Output the results of the mesh-level set intersection
file = File("./LaserBeamCut/Mesh/Domain_Marker.pvd")
file << mesh_cutter.domain_marker()

file = File("./LaserBeamCut/Mesh/Interface_Cells.pvd")
file << mesh_cutter.interface_cells()

file = File("./LaserBeamCut/Mesh/Cut_Cells.pvd")
file << mesh_cutter.cut_cells(0)

file = File("./LaserBeamCut/Mesh/Interior_Facets.pvd")
file << interior_facet_marker

# BOUNDARIES

# No slip and uniform temperature walls
walls = 'near(x[0], 0.0) || near(x[1], 0.0) || near(x[0], 1.0)'

# INTEGRATION MEASUREMENTS

# Define new measures associated with the interior domains and facets
dx = Measure("dx")[mesh_cutter.domain_marker()]
dS_ghost = Measure("dS")[interior_facet_marker]
# Custom integrals to integrate over cut cells
dxc = dx(0) + dxq

# MIXED FINITE ELEMENTS AND FUNCTION SPACE

V = VectorFunctionSpace(mesh, 'P', 1) # velocity
Q = FunctionSpace(mesh, 'P', 1) # pressure
W = MixedFunctionSpace([V, Q])
W_T = FunctionSpace(mesh, 'P', 1) # temperature

# INITIAL CONDITIONS

w_n = interpolate(Expression(("0.", "0.", "0."), degree = 1), W)
u_n, p_n = split(w_n)

theta_n = interpolate(Expression("-1.0", degree = 1), W_T)

# DIRICHLET BOUNDARY CONDITIONS

u_D = Constant((0.0, 0.0))
theta_D = Constant(-1.0)
bc_velocity = DirichletBC(W.sub(0), u_D, walls)
bc_temperature = DirichletBC(W_T, theta_D, walls)

# VARIATIONAL FORMULATION

# Test functions
w = TestFunction(W)
(v, q) = split(w)
theta = TestFunction(W_T)

# Trial functions
w_ = TrialFunction(W)
(u_, p_) = split(w_)
theta_ = TrialFunction(W_T)

# Define the buoyancy coupling term
fB = Ra/Pr*theta_n*Constant((0.0,1.0))

# Define the jumps of the gradients

# Navier Stokes momentum equation
Eq1 = dot(u_-u_n, v)/dt*dx\
+ 2*mu(theta_n)*inner(epsilon(u_), epsilon(v))*dx \
- div(v)*p_*dx \
- dot(fB, v)*dx

# Mass equation
Eq2 = q*div(u_)*dx

# Enthalpy equation
Eq3 = (theta_-theta_n)*theta/dt*dx \
- 1/Ste*(phi(theta_n+DeltaT)-phi(theta_n))*theta/dt*dx \
+ dot(flux,norm)*theta*ds(1)

# Pressure stabilization
Stab_pressure = gamma*p_*q*dxc

# Stabilization terms

# Nitsche term to stabilize Neumann BC may be added but not necessary

# Monolithic system
F1 = Eq1 + Eq2 \
+ Stab_pressure + Ghost_Penalty_Pressure + Ghost_Penalty_Velocity
F2 = Eq3 + Ghost_Penalty_Temperature

# SAVE SOLUTION TO VTK FILES

file_velocity = File('./LaserBeamCut/Velocity/velocity.pvd')
file_pressure = File('./LaserBeamCut/Pressure/pressure.pvd')
file_temperature = File('./LaserBeamCut/Temperature/temperature.pvd')
file_concentration = File('./LaserBeamCut/Concentration/concentration.pvd')

file_velocity_phys = File('./LaserBeamCut/Velocity/velocity_phys.pvd')
file_pressure_phys = File('./LaserBeamCut/Pressure/pressure_phys.pvd')
file_temperature_phys = File('./LaserBeamCut/Temperature/temperature_phys.pvd')

# Time-progress bar
progress_bar = Progress('Time-stepping')
set_log_level(PROGRESS)

# RESOLUTION FOR ALL TIME STEPS
w_ = Function(W)
theta_ = Function(W_T)

for n in range(num_steps):

# Update current time
t += dt

# Define bilinear and linear forms
a1 = lhs(F1)
L1 = rhs(F1)

a2 = lhs(F2)
L2 = rhs(F2)

# Create a cutform
form_a1 = create_dolfin_form(a1)
form_L1 = create_dolfin_form(L1)

form_a2 = create_dolfin_form(a2)
form_L2 = create_dolfin_form(L2)

a1_cut = CutForm(form_a1)
L1_cut = CutForm(form_L1)

a2_cut = CutForm(form_a2)
L2_cut = CutForm(form_L2)

# Define quadrature rule for cut cells
cut_cells_0 = mesh_cutter.cut_cells(0)
cut_cells_0.geometry().dim(),
1)

interface_cells = mesh_cutter.interface_cells()
interface_cells.geometry().dim(),
2);

# Set quadrature rules for cut cells
a1_cut.set_cut_mesh(0, cut_cells_0);

a1_cut.set_cut_mesh(1, interface_cells);

L1_cut.set_cut_mesh(0, cut_cells_0);

L1_cut.set_cut_mesh(1, interface_cells);

a2_cut.set_cut_mesh(0, cut_cells_0);

a2_cut.set_cut_mesh(1, interface_cells);

L2_cut.set_cut_mesh(0, cut_cells_0);

L2_cut.set_cut_mesh(1, interface_cells);

# Assemble system
A1 = cutfem_assemble(a1_cut)
b1 = cutfem_assemble(L1_cut)

A2 = cutfem_assemble(a2_cut)
b2 = cutfem_assemble(L2_cut)

# Apply Dirichlet boundary conditions
bc_velocity.apply(A1,b1)
bc_temperature.apply(A2, b2)

# Deactivate dofs outside of domain
assembler = CutFEMAssembler()
assembler.apply_constraints_cell_domains(A1, a1_cut.form(), W, 2, mesh_cutter.domain_marker())

assembler.apply_constraints_cell_domains(A2, a2_cut.form(), W_T, 2, mesh_cutter.domain_marker())

# Solve problem
solve(A1, w_.vector(), b1)
u_, p_ = w_.split()
w_n.vector()[:] = w_.vector()
solve(A2, theta_.vector(), b2)
theta_n.vector()[:] = theta_.vector()

# Save solution
file_velocity << (u_, t)
file_pressure << (p_, t)
file_temperature << (theta_, t)

# Physical solution
phys_mesh = CutFEMTools_physical_domain_mesh(mesh, mesh_cutter.cut_cells(0),      					mesh_cutter.domain_marker(), 0)

Vphys = VectorFunctionSpace(phys_mesh, 'P', 1)
Qphys = FunctionSpace(phys_mesh, 'P', 1)
W = MixedFunctionSpace([Vphys, Qphys])

w_phys = Function(W)
w_phys.interpolate(w_)
u_phys, p_phys = w_phys.split()

theta_phys = Function(W_T)
theta_phys.interpolate(theta_)

# Save physical solution
file_velocity_phys << (u_phys, t)
file_pressure_phys << (p_phys, t)
file_temperature_phys << (theta_phys, t)

# Post-treatment for the phase field distribution
C = FunctionSpace(mesh, 'P', 1)
concentration = Function(C)
c = TestFunction(C)
Fc = dot((concentration-phi(theta_n)),c)*dxc
solve(Fc == 0, concentration)
Cphys = FunctionSpace(phys_mesh, 'P', 1)
concentration_phys = Function(Cphys)
concentration_phys.interpolate(concentration)
file_concentration << (concentration_phys, t)

# Update the progress bar
progress_bar.update(t/tf)

Applying boundary conditions to linear system.
Applying boundary conditions to linear system.
CompositeFormAssembler: Inserted 2760 ghost_values of size 1.000000e+00 in matrix
CompositeFormAssembler: Inserted 920 ghost_values of size 1.000000e+00 in matrix
Solving linear system of size 31212 x 31212 (PETSc LU solver, petsc).
Solving linear system of size 10404 x 10404 (PETSc LU solver, petsc).
num_cells in phys_mesh=18960
num_vertices in phys_mesh=12520
Still alive!!!
Computed bounding box tree with 40803 nodes for 20402 entities.
No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F'.
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 1.012e-02 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Solving linear system of size 10404 x 10404 (PETSc LU solver, petsc).
Newton iteration 1: r (abs) = 4.032e-17 (tol = 1.000e-10) r (rel) = 3.982e-15 (tol = 1.000e-09)
Newton solver finished in 1 iterations and 1 linear solver iterations.
Applying boundary conditions to linear system.
Applying boundary conditions to linear system.
Segmentation fault (core dumped)

Community: FEniCS Project
1
This doesn't seem like a minimal working example. If you strip your code down to only the necessary part to end up at your problem, it's easier for people to help you. Also, more often than not, I found my problem myself by trying to create a MWE.
written 7 weeks ago by alexdiem