diff --git a/.github/workflows/publish-docs.yml b/.github/workflows/publish-docs.yml index 4cd4892..371f078 100644 --- a/.github/workflows/publish-docs.yml +++ b/.github/workflows/publish-docs.yml @@ -35,12 +35,12 @@ jobs: uses: actions/download-artifact@v4 with: name: webpage - path: "./doc/public" + path: "_build/html" - name: Upload pages artifact uses: actions/upload-pages-artifact@v3 with: - path: "./doc/public" + path: "_build/html" - name: Deploy to GitHub Pages id: deployment diff --git a/doc/demo/demo_plasticity_mohr_coulomb.py b/doc/demo/demo_plasticity_mohr_coulomb.py index 4fe378d..beb7140 100644 --- a/doc/demo/demo_plasticity_mohr_coulomb.py +++ b/doc/demo/demo_plasticity_mohr_coulomb.py @@ -715,7 +715,7 @@ def F_ext(v): if len(points_on_process) > 0: results[i + 1, :] = (-u.eval(points_on_process, cells)[0], load) -print(f"Slope stability factor: {-q.value[-1]*H/c}") +print(f"Slope stability factor: {-q.value[-1] * H / c}") # %% [markdown] # ## Verification diff --git a/doc/demo/demo_plasticity_von_mises.py b/doc/demo/demo_plasticity_von_mises.py index 4e6098f..5ae5c2e 100644 --- a/doc/demo/demo_plasticity_von_mises.py +++ b/doc/demo/demo_plasticity_von_mises.py @@ -150,19 +150,20 @@ # ### Preamble # %% -from mpi4py import MPI +from mpi4py import MPI # noqa: F401 from petsc4py import PETSc import matplotlib.pyplot as plt import numba import numpy as np from demo_plasticity_von_mises_pure_ufl import plasticity_von_mises_pure_ufl -from solvers import LinearProblem from utilities import build_cylinder_quarter, find_cell_by_point import basix import ufl from dolfinx import fem +from dolfinx.fem.petsc import NonlinearProblem +from dolfinx.nls.petsc import NewtonSolver from dolfinx_external_operator import ( FEMExternalOperator, evaluate_external_operators, @@ -203,7 +204,7 @@ # %% k_u = 2 -V = fem.functionspace(mesh, ("Lagrange", k_u, (2,))) +V = fem.functionspace(mesh, ("Lagrange", k_u, (mesh.geometry.dim,))) # Boundary conditions bottom_facets = facet_tags.find(facet_tags_labels["Lx"]) left_facets = facet_tags.find(facet_tags_labels["Ly"]) @@ -245,10 +246,10 @@ def epsilon(v): loading = fem.Constant(mesh, PETSc.ScalarType(0.0)) v = ufl.TestFunction(V) -F = ufl.inner(sigma, epsilon(v)) * dx - ufl.inner(-loading * n, v) * ds(facet_tags_labels["inner"]) +F = ufl.inner(sigma, epsilon(v)) * dx - ufl.inner(loading * -n, v) * ds(facet_tags_labels["inner"]) # Internal state -P_element = basix.ufl.quadrature_element(mesh.topology.cell_name(), degree=k_stress, value_shape=()) +P_element = basix.ufl.quadrature_element(mesh.topology.cell_name(), degree=k_stress) P = fem.functionspace(mesh, P_element) p = fem.Function(P, name="cumulative_plastic_strain") @@ -275,10 +276,10 @@ def epsilon(v): # $\frac{\mathrm{d} \boldsymbol{\sigma}}{\mathrm{d} \boldsymbol{\varepsilon}}$ # must be implemented by the user. In this tutorial, we implement the derivative using the Numba package. # -# First of all, we implement the return-mapping procedure locally in the function -# `_kernel`. It computes the values of the stress tensor, the tangent moduli and -# the increment of cumulative plastic strain at a single Gausse node. For more -# details, visit the [original +# First of all, we implement the return-mapping procedure locally in the +# function `_kernel`. It computes the values of the stress tensor, the tangent +# moduli and the increment of cumulative plastic strain at a single Gausse +# node. For more details, visit the [original # implementation](https://comet-fenics.readthedocs.io/en/latest/demo/2D_plasticity/vonMises_plasticity.py.html) # of this problem for the legacy FEniCS 2019. # @@ -412,91 +413,70 @@ def sigma_external(derivatives): # first loading step, so to initialize the former, we evaluate external operators # for a close-to-zero displacement field. # -# At the same time, we estimate the compilation overhead caused by the first call -# of JIT-ed Numba functions. - -# %% -# We need to initialize `Du` with small values in order to avoid the division by -# zero error -eps = np.finfo(PETSc.ScalarType).eps -Du.x.array[:] = eps - -# %% [markdown] # ### Solving the problem # -# Solving the problem is carried out in a manually implemented Newton solver. - # %% u = fem.Function(V, name="displacement") -du = fem.Function(V, name="Newton_correction") -external_operator_problem = LinearProblem(J_replaced, -F_replaced, Du, bcs=bcs) -# %% -# Defining a cell containing (Ri, 0) point, where we calculate a value of u -# In order to run this program in parallel we need capture the process, to which -# this point is attached -x_point = np.array([[R_i, 0, 0]]) -cells, points_on_process = find_cell_by_point(mesh, x_point) -# %% tags=["scroll-output"] -q_lim = 2.0 / np.sqrt(3.0) * np.log(R_e / R_i) * sigma_0 -num_increments = 20 -max_iterations, relative_tolerance = 200, 1e-8 -load_steps = (np.linspace(0, 1.1, num_increments, endpoint=True) ** 0.5)[1:] -loadings = q_lim * load_steps -results = np.zeros((num_increments, 2)) +class PlasticityProblem(NonlinearProblem): + def form(self, x: PETSc.Vec) -> None: + """This function is called before the residual or Jacobian is + computed. This is usually used to update ghost values, but here + we also use it to evaluate the external operators. -evaluated_operands = evaluate_operands(F_external_operators) -evaluate_external_operators(J_external_operators, evaluated_operands) + Args: + x: The vector containing the latest solution + """ + # The following line is from the standard NonlinearProblem class + x.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD) -for i, loading_v in enumerate(loadings): - loading.value = loading_v - external_operator_problem.assemble_vector() + evaluated_operands = evaluate_operands(F_external_operators) + ((_, sigma_new, dp_new),) = evaluate_external_operators(J_external_operators, evaluated_operands) - residual_0 = residual = external_operator_problem.b.norm() - Du.x.array[:] = 0.0 + # This avoids having to evaluate the external operators of F. + sigma.ref_coefficient.x.array[:] = sigma_new + dp.x.array[:] = dp_new - if MPI.COMM_WORLD.rank == 0: - print(f"\nresidual , {residual} \n increment: {i+1!s}, load = {loading.value}") - for iteration in range(0, max_iterations): - if residual / residual_0 < relative_tolerance: - break - external_operator_problem.assemble_matrix() - external_operator_problem.solve(du) - du.x.scatter_forward() +problem = PlasticityProblem(F_replaced, Du, bcs=bcs, J=J_replaced) - Du.x.petsc_vec.axpy(1.0, du.x.petsc_vec) - Du.x.scatter_forward() +x_point = np.array([[R_i, 0, 0]]) +cells, points_on_process = find_cell_by_point(mesh, x_point) - evaluated_operands = evaluate_operands(F_external_operators) +q_lim = 2.0 / np.sqrt(3.0) * np.log(R_e / R_i) * sigma_0 +num_increments = 20 +load_steps = np.linspace(0, 1.1, num_increments, endpoint=True) ** 0.5 +loadings = q_lim * load_steps +results = np.zeros((num_increments, 2)) - # Implementation of an external operator may return several outputs and - # not only its evaluation. For example, `C_tang_impl` returns a tuple of - # Numpy-arrays with values of `C_tang`, `sigma` and `dp`. - ((_, sigma_new, dp_new),) = evaluate_external_operators(J_external_operators, evaluated_operands) +solver = NewtonSolver(mesh.comm, problem) +solver.max_it = 200 +solver.rtol = 1e-8 +ksp = solver.krylov_solver +opts = PETSc.Options() # type: ignore +option_prefix = ksp.getOptionsPrefix() +opts[f"{option_prefix}ksp_type"] = "preonly" +opts[f"{option_prefix}pc_type"] = "lu" +opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps" +ksp.setFromOptions() - # In order to update the values of the external operator we may directly - # access them and avoid the call of - # `evaluate_external_operators(F_external_operators, evaluated_operands).` - sigma.ref_coefficient.x.array[:] = sigma_new - dp.x.array[:] = dp_new +eps = np.finfo(PETSc.ScalarType).eps - external_operator_problem.assemble_vector() - residual = external_operator_problem.b.norm() +for i, loading_v in enumerate(loadings): + loading.value = loading_v + Du.x.array[:] = eps - if MPI.COMM_WORLD.rank == 0: - print(f" it# {iteration} residual: {residual}") + solver.solve(Du) u.x.petsc_vec.axpy(1.0, Du.x.petsc_vec) u.x.scatter_forward() - # Taking into account the history of loading p.x.petsc_vec.axpy(1.0, dp.x.petsc_vec) sigma_n.x.array[:] = sigma.ref_coefficient.x.array if len(points_on_process) > 0: - results[i + 1, :] = (u.eval(points_on_process, cells)[0], loading.value / q_lim) + results[i, :] = (u.eval(points_on_process, cells)[0], loading.value / q_lim) # %% [markdown] # ### Post-processing @@ -524,6 +504,7 @@ def sigma_external(derivatives): plt.ylabel(r"Applied pressure $q/q_{\text{lim}}$ [-]") plt.legend() plt.grid() + plt.savefig("output.png") plt.show() # %% [markdown] diff --git a/doc/demo/demo_plasticity_von_mises_pure_ufl.py b/doc/demo/demo_plasticity_von_mises_pure_ufl.py index e5e81be..f49bf6f 100644 --- a/doc/demo/demo_plasticity_von_mises_pure_ufl.py +++ b/doc/demo/demo_plasticity_von_mises_pure_ufl.py @@ -149,7 +149,7 @@ def von_mises_expressions(Du, old_sig, old_p): nRes = nRes0 if MPI.COMM_WORLD.rank == 0 and verbose: - print(f"\nIncrement#{i+1!s}: load = {t * q_lim:.3f}, Residual0 = {nRes0:.2e}") + print(f"\nIncrement#{i + 1!s}: load = {t * q_lim:.3f}, Residual0 = {nRes0:.2e}") niter = 0 while nRes / nRes0 > tol and niter < Nitermax: diff --git a/doc/demo/utilities.py b/doc/demo/utilities.py index 5313d6c..4ff404e 100644 --- a/doc/demo/utilities.py +++ b/doc/demo/utilities.py @@ -60,7 +60,7 @@ def build_cylinder_quarter(lc=0.3, R_e=1.3, R_i=1.0): # NOTE: Do not forget to check the leaks produced by the following line of code # [WARNING] yaksa: 2 leaked handle pool objects - mesh, cell_tags, facet_tags = gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, rank=model_rank, gdim=2) + mesh, cell_tags, facet_tags, _, _, _ = gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, rank=model_rank, gdim=2) mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim) mesh.name = "quarter_cylinder"