import os
from pathlib import Path
from mpi4py import MPI
from petsc4py import PETSc
import dolfinx
import numpy as np
import ngsPETSc.utils.fenicsx as ngfx
from packaging.version import Version
from netgen.geom2d import SplineGeometry
from ufl import (
FacetNormal,
Identity,
TestFunction,
TrialFunction,
div,
dot,
ds,
dx,
inner,
lhs,
nabla_grad,
rhs,
sym,
TestFunctions,
split,
variable,
diff,
inner,
grad,
as_vector,
sqrt
)
from basix.ufl import element, mixed_element
from dolfinx import default_real_type, log, plot,mesh,geometry
from dolfinx.fem import(
Function,
functionspace,
Expression,
Constant,
Function,
extract_function_spaces,
functionspace,
assemble_scalar,
dirichletbc,
form,
locate_dofs_geometrical,
)
from dolfinx.fem.petsc import NonlinearProblem, LinearProblem
from dolfinx.io import XDMFFile
import pyvista as pv
import pyvista
try:
import pyvista as pv
import pyvistaqt as pvqt
have_pyvista = True
except ModuleNotFoundError:
print("pyvista and pyvistaqt are required to visualise the solution")
have_pyvista = False
from dolfinx.plot import vtk_mesh
from dolfinx import fem
import ufl
# Save all logging to file
log.set_output_file("log.txt")
rank = MPI.COMM_WORLD.rank
######ORIZW TIS PARAMETROUS POU EINAI STATHERES
r0=5
e4=0.05
Delta=0.65
W0=1
D=1
tau0=1
l=1.5957*D
lamda =1.5957*D
dt =0.016 ## XRONIKO BIMA
t=0
T= 150 # XRONIKI PERIODOS
MPI.COMM_WORLD.Barrier()
start_time = MPI.Wtime()
######################### MESH NETGEN
geo=SplineGeometry()
pnts= [(-400,-400),(400,-400),(400,400),(-400,400)]
p1,p2,p3,p4 = [geo.AppendPoint(*pnt) for pnt in pnts]
curves = [
[["line",p1,p2],"line"],
[["line",p2,p3],"line"],
[["line",p3,p4],"line"],
[["line",p4,p1],"line"],
]
for c,bc in curves :
geo.Append(c,bc=bc)
#####################################
geoModel = ngfx.GeometricModel(geo, MPI.COMM_WORLD)
domain, (ct, ft), region_map = geoModel.model_to_mesh(gdim=2, hmax=15)
########################################################
def mark_cells_radius_DG(domain, R_ref, center=(0.0, 0.0)):
W = fem.functionspace(domain, ("DG", 0))
x = ufl.SpatialCoordinate(domain)
cx, cy = center
r = sqrt((x[0] - cx)**2 + (x[1] - cy)**2)
indicator = ufl.conditional(ufl.le(r, R_ref), 1.0, 0.0)
markers = fem.Function(W)
ip = W.element.interpolation_points
if Version(dolfinx.__version__) < Version("0.10.0"):
ip = ip()
markers.interpolate(fem.Expression(indicator, ip))
# --- IMPORTANT: keep only owned cells (exclude ghosts) ---
tdim = domain.topology.dim
n_owned = domain.topology.index_map(tdim).size_local
arr_owned = markers.x.array[:n_owned]
cells = np.flatnonzero(arr_owned > 0.5).astype(np.int32)
return cells
mesh_id = 0
out_folder = Path("code_3")
out_folder.mkdir(parents=True, exist_ok=True)
# γράψε αρχικό mesh
# geometric initial refinement
n_init_ref = 8
R_ref = 6
tdim = domain.topology.dim
# καλό να υπάρχει connectivity που χρειάζεται το compute_incident_entities
domain.topology.create_connectivity(tdim, tdim-1)
for k in range(n_init_ref):
cells_local = mark_cells_radius_DG(domain, R_ref, center=(0.0, 0.0))
# 100% ασφάλεια: μόνο owned (όχι ghost) ΚΑΙ σωστό dtype
n_owned = domain.topology.index_map(tdim).size_local
cells_local = cells_local[cells_local < n_owned].astype(np.int32, copy=False)
# Debug check (αν θες)
# assert cells_local.size == 0 or (cells_local.min() >= 0 and cells_local.max() < n_owned)
domain, (ct, ft) = geoModel.refineMarkedElements(tdim, cells_local)
import sys
comm = MPI.COMM_WORLD
rank = comm.rank
tdim = domain.topology.dim
n_local = domain.topology.index_map(tdim).size_local
all_locals = comm.gather(n_local, root=0)
comm.Barrier()
if rank == 0:
arr = np.array(all_locals)
print("Local cells per rank:", arr.tolist())
print("min/mean/max:", int(arr.min()), float(arr.mean()), int(arr.max()))
print("max/mean:", float(arr.max()/arr.mean()))
sys.stdout.flush()
###################### dokimazw na plotarw to mesh na dw pws fainetai pyhto
grid = pyvista.UnstructuredGrid(*vtk_mesh(domain))
grid.cell_data["ct"] = ct.values
plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True, scalars="ct", show_scalar_bar=True)
plotter.view_xy()
if not pyvista.OFF_SCREEN:
plotter.show()
######################################
import numpy as np
from dolfinx import geometry
from mpi4py import MPI
class TipTrackerRoot:
def __init__(self, mesh, x_min, x_max, y0=0.0,
x0_guess=None, # αρχική εικασία για tip
tol_x=1e-10, tol_f=1e-12, max_it=200,
scan_n=60, expand_factor=1.6,
out_path="tip_velocity.csv", comm=None):
self.mesh = mesh
self.comm = comm if comm is not None else MPI.COMM_WORLD
self.rank = self.comm.rank
self.x_min = float(x_min)
self.x_max = float(x_max)
self.y0 = float(y0)
self.tol_x = float(tol_x)
self.tol_f = float(tol_f)
self.max_it = int(max_it)
self.scan_n = int(scan_n)
self.expand_factor = float(expand_factor)
self.bb = geometry.bb_tree(self.mesh, self.mesh.topology.dim)
self.x_tip_old = None
self.x_guess = x0_guess
self.f = None
if self.rank == 0:
self.f = open(out_path, "w", buffering=1)
self.f.write("step,t,x_tip,v_tip\n")
def _eval_phi_at(self, phi, x):
"""
Αξιολόγηση φ στο point (x, y0) σε MPI με ασφάλεια:
- κάθε rank δοκιμάζει αν το point πέφτει σε local cell
- όσοι βρίσκουν τιμή συνεισφέρουν, μετά κάνουμε average
(σε partition boundaries μπορεί να βρεθεί από 2 ranks, αλλά η τιμή είναι ίδια για CG).
"""
pt = np.array([[float(x), self.y0, 0.0]], dtype=np.float64)
candidates = geometry.compute_collisions_points(self.bb, pt)
colliding = geometry.compute_colliding_cells(self.mesh, candidates, pt)
local_found = 0
local_val = 0.0
links = colliding.links(0)
if len(links) > 0:
cell = np.array([links[0]], dtype=np.int32)
local_val = float(phi.eval(pt, cell).reshape(-1)[0])
local_found = 1
found = self.comm.allreduce(local_found, op=MPI.SUM)
val_sum = self.comm.allreduce(local_val, op=MPI.SUM)
if found == 0:
return np.nan
return val_sum / found
def _find_bracket(self, phi, x_center):
"""
Βρίσκει [a,b] με sign change γύρω από x_center.
1) προσπαθεί με expanding window
2) αν αποτύχει, κάνει coarse scan σε όλο [x_min, x_max]
"""
# clamp center
x_center = float(np.clip(x_center, self.x_min, self.x_max))
# αρχικό μικρό παράθυρο
width = max(1e-6, 0.01 * (self.x_max - self.x_min))
for _ in range(40):
a = max(self.x_min, x_center - width)
b = min(self.x_max, x_center + width)
fa = self._eval_phi_at(phi, a)
fb = self._eval_phi_at(phi, b)
if np.isfinite(fa) and np.isfinite(fb) and fa * fb < 0:
return a, b, fa, fb
width *= self.expand_factor
# fallback: coarse scan across full interval
xs = np.linspace(self.x_min, self.x_max, self.scan_n)
fs = np.array([self._eval_phi_at(phi, x) for x in xs], dtype=np.float64)
ok = np.isfinite(fs)
xs, fs = xs[ok], fs[ok]
if xs.size < 2:
return None
s = np.sign(fs)
idx = np.where(s[:-1] * s[1:] < 0)[0]
if idx.size == 0:
return None
j = idx[0] # πρώτο crossing (πιο αριστερό)
return float(xs[j]), float(xs[j+1]), float(fs[j]), float(fs[j+1])
def _bisect(self, phi, a, b, fa, fb):
"""
Bisection μέχρι tol_x ή tol_f.
"""
if fa * fb > 0:
return None
for _ in range(self.max_it):
m = 0.5 * (a + b)
fm = self._eval_phi_at(phi, m)
if not np.isfinite(fm):
# αν πέσει σε “κενό” λόγω γεωμετρίας, μικρή μετατόπιση
m = np.nextafter(m, b)
fm = self._eval_phi_at(phi, m)
if np.isfinite(fm) and abs(fm) < self.tol_f:
return float(m)
if abs(b - a) < self.tol_x:
return float(0.5 * (a + b))
# keep bracket
if fa * fm < 0:
b, fb = m, fm
else:
a, fa = m, fm
return float(0.5 * (a + b))
def update(self, phi, step, t, dt):
# initial guess
if self.x_tip_old is None:
x_center = self.x_guess if self.x_guess is not None else 0.5*(self.x_min+self.x_max)
else:
x_center = self.x_tip_old
br = self._find_bracket(phi, x_center)
if br is None:
x_tip = np.nan
else:
a, b, fa, fb = br
x_tip = self._bisect(phi, a, b, fa, fb)
# velocity
if self.x_tip_old is None or not np.isfinite(x_tip):
v_tip = 0.0
else:
v_tip = (x_tip - self.x_tip_old) / dt
self.x_tip_old = x_tip
if self.rank == 0 and self.f is not None:
self.f.write(f"{step},{t},{x_tip},{v_tip}\n")
return x_tip, v_tip
def close(self):
if self.rank == 0 and self.f is not None:
self.f.close()
##################################
def tau(phi):
gradphi=grad(phi)
absgradphi = sqrt(inner( gradphi, gradphi)+1e-12)
abs4gradphi=absgradphi**4
sum_grad_phi4= gradphi[0]**4 + gradphi[1]**4
whole = sum_grad_phi4/abs4gradphi
a=(1-3*e4)*(1+((4*e4)/(1-3*e4))*whole)
return tau0*(a)**2
def W(phi):
gradphi=grad(phi)
absgradphi = sqrt(inner( gradphi, gradphi)+1e-12)
abs4gradphi=absgradphi**4
sum_grad_phi4= gradphi[0]**4 + gradphi[1]**4
whole = sum_grad_phi4/abs4gradphi
a=(1-3*e4)*(1+((4*e4)/(1-3*e4))*whole)
return W0*a
def omega(phi):
gradphi=grad(phi)
absgradphi = sqrt(inner( gradphi, gradphi)+1e-12)
abs4gradphi=absgradphi**4
V1= gradphi[0]*(( gradphi[1]**2)*(( gradphi[0]**2)-( gradphi[1]**2)))
V2= gradphi[1]*(( gradphi[0]**2)*(( gradphi[1]**2)-( gradphi[0]**2)))
V=as_vector([V1,V2])
return ((16*e4*W0)/abs4gradphi)*V
def f(phi,u):
return (phi-l*u*(1-phi**2))*(1-phi**2)
def build_system(domain, dt, D, petsc_options,
old_sol=None, old_sol0=None,
init_phi=None, init_u=None,
petsc_prefix="Crystal_growth"):
"""
Χτίζει: ME, sol, sol0, problem πάνω στο current domain (mesh).
- An old_sol/old_sol0 δοθούν: μεταφέρει λύση nonmatching (για mesh change).
- Αλλιώς, an init_phi/init_u δοθούν: βάζει τα initial initial conditions .
"""
# 1) Spaces + functions
P1 = element("Lagrange", domain.basix_cell(), 1, dtype=default_real_type)
ME = fem.functionspace(domain, mixed_element([P1, P1]))
sol = fem.Function(ME)
sol0 = fem.Function(ME)
# 2) είτε μεταφορά λύσης (mesh changed) είτε initial conditions
if old_sol is not None and old_sol0 is not None:
# Old components
old_phi = old_sol.sub(0)
old_u = old_sol.sub(1)
old_phi0 = old_sol0.sub(0)
old_u0 = old_sol0.sub(1)
# Collapsed spaces + maps
Vphi_new, map_phi = ME.sub(0).collapse()
Vu_new, map_u = ME.sub(1).collapse()
phi_new = fem.Function(Vphi_new)
u_new = fem.Function(Vu_new)
phi0_new = fem.Function(Vphi_new)
u0_new = fem.Function(Vu_new)
# Nonmatching interpolation data
tdim = domain.topology.dim
num_local_cells = domain.topology.index_map(tdim).size_local
new_cells = np.arange(num_local_cells, dtype=np.int32)
idata_phi = fem.create_interpolation_data(Vphi_new, old_phi.function_space, new_cells)
idata_u = fem.create_interpolation_data(Vu_new, old_u.function_space, new_cells)
# Interpolate old -> new
phi_new.interpolate_nonmatching(old_phi, new_cells, idata_phi)
phi0_new.interpolate_nonmatching(old_phi0, new_cells, idata_phi)
u_new.interpolate_nonmatching(old_u, new_cells, idata_u)
u0_new.interpolate_nonmatching(old_u0, new_cells, idata_u)
# Write back into mixed sol/sol0
sol.sub(0).x.array[map_phi] = phi_new.x.array
sol.sub(1).x.array[map_u] = u_new.x.array
sol0.sub(0).x.array[map_phi] = phi0_new.x.array
sol0.sub(1).x.array[map_u] = u0_new.x.array
sol.x.scatter_forward()
sol0.x.scatter_forward()
else:
# initial conditions (στην αρχή)
if init_phi is not None:
sol.sub(0).interpolate(init_phi)
if init_u is not None:
sol.sub(1).interpolate(init_u)
sol.x.scatter_forward()
sol0.x.array[:] = sol.x.array
sol0.x.scatter_forward()
# 3) Weak form (πάντα ξαναχτίζεται στο νέο ME/sol/sol0)
q, v = ufl.TestFunctions(ME)
phi, u = ufl.split(sol)
phi0, u0 = ufl.split(sol0)
F1 = (
ufl.inner(u, v) * ufl.dx
- ufl.inner(u0, v) * ufl.dx
+ dt * D * ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
- 0.5 * ufl.inner(phi, v) * ufl.dx
+ 0.5 * ufl.inner(phi0, v) * ufl.dx
)
F2 = (
ufl.inner(tau(phi0) * phi, q) * ufl.dx
- ufl.inner(tau(phi0) * phi0, q) * ufl.dx
+ dt * ufl.inner(W(phi0)**2 * ufl.grad(phi), ufl.grad(q)) * ufl.dx
- dt * ufl.inner(f(phi0, u0), q) * ufl.dx
+ dt * ufl.inner(W(phi0) * omega(phi0), ufl.grad(q)) * ufl.dx
)
F = F1 + F2
# 4) Problem
problem = NonlinearProblem(F, sol,
petsc_options_prefix=petsc_prefix,
petsc_options=petsc_options)
return ME, sol, sol0, problem
def mark_cells_multiband_buffered(phi: fem.Function,
core=0.15, h0=0.18,
out1=0.40, in1=0.20, h1=0.30,
out2=0.75, in2=0.35, h2=0.50,
out3=0.95, in3=0.50, h3=0.80):
mesh = phi.function_space.mesh
W = fem.functionspace(mesh, ("DG", 0))
ip = W.element.interpolation_points
if Version(dolfinx.__version__) < Version("0.10.0"):
ip = ip()
phi_dg = fem.Function(W)
phi_dg.interpolate(fem.Expression(phi, ip))
phi_dg.x.scatter_forward()
hfun = fem.Function(W)
hfun.interpolate(fem.Expression(ufl.CellDiameter(mesh), ip))
hfun.x.scatter_forward()
tdim = mesh.topology.dim
n_owned = mesh.topology.index_map(tdim).size_local
p = phi_dg.x.array[:n_owned]
h = hfun.x.array[:n_owned]
marked = (
((np.abs(p) < core) & (h > h0)) |
((p > -out1) & (p < in1) & (h > h1)) |
((p > -out2) & (p < in2) & (h > h2)) |
((p > -out3) & (p < in3) & (h > h3))
)
return np.flatnonzero(marked).astype(np.int32)
def mark_cells_multiband_outside(phi: fem.Function,
eps0=0.10, h0=0.12, # very fine around phi=0 on both sides
neg1=0.30, h1=0.22, # first outer band: -neg1 < phi < -eps0
neg2=0.60, h2=0.40, # second outer band: -neg2 < phi < -neg1
neg3=0.90, h3=0.70): # third outer band: -neg3 < phi < -neg2
mesh = phi.function_space.mesh
W = fem.functionspace(mesh, ("DG", 0))
ip = W.element.interpolation_points
if Version(dolfinx.__version__) < Version("0.10.0"):
ip = ip()
# cellwise phi
phi_dg = fem.Function(W)
phi_dg.interpolate(fem.Expression(phi, ip))
phi_dg.x.scatter_forward()
# cell diameter
hfun = fem.Function(W)
hfun.interpolate(fem.Expression(ufl.CellDiameter(mesh), ip))
hfun.x.scatter_forward()
tdim = mesh.topology.dim
n_owned = mesh.topology.index_map(tdim).size_local
p = phi_dg.x.array[:n_owned]
h = hfun.x.array[:n_owned]
marked = (
((p > -eps0) & (p < eps0) & (h > h0)) | # symmetric thin band around phi=0
((p <= -eps0) & (p > -neg1) & (h > h1)) | # outside only
((p <= -neg1) & (p > -neg2) & (h > h2)) | # farther outside
((p <= -neg2) & (p > -neg3) & (h > h3)) # even farther outside
)
cells = np.flatnonzero(marked).astype(np.int32)
return cells
def expand_marked_cells(mesh, cells, n_layers=1):
tdim = mesh.topology.dim
fdim = tdim - 1
mesh.topology.create_connectivity(tdim, fdim)
mesh.topology.create_connectivity(fdim, tdim)
c_to_f = mesh.topology.connectivity(tdim, fdim)
f_to_c = mesh.topology.connectivity(fdim, tdim)
n_owned = mesh.topology.index_map(tdim).size_local
marked = set()
for c in np.asarray(cells, dtype=np.int64).ravel():
if 0 <= c < n_owned:
marked.add(int(c))
for _ in range(n_layers):
new_marked = set(marked)
for c in marked:
for f in c_to_f.links(c):
for cn in f_to_c.links(f):
if 0 <= int(cn) < n_owned:
new_marked.add(int(cn))
marked = new_marked
return np.array(sorted(marked), dtype=np.int32)
#ΔΗΛΩΣΗ ΤΟΥ ΤΙ OPTIONS ΘΑ ΕΧΕΙ Ο SOLVER #######
use_superlu = PETSc.IntType == np.int64
sys = PETSc.Sys()
if sys.hasExternalPackage("mumps") and not use_superlu:
linear_solver = "mumps"
elif sys.hasExternalPackage("superlu_dist"):
linear_solver = "superlu_dist"
else:
linear_solver = "petsc"
petsc_options = {
"snes_type": "newtonls",
"snes_linesearch_type": "basic",
"snes_atol": 1e-6,
"snes_rtol": 1e-6,
"snes_monitor": None,
"ksp_error_if_not_converged": True,
"ksp_type": "gmres",
"ksp_rtol": 1e-8,
"ksp_monitor": None,
"pc_type": "hypre",
"pc_hypre_type": "boomeramg",
"pc_hypre_boomeramg_max_iter": 1,
"pc_hypre_boomeramg_cycle_type": "v",
}
###################################################ARXIKI SINTHIKI GIA TO PHI PARAGONTAS PHI ######
def Initial_condition_PHI(x):
r1 = np.sqrt(x[0]**2 + x[1]**2)
return -np.tanh((r1 - r0) / np.sqrt(2))
##############ARXIKI ADIASTATIS SIGKENTRWSIS%%%%%%%%%%%%%%%%%
def Initial_condition_U(x):
r2 = np.sqrt(x[0]**2 + x[1]**2)
u_init = np.where(r2<=r0,
0,
-Delta * (1 - np.exp(-(r2 -r0)/1))
)
return u_init
######################################################ARXIKI DILWSI TOY PROBLIMATOS
ME, sol, sol0, problem = build_system(
domain, dt, D, petsc_options,
init_phi=Initial_condition_PHI,
init_u=Initial_condition_U
)
out_folder = Path("code_3")
out_folder.mkdir(parents=True, exist_ok=True)
file = XDMFFile(MPI.COMM_WORLD,out_folder / f"output_mesh{mesh_id}.xdmf", "w")
file.write_mesh(domain)
log_path = out_folder / "residuals.txt"
if rank == 0:
f_log = open(log_path, "w", buffering=1) # line-buffered
f_log.write("step,t,snes_reason,snes_its,snes_norm,ksp_its,ksp_rnorm\n")
else:
f_log = None
###############################
step = 0
# --- Initial refinement passes (before time loop) ---
##############################arxiko refinement###############
gdim = domain.geometry.dim
# ∇phi (vector)
def rebuild_postprocess_objects(domain, tip_tracker=None):
gdim = domain.geometry.dim
Vgrad = fem.functionspace(domain, ("CG", 1, (gdim,)))
grad_u = fem.Function(Vgrad, name="grad_u")
grad_phi = fem.Function(Vgrad, name="grad_phi")
if tip_tracker is not None:
tip_tracker.mesh = domain
tip_tracker.bb = geometry.bb_tree(domain, domain.topology.dim)
return Vgrad, grad_u, grad_phi
def keep_owned_cells(mesh, cells):
tdim = mesh.topology.dim
imap = mesh.topology.index_map(tdim)
n_owned = imap.size_local
cells = np.asarray(cells, dtype=np.int64).ravel()
cells = cells[(cells >= 0) & (cells < n_owned)]
cells = np.unique(cells)
return cells.astype(np.int32, copy=False)
tip_tracker = TipTrackerRoot(
mesh=domain,
x_min=-400.0,
x_max=0.0,
y0=0.0,
x0_guess=-r0, # προαιρετικό: αρχική εικασία
tol_x=1e-10, tol_f=1e-12, # πολύ αυστηρά
out_path=str(out_folder / "tip_velocity.csv"),
comm=MPI.COMM_WORLD
)
Vgrad, grad_u, grad_phi = rebuild_postprocess_objects(domain, tip_tracker)
while t < T:
t += dt
step += 1
problem.solve()
sol.x.scatter_forward()
# diagnostics
converged_reason = problem.solver.getConvergedReason()
num_iterations = problem.solver.getIterationNumber()
snes = problem.solver
ksp = snes.getKSP()
snes_reason = snes.getConvergedReason()
snes_its = snes.getIterationNumber()
# SNES residual norm
snes_norm = snes.getFunctionNorm()
# KSP stats (τελευταίο linear solve)
ksp_its = ksp.getIterationNumber()
ksp_rnorm = ksp.getResidualNorm()
if step % 99 == 0:
old_sol_fine = sol.copy()
old_sol0_fine = sol0.copy()
geoModel = ngfx.GeometricModel(geo, MPI.COMM_WORLD)
domain, (ct, ft), region_map = geoModel.model_to_mesh(gdim=2, hmax=10.0)
max_rebuild_ref = 10
for k in range(max_rebuild_ref):
ME_tmp, sol_tmp, sol0_tmp, _ = build_system(
domain, dt, D, petsc_options,
old_sol=old_sol_fine, old_sol0=old_sol0_fine
)
cells = mark_cells_multiband_buffered(
sol_tmp.sub(0),
core=0.15, h0=0.18,
out1=0.40, in1=0.20, h1=0.30,
out2=0.75, in2=0.35, h2=0.50,
out3=0.95, in3=0.50, h3=0.80)
cells = expand_marked_cells(domain, cells, n_layers=1)
cells = keep_owned_cells(domain, cells)
local_marked = cells.size
global_marked = comm.allreduce(local_marked, op=MPI.SUM)
if rank == 0:
print(f"[rebuild pass {k+1}] global_marked={global_marked}", flush=True)
if global_marked < 10:
break
tdim = domain.topology.dim
n_owned = domain.topology.index_map(tdim).size_local
n_all = domain.topology.index_map(tdim).size_local + domain.topology.index_map(tdim).num_ghosts
if rank == 0:
print(f"About to refine: cells.size={cells.size}, dtype={cells.dtype}", flush=True)
if cells.size > 0:
print(f"cells.min()={cells.min()}, cells.max()={cells.max()}", flush=True)
print(f"n_owned={n_owned}, n_all(local+ghost)={n_all}", flush=True)
assert cells.dtype == np.int32
assert np.all(cells >= 0)
assert np.all(cells < n_owned), f"Found non-owned cells: max={cells.max()} n_owned={n_owned}"
assert len(np.unique(cells)) == len(cells), "Duplicate cell ids in refine list"
domain, (ct, ft) = geoModel.refineMarkedElements(domain.topology.dim, cells)
mesh_id += 1
file.close()
file = XDMFFile(MPI.COMM_WORLD, out_folder / f"output_mesh{mesh_id}.xdmf", "w")
file.write_mesh(domain)
ME, sol, sol0, problem = build_system(
domain, dt, D, petsc_options,
old_sol=old_sol_fine, old_sol0=old_sol0_fine
)
Vgrad, grad_u, grad_phi = rebuild_postprocess_objects(domain, tip_tracker)
if step % 50 == 0:
cells = mark_cells_multiband_buffered(
sol.sub(0),
core=0.15, h0=0.2,
out1=0.40, in1=0.20, h1=0.30,
out2=0.75, in2=0.35, h2=0.60,
out3=0.95, in3=0.50, h3=1.00)
cells = expand_marked_cells(domain, cells, n_layers=1)
cells = keep_owned_cells(domain, cells)
local_marked = cells.size
global_marked = comm.allreduce(local_marked, op=MPI.SUM)
if rank == 0:
print(f"[step {step}] global_marked={global_marked}", flush=True)
if global_marked >= 20:
old_sol = sol.copy()
old_sol0 = sol0.copy()
print(f"Refining {global_marked} cells at step {step}", flush=True)
domain, (ct, ft) = geoModel.refineMarkedElements(domain.topology.dim, cells)
domain.topology.create_connectivity(domain.topology.dim, domain.topology.dim - 1)
domain.topology.create_connectivity(domain.topology.dim - 1, domain.topology.dim)
mesh_id += 1
file.close()
file = XDMFFile(MPI.COMM_WORLD, out_folder / f"output_mesh{mesh_id}.xdmf", "w")
file.write_mesh(domain)
ME, sol, sol0, problem = build_system(
domain, dt, D, petsc_options,
old_sol=old_sol, old_sol0=old_sol0
)
Vgrad, grad_u, grad_phi = rebuild_postprocess_objects(domain, tip_tracker)
phi = sol.sub(0)
u = sol.sub(1)
# Γράψιμο ΜΕΣΑ στο loop
if rank == 0:
f_log.write(f"{step},{t},{snes_reason},{snes_its},{snes_norm},{ksp_its},{ksp_rnorm}\n")
# 1) διάνυσμα gradient
phi = sol.sub(0)
u = sol.sub(1)
tip_tracker.update(phi, step=step, t=t, dt=dt)
ip_grad = Vgrad.element.interpolation_points
if Version(dolfinx.__version__) < Version("0.10.0"):
ip_grad = ip_grad()
grad_phi.interpolate(fem.Expression(ufl.grad(phi), ip_grad))
grad_u.interpolate(fem.Expression(ufl.grad(u), ip_grad))
sol0.x.array[:] = sol.x.array
sol0.x.scatter_forward()
grad_phi.x.scatter_forward()
grad_u.x.scatter_forward()
# Write output
if step % 1==0:
phi_out = sol.sub(0)
u_out = sol.sub(1)
file.write_function(phi_out, t)
file.write_function(u_out, t)
file.write_function(grad_u,t)
file.write_function(grad_phi,t)
file.flush()
tip_tracker.close()
file.close()
if rank == 0 and f_log is not None:
f_log.close()
if rank == 0:
print(f"Simulation completed: {step} steps, final time t={t}")
MPI.COMM_WORLD.Barrier()
end_time = MPI.Wtime()
print(f"Elapsed time: {end_time-start_time} s.",flush=True) #remove later
if rank == 0:
print(f"Simulation completed: {step} steps, final time t={t}")
MPI.COMM_WORLD.Barrier()
end_time = MPI.Wtime()
print(f"Elapsed time: {end_time-start_time} s.",flush=True) #remove later
dolfinx.common.list_timings(MPI.COMM_WORLD)
Hi, I am running a phase-field simulation where I do adaptive mesh refinement and coarsening during the time loop by rebuilding a new mesh and then interpolating the old solution onto the new mesh .
More specifically, I periodically create a new mesh, transfer the solution with interpolation / nonmatching interpolation, and continue the solve on the updated mesh. However, after around 2300+ time steps , the code consistently crashes with a PETSc segmentation fault during the mesh adaptation stage.
The error I keep getting is:
[0]PETSC ERROR: ------------------------------------------------------------------------
[0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range
[0]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger
[0]PETSC ERROR: or see https://petsc.org/release/faq/#valgrind and https://petsc.org/release/faq/
[0]PETSC ERROR: configure using --with-debugging=yes, recompile, link, and run
[0]PETSC ERROR: to get more information on the crash.
[0]PETSC ERROR: Run with -malloc_debug to check if memory corruption is causing the crash.
Abort(59) on node 0 (rank 0 in comm 0): application called MPI_Abort(MPI_COMM_WORLD, 59) - process 0
The nonlinear solve itself converges normally before that, and the crash appears during the remeshing / refinement part, not during the Newton solve. This works for many iterations, but after roughly 2300 time steps the code consistently crashes with a PETSc SEGV .
Has anyone seen this kind of issue when repeatedly doing adaptive remeshing with interpolation onto a new mesh?