Thank you so much Andrash and Dokken,
i will try to find time to test it out this week, this addition to the package will be very usefull for simulating the effect of a polarized gate near non-easy to define materials.
i hope to be able to make the MWE work by myself
Thanks again!
Hi, i am a quite new user of linux, and i have some trouble understanding dockering the nightly build of Dolfinx in my python env, if you can point me on what i did wrong or should change, i’ll be very thankfull as it’s been hours i’m on it :
some context :
I use virtualbox to have a VirtualMachine with ubuntu 22.04, in which i used anaconda to create a python environment named “fenicsx”.
what i did so far :
in ubuntu terminal using sudo i installed docker, and give my user the permissions on docker group. (following How to fix docker: Got permission denied issue - Stack Overflow)
Back to “fenicsx” environment, i pip installed docker and the 0.9.0dev0 version of external operator pkg. After that i uninstalled fenics-dolfinx 0.8.0. and i still have basix ufl ffcx
Now i’m trying to launch the test the file of Dokken from the link you provided Andrash, using the following command at the beginning of the test file,
import docker
client = docker.from_env()
client.images.pull('dolfinx/dolfinx:nightly')# i comment it after the first use
client.images.list()
client.containers.run('dolfinx/dolfinx:nightly')
and this at the end of the script :
test_external_operator_codim_1(1)
i get a “no module named doflinx” error message
do i miss some code for docker to actually make the nightly build of dolfinx work with my python environment?
Thanks for reading
I wouldn’t use the python “docker” package to run docker.
I would run in from the terminal (you should read a general docker tutorial): https://youtu.be/pTFZFxd4hOI?feature=shared&t=1845
However, if you insist of running your code through python, then you could done something like:
import argparse
import docker
import os
from pathlib import Path
parser = argparse.ArgumentParser(
description="Run a minimal working example",
formatter_class=argparse.ArgumentDefaultsHelpFormatter,
)
parser.add_argument(
"--filename", type=Path, required=True, help="Name of the file to run"
)
if __name__ == "__main__":
args = parser.parse_args()
client = docker.from_env()
client.images.pull(
"ghcr.io/fenics/dolfinx/dolfinx:nightly"
) # i comment it after the first use
client.images.list()
lines = client.containers.run(
"ghcr.io/fenics/dolfinx/dolfinx:nightly",
f"python3 /tmp/{args.filename}",
volumes={os.getcwd(): {"bind": "/tmp/", "mode": "rw"}},
stream=True,
)
for line in lines:
print(line)
which can be run with
python3 run_docker_script.py --filename name_of_script_to_run.py
where you have your script in name_of_script_to_run.py.
1 Like
Thanks for the reply, i now use docker in command window as i couldn’t get it totally right with python.
i will soon try to test the new external_ope_pkg with the test file and work on a MWE, but i noticed something that is bothering me :
The docker image of dolfinx nighty build is shipped with python 3.12.3, unfortunately one of the packages i must use, pybinding, cannot be installed/build on python 3.11 and onwards.
Is there any way for me to have dolfinx 0.9.0 dev0 with 3.10 python environment? (i have pybinding working in 3.10.13 with dolfinx 0.8.0).
Also, will there soon be an update to dolfinx 0.9.0 that can be installed easily via conda install, and that work with the newer external operator package?
Thanks for all the replies so far
Have a nice day
You would have to install dolfinx manually to get this working.
There will be a 0.9.0 release this fall, but I do not know exactly when this will happen.
1 Like
Hi again,
i made some progress in the implementation of the external op pkg into a custom Newton solver. As a reminder i want to have a Neumann BC where ‘g’ is dependant on the potential at the border. Where the way it depends on it, is calculated with some other python script, hence the use of the external operator pkg.
So in the Docker, with nighlty Dolfinx 0.9.0.dev0, after installing the external operator pkg, i was able to test the codim1 test function of Dokken provided in the link by Latyshev. It seems to work fine, no error occured when asking : test_external_operator_codim_1(2)
After that i went on changing the Newton solver i previously made.
The Newton solver does converge and in a nicer way than before, but not at the same values (for correction norm and potential) as when using g from UFL. see image below :
with the external operator pkg :
with g as UFL :
Quite the difference for the potential at the bottom border, external operator pkg give a flat 0.618V at the bottom. i’m struggling to understand where is my mistake?
I didn’t know what was codim 0 and 1 , so i asked chatgpt about what it means in fenicsx, so normally sonce i’m working on a border, i should look at and use the codim1 test function for my usecase, and that’s what i did, but at this point i’m not sure of anything anymore. If you can find some time this week or next week to look at the MWE, i’ll be gratefull.
Thanks again
this is the MWE Newton solver with external operator pkg :
"""
Created on Tue Jun 18 13:02:24 2024
In 2d planar square system, solve laplacian(uh) = 0 in Omega, uh = 1 in d_omega_D (top border), g = -uh² in d_omega_N (Neumann bottom border)
g use external operator package
"""
import dolfinx
from dolfinx import default_scalar_type, log, fem
from dolfinx.fem import (Constant, Function, functionspace,petsc,
assemble_scalar, dirichletbc, form, locate_dofs_geometrical,locate_dofs_topological)
from dolfinx.mesh import create_unit_square, locate_entities_boundary, meshtags, locate_entities
from dolfinx.plot import vtk_mesh
from mpi4py import MPI
import ufl
from ufl import * #SpatialCoordinate, TestFunction, TrialFunction, dot, ds, dx, grad, inner, Measure
from petsc4py.PETSc import ScalarType
from petsc4py import PETSc
import numpy as np
import pyvista
import matplotlib.pyplot as plt
from dolfinx_external_operator import (
FEMExternalOperator,
evaluate_external_operators,
evaluate_operands,
replace_external_operators)
import basix
def g_impl(uuh):
output = np.zeros(np.size(uuh))
output = -1*pow(uuh,2) # NUMPY
return output.reshape(-1)# The output must be returned flattened to one dimension
def dgdu_impl(uuh):
output1 = np.zeros(np.size(uuh))
output1 = -2*uuh # NUMPY
return output1.reshape(-1) # The output must be returned flattened to one dimension
def g_external(derivatives): # the switch between g and dg/du
if derivatives == (0,): # no derivation, the function itself
return g_impl
elif derivatives == (1,): # the derivative with respect to the operand `uh`
return dgdu_impl
else:
return NotImplementedError
##################### Code #######################################
quadrature_degree=2
mesh = create_unit_square(MPI.COMM_WORLD, 5, 5)
mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
boundaries = [(1, lambda x: x[1]<=0.001)] # I want the bottom border of the unit square
facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
for (marker, locator) in boundaries:
facets = locate_entities(mesh, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full_like(facets, marker))
ext_facets = np.hstack(facet_indices).astype(np.int32)
# print('facet_indices low border'),print(facet_indices),print(' ')
V = functionspace(mesh, ("CG", 1))
uh = Function(V)
v = TestFunction(V)
a = inner(grad(uh), grad(v)) * dx
x = SpatialCoordinate(mesh)
## From Dokken test file :
submesh, sub_to_parent, _, _ = dolfinx.mesh.create_submesh(mesh, mesh.topology.dim - 1, ext_facets)
num_entities = mesh.topology.index_map(mesh.topology.dim - 1).size_local
parent_to_sub = np.full(num_entities, -1, dtype=np.int32)
parent_to_sub[sub_to_parent] = np.arange(len(sub_to_parent), dtype=np.int32)
entity_maps = {submesh: parent_to_sub}
Qe = basix.ufl.quadrature_element(submesh.basix_cell(), degree=quadrature_degree, value_shape=())
Q = dolfinx.fem.functionspace(submesh, Qe)
####### g with external operator package :
g = FEMExternalOperator(uh, function_space=Q, external_function=g_external) # g is now Symbolic not numpy involved
####### ds fabrication for Neumann BC, bottom of square ####
ft = dolfinx.mesh.meshtags(mesh, mesh.topology.dim - 1, ext_facets, np.full_like(ext_facets, 1))
ds = ufl.Measure(
"ds", domain=mesh, subdomain_data=ft, subdomain_id=1, metadata={"quadrature_degree": quadrature_degree}
)
##################### Dirichlet BC on top +1V #####
dofs_up = locate_dofs_geometrical(V, lambda x: x[1]>=0.999)
u_up = Function(V)
u_up.interpolate(lambda x: x[0]*0+ 1)
bc_up = dirichletbc(u_up, dofs_up)
bcs = [bc_up]
##########################################################
du = dolfinx.fem.Function(V)
maxiter = 7
i=0
correct_list=[]
L1 = inner(g,v) * ds # g symbolic used, this should be Neumann BC
F=a-L1
while i < maxiter:
######### Assemble Jacobian J and residual F ############
F_replaced, F_external_operators = replace_external_operators(F)
evaluated_operands = evaluate_operands(F_external_operators, entity_maps=entity_maps)
_ = evaluate_external_operators(F_external_operators, evaluated_operands)
F_compiled = dolfinx.fem.form(F_replaced, entity_maps=entity_maps)
J = derivative(F, uh) # define derivative
J_expanded = ufl.algorithms.expand_derivatives(J)
J_replaced, J_external_operators = replace_external_operators(J_expanded)
evaluated_operands = evaluate_operands(J_external_operators, entity_maps=entity_maps)
_ = evaluate_external_operators(J_external_operators, evaluated_operands)
J_compiled = dolfinx.fem.form(J_replaced, entity_maps=entity_maps)
#########################################################
A = dolfinx.fem.petsc.create_matrix(J_compiled)
L = dolfinx.fem.petsc.create_vector(F_compiled)
solver = PETSc.KSP().create(mesh.comm)
solver.setOperators(A)
dolfinx.fem.petsc.assemble_matrix(A, J_compiled, bcs=bcs)
A.assemble()
dolfinx.fem.petsc.assemble_vector(L, F_compiled)
L.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
L.scale(-1)
dolfinx.fem.petsc.apply_lifting(L, [J_compiled], [bcs], x0=[uh.vector], scale=1)
dolfinx.fem.petsc.set_bc(L, bcs, uh.vector, 1.0)
L.ghostUpdate(addv=PETSc.InsertMode.INSERT_VALUES, mode=PETSc.ScatterMode.FORWARD)
solver.solve(L, du.vector)
du.x.scatter_forward()
uh.x.array[:] += du.x.array
correction_norm = du.vector.norm(0) # Compute norm of update
correct_list.append(correction_norm)
print(f"Iteration {i}: Correction norm {correction_norm}")
i += 1
if correction_norm < 1e-8:
break
#### PLOTTING BELOW
# plt.figure(1)
# plt.semilogy(range(i),correct_list,'-+')
# plt.title('correction norm (iteration)')
# plt.savefig('correction_norm.pdf')
# from dolfinx import plot
# plt.figure(2) # get uh at y=y_bottom for all x
# cells, types, x = plot.vtk_mesh(V)
# UH = uh.x.array.real
# y=x[:,1]
# x2=x[:,0]
# ind=np.where(y<0.01)[0] #
# UHX = UH[ind]
# X = x[ind,0]
# sort=np.argsort(X) # index in increasing order for X
# X=X[sort]### applied to X
# UHX=UHX[sort] ### applied to UHX
# plt.plot(X,UHX,'-')
# plt.title('Potentiel(x) y=y$_{bottom}$')
# plt.savefig('pot_Neum_BC.pdf')
this is the MWE Nweton solver with g UFL :
"""
Created on Tue Jun 18 13:02:24 2024
In 2d planar square system, solve laplacian(uh) = 0 in Omega, uh = 1 in d_omega_D (top border), g = -uh² in d_omega_N (Neumann bottom border)
g use external operator package
"""
import dolfinx
from dolfinx.fem import (Constant, Function, functionspace,petsc,
assemble_scalar, dirichletbc, form, locate_dofs_geometrical,locate_dofs_topological)
from dolfinx.mesh import create_unit_square, locate_entities_boundary, meshtags, locate_entities
from dolfinx.plot import vtk_mesh
from mpi4py import MPI
import ufl
from ufl import * #SpatialCoordinate, TestFunction, TrialFunction, dot, ds, dx, grad, inner, Measure
from petsc4py.PETSc import ScalarType
from petsc4py import PETSc
import numpy as np
import pyvista
import matplotlib.pyplot as plt
##################### Code #######################################
mesh = create_unit_square(MPI.COMM_WORLD, 5, 5)
V = functionspace(mesh, ("CG", 1))
uh = Function(V)
v = TestFunction(V)
a = inner(grad(uh), grad(v)) * dx
x = SpatialCoordinate(mesh)
##################### ds fabrication for Neumann BC, bottom of square ####
boundaries = [(1, lambda x: x[1]<=0.001)]
facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
for (marker, locator) in boundaries:
facets = locate_entities(mesh, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
ds = Measure("ds", domain=mesh, subdomain_data=facet_tag)
##################### Dirichlet BC on top +1V #####
dofs_up = locate_dofs_geometrical(V, lambda x: x[1]>=0.999)
u_up = Function(V)
u_up.interpolate(lambda x: x[0]*0+ 1)
bc_up = dirichletbc(u_up, dofs_up)
bcs = [bc_up]
####### g(uh ufl) :
g = -1*pow(uh,2) ### Full UFL formulation for comparaison purpose
#########################################################
du = dolfinx.fem.Function(V)
maxiter = 10
i=0
correct_list=[]
L1 = inner(g,v) * ds # g symbolic used, this should be Neumann BC
F=a-L1
while i < maxiter:
residual = dolfinx.fem.form(F)
J = ufl.derivative(F, uh)
jacobian = dolfinx.fem.form(J)
############################################
A = dolfinx.fem.petsc.create_matrix(jacobian)
L = dolfinx.fem.petsc.create_vector(residual)
solver = PETSc.KSP().create(mesh.comm)
solver.setOperators(A)
dolfinx.fem.petsc.assemble_matrix(A, jacobian, bcs=bcs)
A.assemble()
dolfinx.fem.petsc.assemble_vector(L, residual)
L.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
L.scale(-1)
# Compute b - J(u_D-u_(i-1))
dolfinx.fem.petsc.apply_lifting(L, [jacobian], [bcs], x0=[uh.vector], scale=1)
# Set du|_bc = u_{i-1}-u_D
dolfinx.fem.petsc.set_bc(L, bcs, uh.vector, 1.0)
L.ghostUpdate(addv=PETSc.InsertMode.INSERT_VALUES, mode=PETSc.ScatterMode.FORWARD)
# Solve linear problem
solver.solve(L, du.vector)
du.x.scatter_forward()
# Update u_{i+1} = u_i + delta u_i
uh_old=uh.copy()
uh.x.array[:] += du.x.array
# Compute norm of update
correction_norm = du.vector.norm(0)
correct_list.append(correction_norm)
print(f"Iteration {i}: Correction norm {correction_norm}")
i += 1
if correction_norm < 1e-8:
break
### PLOTTING BELOW
# plt.figure(1)
# plt.semilogy(range(i),correct_list,'-+')
# plt.title('correction norm (iteration)')
# plt.figure(2)
# pyvista_cells, cell_types, geometry = vtk_mesh(V)
# grid = pyvista.UnstructuredGrid(pyvista_cells, cell_types, geometry)
# grid.point_data["uh"] = uh.x.array
# grid.set_active_scalars("uh")
# plotter = pyvista.Plotter()
# plotter.add_text("uh", position="upper_edge", font_size=14, color="black")
# plotter.add_mesh(grid, show_edges=True)
# plotter.view_xy()
# if not pyvista.OFF_SCREEN: ## this makes dolfinx 0.9.0dev0 : PETSC ERROR:Segementation violation, Caught signal number 11 SEGV
# plotter.show()
# else:
# figure = plotter.screenshot("multiple_dirichlet.png")
# print('out of loop10')
# from dolfinx import plot
# plt.figure(2) # get uh at y=y_bottom for all x
# cells, types, x = plot.vtk_mesh(V)
# UH = uh.x.array.real
# y=x[:,1]
# x2=x[:,0]
# ind=np.where(y<0.01)[0] #
# UHX = UH[ind]
# X = x[ind,0]
# sort=np.argsort(X) # index in increasing order for X
# X=X[sort]### applied to X
# UHX=UHX[sort] ### applied to UHX
# plt.plot(X,UHX,'-')
# plt.title('Potentiel(x) y=y$_{bottom}$')
Hi again,
i was wondering why the potential solution with external op pkg looked better than with UFL. I mean in UFL MWE, the corners of the unit square have a big impact on the potential.
So i decided to redo the MWE on Pdetool in MATLAB, and I got the exact same solution as with the external operator package.
I think i can assume now that the MWE with ext op pkg does give the right solution, and that my UFL MWE script is actually incorrect.
I think the Neumann BC on bottom is not correctly applied to the bottom of the square for the ufl mwe If you can express your thought on this point, that would be nice.
In conclusion, i think we achieved the goal of this thread, thanks again Dokken and Latyshev.
Consider the following, where I have combined the two methods, and have ensured that you use the same the same quadrature rule for both gs and enforce them on the same boundary, yielding:
"""
Created on Tue Jun 18 13:02:24 2024
In 2d planar square system, solve laplacian(uh) = 0 in Omega, uh = 1 in d_omega_D (top border), g = -uh² in d_omega_N (Neumann bottom border)
g use external operator package
"""
import dolfinx
import basix.ufl
from dolfinx.fem import (Function, functionspace,
dirichletbc, locate_dofs_geometrical)
from dolfinx.mesh import create_unit_square, locate_entities
import dolfinx.fem.petsc
from mpi4py import MPI
import ufl
from ufl import *
from petsc4py import PETSc
import numpy as np
import argparse
parser = argparse.ArgumentParser()
impl = parser.add_mutually_exclusive_group(required=True)
impl.add_argument("--ufl", action="store_true")
impl.add_argument("--external", action="store_true")
args = parser.parse_args()
if args.ufl:
use_ufl = True
elif args.external:
use_ufl = False
else:
raise ValueError("Invalid implementation")
##################### Code #######################################
mesh = create_unit_square(MPI.COMM_WORLD, 5, 5)
V = functionspace(mesh, ("CG", 1))
uh = Function(V)
v = TestFunction(V)
a = inner(grad(uh), grad(v)) * dx
x = SpatialCoordinate(mesh)
##################### ds fabrication for Neumann BC, bottom of square ####
boundaries = [(1, lambda x: x[1]<=0.001)]
facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
for (marker, locator) in boundaries:
facets = locate_entities(mesh, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tags = dolfinx.mesh.meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
##################### Dirichlet BC on top +1V #####
dofs_up = locate_dofs_geometrical(V, lambda x: x[1]>=0.999)
u_up = Function(V)
u_up.interpolate(lambda x: x[0]*0+ 1)
bc_up = dirichletbc(u_up, dofs_up)
bcs = [bc_up]
####### g(uh ufl) :
quadrature_degree = 2
if use_ufl:
g = -1*pow(uh,2) ### Full UFL formulation for comparaison purpose
else:
print("applying external operator")
from dolfinx_external_operator import (
FEMExternalOperator,
evaluate_external_operators,
evaluate_operands,
replace_external_operators)
def g_impl(uuh):
output = np.zeros(np.size(uuh))
output = -1*pow(uuh,2) # NUMPY
return output.reshape(-1)# The output must be returned flattened to one dimension
def dgdu_impl(uuh):
output1 = np.zeros(np.size(uuh))
output1 = -2*uuh # NUMPY
return output1.reshape(-1) # The output must be returned flattened to one dimension
def g_external(derivatives): # the switch between g and dg/du
if derivatives == (0,): # no derivation, the function itself
return g_impl
elif derivatives == (1,): # the derivative with respect to the operand `uh`
return dgdu_impl
else:
return NotImplementedError
submesh, sub_to_parent, _, _ = dolfinx.mesh.create_submesh(mesh, facet_tags.dim, facet_tags.find(1))
num_entities = mesh.topology.index_map(mesh.topology.dim - 1).size_local
parent_to_sub = np.full(num_entities, -1, dtype=np.int32)
parent_to_sub[sub_to_parent] = np.arange(len(sub_to_parent), dtype=np.int32)
entity_maps = {submesh: parent_to_sub}
Qe = basix.ufl.quadrature_element(submesh.basix_cell(), degree=quadrature_degree, value_shape=())
Q = dolfinx.fem.functionspace(submesh, Qe)
####### g with external operator package :
g = FEMExternalOperator(uh, function_space=Q, external_function=g_external) # g is now Symbolic not numpy involved
#########################################################
du = dolfinx.fem.Function(V)
maxiter = 10
i=0
correct_list=[]
ds = ufl.Measure(
"ds", domain=mesh, subdomain_data=facet_tags, subdomain_id=1, metadata={"quadrature_degree": quadrature_degree}
)
L1 = inner(g,v) * ds # g symbolic used, this should be Neumann BC
F=a-L1
J = derivative(F, uh)
if use_ufl:
jacobian = dolfinx.fem.form(J)
residual = dolfinx.fem.form(F)
else:
F_replaced, F_external_operators = replace_external_operators(F)
J_expanded = ufl.algorithms.expand_derivatives(J)
J_replaced, J_external_operators = replace_external_operators(J_expanded)
jacobian = dolfinx.fem.form(J_replaced, entity_maps=entity_maps)
residual = dolfinx.fem.form(F_replaced, entity_maps=entity_maps)
A = dolfinx.fem.petsc.create_matrix(jacobian)
L = dolfinx.fem.petsc.create_vector(residual)
solver = PETSc.KSP().create(mesh.comm)
solver.setOperators(A)
print(f"Implementation UFL:{use_ufl}")
while i < maxiter:
############################################
if not use_ufl:
evaluated_operands = evaluate_operands(J_external_operators, entity_maps=entity_maps)
_ = evaluate_external_operators(J_external_operators, evaluated_operands)
evaluated_operands = evaluate_operands(F_external_operators, entity_maps=entity_maps)
_ = evaluate_external_operators(F_external_operators, evaluated_operands)
A.zeroEntries()
dolfinx.fem.petsc.assemble_matrix(A, jacobian, bcs=bcs)
A.assemble()
with L.localForm() as loc:
loc.set(0)
dolfinx.fem.petsc.assemble_vector(L, residual)
L.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
L.scale(-1)
# Compute b - J(u_D-u_(i-1))
dolfinx.fem.petsc.apply_lifting(L, [jacobian], [bcs], x0=[uh.vector], scale=1)
# Set du|_bc = u_{i-1}-u_D
dolfinx.fem.petsc.set_bc(L, bcs, uh.vector, 1.0)
L.ghostUpdate(addv=PETSc.InsertMode.INSERT_VALUES, mode=PETSc.ScatterMode.FORWARD)
# Solve linear problem
solver.solve(L, du.vector)
du.x.scatter_forward()
# Update u_{i+1} = u_i + delta u_i
uh.x.array[:] += du.x.array
# Compute norm of update
correction_norm = du.vector.norm(0)
correct_list.append(correction_norm)
print(f"Iteration {i}: Correction norm {correction_norm}")
i += 1
if correction_norm < 1e-8:
break
if use_ufl:
implementation = "ufl"
else:
implementation = "external"
with dolfinx.io.XDMFFile(mesh.comm, f"u_{implementation}.xdmf", "w") as xdmf:
xdmf.write_mesh(mesh)
xdmf.write_function(uh)
1 Like
Thanks for the nice script Dokken,
so it was indeed my mistake in ufl. Hopefully now i have everything woking well. Have a nice day