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))
# Symmetric gradient
def epsilon(u):
    return sym(nabla_grad(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

# Variational problem to compute the ls normal as n=grad(ls)/|grad(ls)|
Vn = VectorFunctionSpace(mesh, "CG", 1)
Fn = FunctionSpace(mesh, "CG", 2)
n_K = TrialFunction(Vn)
vn = TestFunction(Vn)

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

normal = sqrt(dot(grad(ls),grad(ls)))
an = normal*dot(n_K,vn)*dx
Ln = dot(grad(ls),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
dxq = dc(0, metadata={"num_cells": 1})
dsq = dc(1, metadata={"num_cells": 1})
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 
jump_grad_u = jump(grad(u_), norm)
jump_grad_v = jump(grad(v), norm)

jump_grad_p = dot(grad(p_)('+'),norm('+'))-dot(grad(p_)('-'),norm('-'))
jump_grad_q = dot(grad(q)('+'),norm('+'))-dot(grad(q)('-'),norm('-'))

jump_grad_temp = dot(grad(theta_)('+'),norm('+'))-dot(grad(theta_)('-'),norm('-'))
jump_grad_theta = dot(grad(theta)('+'),norm('+'))-dot(grad(theta)('-'),norm('-'))

# Navier Stokes momentum equation
Eq1 = dot(u_-u_n, v)/dt*dx\
    + dot(dot(grad(u_), u_n), v)*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 \
    + 1/Pr*dot(nabla_grad(theta_), nabla_grad(theta))*dx \
    - dot(theta_*u_n, nabla_grad(theta))*dx\
    + dot(flux,norm)*theta*ds(1)

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

# Stabilization terms
Ghost_Penalty_Temperature = avg(gamma_ghost_temperature)*avg(h)*avg(k)*jump_grad_temp*jump_grad_theta*(dS(1)+dS(3)) # Ghost penalty for temperature
Ghost_Penalty_Pressure = avg(gamma_ghost_pressure)*avg(1/mu(theta_n))*avg(h**3)*dot(jump_grad_p,jump_grad_q)*(dS(1)+dS(3)+dS(0)) # Ghost penalty and CIP for pressure
Ghost_Penalty_Velocity = avg(gamma_ghost_velocity)*avg(mu(theta_n))*avg(h)*dot(jump_grad_u,jump_grad_v)*(dS(1)+dS(3)) # Ghost penalty for velocity

# 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)
    quadrature_cut_cells_0 = Quadrature(cut_cells_0.type().cell_type(),
                                    cut_cells_0.geometry().dim(),
                                    1)

    interface_cells = mesh_cutter.interface_cells()
    quadrature_interface = Quadrature(interface_cells.type().cell_type(),
                                  interface_cells.geometry().dim(),
                                  2);

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

    a1_cut.set_quadrature(1, quadrature_interface);
    a1_cut.set_cut_mesh(1, interface_cells);

    L1_cut.set_quadrature(0, quadrature_cut_cells_0);
    L1_cut.set_cut_mesh(0, cut_cells_0);

    L1_cut.set_quadrature(1, quadrature_interface);
    L1_cut.set_cut_mesh(1, interface_cells);
	
    a2_cut.set_quadrature(0, quadrature_cut_cells_0);
    a2_cut.set_cut_mesh(0, cut_cells_0);

    a2_cut.set_quadrature(1, quadrature_interface);
    a2_cut.set_cut_mesh(1, interface_cells);

    L2_cut.set_quadrature(0, quadrature_cut_cells_0);
    L2_cut.set_cut_mesh(0, cut_cells_0);

    L2_cut.set_quadrature(1, quadrature_interface);
    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())
    A1.apply("add")
	
    assembler.apply_constraints_cell_domains(A2, a2_cut.form(), W_T, 2, mesh_cutter.domain_marker())
    A2.apply("add")
 
    # 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  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »