Update comments and explanations in PhaseField examples

examples/PhaseField/plot_2000.py

@@ -4,8 +4,11 @@
 Crack surface density functional
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
-In the following example, we consider the boundary value problem of the phase-field model, which represents the crack surface density functional, thus providing a continuous approximation of the discontinuous crack  (:ref:`theory_phase_field`). Due to the symmetry of the problem, only the left half of the bar is considered. Therefore, a boundary condition is applied at the left end of this half bar, as illustrated in the following diagrams.
-For a one-dimensional simulation, the boundary condition $\\phi=1$ is set at the point located on the left end.
+In this example, we examine the boundary value problem related to the representation of a crack surface. This approach models the crack surface density functional, offering a continuous approximation of an otherwise discontinuous crack surface (:ref:`theory_phase_field`).
+
+Due to the symmetry of the problem, we consider only the left half of the bar. A boundary condition is therefore applied at the left end of this half-bar, as shown in the diagrams below.
+
+For a one-dimensional simulation, the boundary condition $\phi=1$ is imposed at the left end of this segment.
 
 .. code-block::
 
@@ -75,56 +78,116 @@
 ###############################################################################
 # Parameters definition
 # ---------------------
-# A length scale parameter of  $l = 4$  is defined. The phase field is saved in a VTU file.
-# All the results are stored in a folder named results_folder_name = name.
-Data = Input(l=4.0,
-             save_solution_xdmf=False,
-             save_solution_vtu=True,
-             results_folder_name="2000_General")
+# First, we define an input class, which contains all the parameters needed for the setup 
+# and results of the simulation.
+#
+# The first term, $l$, specifies the length scale parameter for the problem, with $l = 4.0$.
+#
+# The next two options, `save_solution_xdmf` and `save_solution_vtu`, determine the file format 
+# used to save the phase-field results (.xdmf or .vtu), which can then be visualized using 
+# tools like ParaView or pvista. Both parameters are boolean values (`True` or `False`). 
+# In this case, we set `save_solution_vtu=True` to save the phase-field results in .vtu format.
+#
+# Lastly, `results_folder_name` specifies the name of the folder where all results 
+# and log information will be saved. If the folder does not exist, `phasefieldx` will create it. 
+# However, if the folder already exists, any previous data in it will be removed, and a new 
+# empty folder will be created in its place.
+#
+Data = Input(
+    l=4.0,
+    save_solution_xdmf=False,
+    save_solution_vtu=True,
+    results_folder_name="2000_General"
+)
+
 
 
 ###############################################################################
-# Mesh definition
+# Mesh Definition
 # ---------------
+# We define the mesh parameters and set up a two-dimensional simulation. This setup 
+# supports various dimensions (‘1D’, ‘2D’, or ‘3D’) by creating a mesh that consists 
+# of either a line, rectangle, or box, depending on the selected dimension, with 
+# corresponding line, quadrilateral, or hexahedral elements.
+#
+# - `divx`, `divy`, and `divz` specify the number of divisions along the x, y, and z 
+#   axes, respectively. Here, `divx=100`, `divy=1`, and `divz=1` are set to divide 
+#   the x-axis primarily, as needed for a 2D or 1D mesh.
+#
+# - `lx`, `ly`, and `lz` define the physical dimensions of the domain in the x, y, and z 
+#   directions. In this example, we set `lx=5.0`, `ly=1.0`, and `lz=1.0`.
+#
+# Specify the simulation dimension with the `dimension` variable (`"1d"`, `"2d"`, or `"3d"`).
+# Here, we choose `"2d"`.
 divx, divy, divz = 100, 1, 1
 lx, ly, lz = 5.0, 1.0, 1.0
-
-
-###############################################################################
-# A two-dimensional simulation is considered
-# (it is also possible to select ‘1D’ and ‘3D’ simulations)
 dimension = "2d"
 
+# Mesh creation based on specified dimension
 if dimension == "1d":
-    msh = dolfinx.mesh.create_interval(mpi4py.MPI.COMM_WORLD,
-                                       divx,
-                                       np.array([0.0, lx]))
+    # Creates a 1D mesh, which consists of a line divided into `divx` line elements, 
+    # extending from 0 to lx along the x-axis.
+    msh = dolfinx.mesh.create_interval(
+        mpi4py.MPI.COMM_WORLD,
+        divx,
+        np.array([0.0, lx])
+    )
 
 elif dimension == "2d":
-    msh = dolfinx.mesh.create_rectangle(mpi4py.MPI.COMM_WORLD,
-                                        [np.array([0.0, 0.0]),
-                                         np.array([lx, ly])],
-                                        [divx, divy],
-                                        cell_type=dolfinx.mesh.CellType.quadrilateral)
+    # Creates a 2D mesh, which consists of a rectangle covering [0, 0] to [lx, ly] with `divx` 
+    # divisions in x and `divy` divisions in y, using quadrilateral elements.
+    msh = dolfinx.mesh.create_rectangle(
+        mpi4py.MPI.COMM_WORLD,
+        [np.array([0.0, 0.0]), np.array([lx, ly])],
+        [divx, divy],
+        cell_type=dolfinx.mesh.CellType.quadrilateral
+    )
 
 elif dimension == "3d":
-    msh = dolfinx.mesh.create_box(mpi4py.MPI.COMM_WORLD,
-                                  [np.array([0.0, 0.0, 0.0]),
-                                   np.array([lx, ly, lz])],
-                                  [divx, divy, divz],
-                                  cell_type=dolfinx.mesh.CellType.hexahedron)
-
-
-###############################################################################
-# The left part is denoted and shown to impose the boundary conditions
+    # Creates a 3D mesh, which consists of a box extending from [0, 0, 0] to [lx, ly, lz], 
+    # divided into `divx`, `divy`, and `divz` parts along the x, y, and z axes, respectively, 
+    # with hexahedral elements.
+    msh = dolfinx.mesh.create_box(
+        mpi4py.MPI.COMM_WORLD,
+        [np.array([0.0, 0.0, 0.0]), np.array([lx, ly, lz])],
+        [divx, divy, divz],
+        cell_type=dolfinx.mesh.CellType.hexahedron
+    )
+
+
+# Left Boundary Identification
+# ----------------------------
+# This function identifies points on the left side of the domain where the boundary 
+# condition will be applied. Specifically, it returns `True` for points where `x=0`,
+# and `False` otherwise. This allows us to selectively apply boundary conditions 
+# only to this part of the mesh.
 def left(x):
     return np.equal(x[0], 0)
 
 
+# `fdim` represents the dimension of the boundary facets on the mesh, which is one 
+# less than the mesh's overall dimensionality (`msh.topology.dim`). For example, 
+# if the mesh is 2D, `fdim` will be 1, representing 1D boundary edges.
 fdim = msh.topology.dim - 1
 
+# Using the `left` function, we locate the facets on the left side of the mesh 
+# where `x=0`. The `locate_entities_boundary` function returns an array of facet 
+# indices that represent the identified boundary entities.
 left_facet_marker = dolfinx.mesh.locate_entities_boundary(msh, fdim, left)
+
+
+# `get_ds_bound_from_marker` is a function that generates a measure for integrating 
+# boundary conditions specifically on the facets identified by `left_facet_marker`. 
+# This measure is assigned to `ds_left` and will be used for applying boundary 
+# conditions on the left side.
 ds_left = get_ds_bound_from_marker(left_facet_marker, msh, fdim)
+
+
+# `ds_list` is an array that stores boundary condition measures and associated 
+# names for each boundary to facilitate result-saving processes. Each entry in 
+# `ds_list` is an array in the form `[ds_, "name"]`, where `ds_` is the boundary 
+# condition measure, and `"name"` is a label for saving purposes. Here, `ds_left` 
+# is labeled as `"left"` for clarity when saving results.
 ds_list = np.array([
                    [ds_left, "left"],
                    ])
@@ -139,20 +202,53 @@ def left(x):
 
 
 ###############################################################################
-# Boundary Condition
-# ------------------
-# The boundary conditions of $\phi=1$ is set on the left part of the mesh.
+# Boundary Condition Setup for Scalar Field $\phi$
+# ------------------------------------------------
+# We define and apply a Dirichlet boundary condition for the scalar field $\phi$
+# on the left side of the mesh, setting $phi = 1$ on this boundary. This setup is 
+# for a simple, static linear problem, meaning the boundary conditions and loading 
+# are constant and do not change throughout the simulation.
+#
+# - `bc_phi` is a function that creates a Dirichlet boundary condition on a specified 
+#   facet of the mesh for the scalar field $\phi$.
+# - `bcs_list_phi` is a list that stores all the boundary conditions for $\phi$, 
+#   facilitating easy management and extension of conditions if needed.
+# - `update_boundary_conditions` and `update_loading` are set to `None` as they are 
+#   unused in this static case with constant boundary conditions and loading.
+
 bc_left = bc_phi(left_facet_marker, V_phi, fdim, value=1.0)
 bcs_list_phi = [bc_left]
 update_boundary_conditions = None
 update_loading = None
 
 
 ###############################################################################
-# Call the solver
-# ---------------
-# We set a final time of t=1 and a time step of \Delta t=1.
-# Note that this is a linear boundary value problem.
+# Solver Call for a Static Linear Problem
+# ---------------------------------------
+# We define the parameters for a simple, static linear boundary value problem 
+# with a final time `t = 1.0` and a time step `Δt = 1.0`. Although this setup 
+# includes time parameters, they are primarily used for structural consistency 
+# with a generic solver function and do not affect the result, as the problem 
+# is linear and time-independent.
+#
+# Parameters:
+# - `final_time`: The end time for the simulation, set to 1.0.
+# - `dt`: The time step for the simulation, set to 1.0. In a static context, this
+#   only provides uniformity with dynamic cases but does not change the results.
+# - `path`: Optional path for saving results; set to `None` here to use the default.
+# - `quadrature_degree`: Defines the accuracy of numerical integration; set to 2 
+#   for this problem.
+#
+# Function Call:
+# The `solve` function is called with:
+# - `Data`: Simulation data and parameters.
+# - `msh`: Mesh of the domain.
+# - `V_phi`: Function space for `phi`.
+# - `bcs_list_phi`: List of boundary conditions.
+# - `update_boundary_conditions`, `update_loading`: Set to `None` as they are 
+#   unused in this static problem.
+# - `ds_list`: Boundary measures for integration on specified boundaries.
+# - `dt` and `final_time` to define the static solution time window.
 
 final_time = 1.0
 dt = 1.0
@@ -215,35 +311,70 @@ def left(x):
 
 
 ###############################################################################
-# Plot: Energy values
-# -------------------
-# The energy values are compared with the analytic ones.
-l_array = np.linspace(0.1, lx, 100)
-a_div_l = lx / l_array
-tanh_a_div_l = np.tanh(a_div_l)
+# Plot: Energy Values Comparison
+# -------------------------------
+# In this section, we compare the energy values obtained from the simulation 
+# with their corresponding analytic expressions.
+#
+# The energy components are calculated for the scalar field `phi` and its gradient. 
+# We compare the following energy terms:
+# - `W_phi`: The energy associated with the scalar field `phi`.
+#            $W_{\phi} = \frac{1}{2l} \int_{-a}^{a}   \left[ e^{-|x|/l} + \frac{1}{e^{\frac{2a}{l}}+1} 2 \sinh \left( \frac{|x|}{l} \right) \right]^2 dx$
+# - `W_gradphi`: The energy associated with the gradient of the scalar field.
+#                $W_{\nabla \phi}(\phi) = \frac{l}{2} \int_{-a}^{a}   \left[ \frac{-\text{sign}(x)}{l} e^{-|x|/l} + \frac{1}{e^{\frac{2a}{l}} +1} \frac{\text{sign}(x)}{l} 2 \cosh\left(\frac{|x|}{l}\right) \right]^2 dx$
+# - `W`: The total energy, which includes contributions from both the field and its gradient.
+#        $W = \tanh \left( \frac{a}{l} \right)$
+#
+# Theoretical expressions are computed for these energies as functions of the 
+# length scale parameter `l`. Additionally, energy values from the simulation 
+# are overlaid on the theoretical curve for comparison.
+
+l_array = np.linspace(0.1, lx, 100)  # Create an array of length scale values
+a_div_l = lx / l_array  # Compute the ratio of `lx` to the length scale
+tanh_a_div_l = np.tanh(a_div_l)  # Compute the hyperbolic tangent of the ratio
+
+# Compute the theoretical energy for the scalar field `phi`
+# The expression represents the energy of `phi` based on the length scale parameter `l`.
 W_phi = 0.5 * tanh_a_div_l + 0.5 * a_div_l * (1.0 - tanh_a_div_l**2)
+
+# Compute the theoretical energy for the gradient of `phi`
+# This energy is based on the gradient of `phi` and the length scale parameter `l`.
 W_gradphi = 0.5 * tanh_a_div_l - 0.5 * a_div_l * (1.0 - tanh_a_div_l**2)
-W = np.tanh(a_div_l)
 
+# Compute the total theoretical energy
+# This energy is related to the hyperbolic tangent of the ratio `a/l`.
+W = np.tanh(a_div_l)
 
 ###############################################################################
-# The next graph shows the theoretical energy versus the length scale parameter,
-# as well as those corresponding to the length scale parameter used in the
-# simulation. It is noted that the solution coincides with the theoretical one.
-fig, energy = plt.subplots()
-
-plt.rc('text', usetex=True)  # Enable LaTeX rendering
-energy.plot(l_array, W_phi, 'r-', label=r'$W_{\phi}$')
-energy.plot(l_array, W_gradphi, 'b-', label=r'$W_{\nabla \phi}$')
-energy.plot(l_array, W, 'k-', label=r'$W$')
-
+# Plotting Energy Comparison
+# ---------------------------
+# The following plot compares the theoretical energy values (`W_phi`, `W_gradphi`, 
+# and `W`) with the energy values obtained from the simulation for various 
+# length scale parameters `l`. The simulation results are plotted as points 
+# on top of the theoretical curves to assess the agreement between theory and 
+# the computed values.
+
+fig, energy = plt.subplots()  # Create a figure for plotting energy
+
+# Enable LaTeX rendering for plot labels
+plt.rc('text', usetex=True)
+
+# Plot the theoretical energy values
+energy.plot(l_array, W_phi, 'r-', label=r'$W_{\phi}$')  # Energy for `phi`
+energy.plot(l_array, W_gradphi, 'b-', label=r'$W_{\nabla \phi}$')  # Energy for gradient of `phi`
+energy.plot(l_array, W, 'k-', label=r'$W$')  # Total energy
+
+# Plot the energy values obtained from the simulation
+# The values are scaled by 2 and correspond to the total energy components.
 energy.plot(Data.l, 2 * S.energy_files['total.energy']["gamma_phi"][0], 'k*', label=S.label)
 energy.plot(Data.l, 2 * S.energy_files['total.energy']["gamma_gradphi"][0], 'k*')
 energy.plot(Data.l, 2 * S.energy_files['total.energy']["gamma"][0], 'k*')
 
-energy.set_xlabel('l')
-energy.set_ylabel('energy')
-energy.grid(color='k', linestyle='-', linewidth=0.3)
-energy.legend()
+# Set plot labels and grid
+energy.set_xlabel('Length Scale Parameter $l$')  # X-axis label
+energy.set_ylabel('Energy')  # Y-axis label
+energy.grid(color='k', linestyle='-', linewidth=0.3)  # Grid settings
+energy.legend()  # Display the legend
+
 
 plt.show()