PhaseFieldX Reaction Forces: Available Components and Computation Methods

PhaseFieldX computes reaction forces as vectors R = (Rx, Ry, Rz) at constrained Dirichlet boundaries by extracting the negative residual from the assembled system and summing components according to the problem dimension.

The castillonmiguel/phasefieldx repository provides functionality for calculating reaction forces that develop at constrained boundaries during phase-field simulations. These forces represent the structural response at Dirichlet boundaries and are essential for verifying equilibrium and analyzing boundary interactions. Understanding how to access and compute these values enables accurate post-processing of mechanical behavior in fracture and multi-physics simulations.

Available Reaction Force Components

PhaseFieldX returns reaction forces as Cartesian vectors for each monitored boundary. During configuration, the solver records forces for boundaries specified in the bcs_list_u_names argument.

The available components include:

• Vector components: Each boundary produces a reaction vector R = (Rx, Ry, Rz), containing only the components relevant to the simulation dimension (1D, 2D, or 3D).
• File storage: Forces are written to plain-text files with the .reaction extension, such as bottom.reaction, top.reaction, or bottom_left.reaction.
• Data structure: These files contain tabulated data with the header #step\tRx\tRy\tRz and are loaded into pandas DataFrame objects within the Result class according to src/phasefieldx/PostProcessing/ReferenceResult.py.

Access the data through the Result object:

# Access the y-component of reaction force at the bottom boundary

bottom_Ry = result.reaction_files['bottom.reaction']["Ry"]

How Reaction Forces Are Computed in PhaseFieldX

The core computation resides in calculate_reaction_forces within src/phasefieldx/Reactions/reactions_forces.py. The algorithm derives reaction forces from the finite element residual vector, converting internal stresses into equivalent nodal forces at constrained boundaries.

Step 1 – Assembling the Residual Vector

The process begins by constructing the residual vector F for the current solution field u:

residual_vector = dolfinx.fem.petsc.create_vector(V)
dolfinx.fem.petsc.assemble_vector(residual_vector, F_form)

Step 2 – Applying Dirichlet Boundary Conditions

To incorporate constrained degrees of freedom, the solver lifts Dirichlet contributions into the residual:

dolfinx.fem.petsc.apply_lifting(residual_vector,
                                [J_form], [bcs],
                                x0=[u.x.petsc_vec], alpha=1.0)
dolfinx.fem.petsc.set_bc(residual_vector, bcs,
                          u.x.petsc_vec, alpha=1.0)

Step 3 – Synchronization and Sign Inversion

For parallel consistency, ghost values synchronize across processes. The sign inversion enforces the physical interpretation that reaction = –residual:

residual_vector.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
residual_vector.scale(-1.0)

Step 4 – Component Summation by Dimension

The final step sums residual components belonging to constrained degrees of freedom. The PETSc vector stores data in an interleaved format [rx₀, ry₀, rz₀, rx₁, ry₁, ...]. The function slices this array based on dimension:

• 1D: Rx = sum(residual[0::1])
• 2D: Rx = sum(residual[0::2]), Ry = sum(residual[1::2])
• 3D: Rx = sum(residual[0::3]), Ry = sum(residual[1::3]), Rz = sum(residual[2::3])

This produces a NumPy array [Rx, Ry, Rz] (unused dimensions remain zero) returned to the caller. The solver module src/phasefieldx/Element/Phase_Field_Fracture/solver/solver_ener_variational.py writes these values to the .reaction file at line 431.

Accessing Reaction Force Data in Python

After running a simulation, access pre-computed reaction forces through the result object:

# Load reaction forces from the top boundary

top_data = S.reaction_files['top.reaction']
Rx_history = top_data["Rx"]
Ry_history = top_data["Ry"]
print(f"Final top reaction: Rx={Rx_history.iloc[-1]:.4e}, Ry={Ry_history.iloc[-1]:.4e}")

For custom solvers, invoke the low-level routine directly:

from phasefieldx.Reactions.reactions_forces import calculate_reaction_forces

R = calculate_reaction_forces(
    J_u_form, F_u_form,
    [bc_list_u[i]],  # Single DirichletBC in list

    u, V_u,
    msh.topology.dim  # Spatial dimension (2 or 3)

)

# R contains [Rx, Ry, Rz] as NumPy array

Summary

• PhaseFieldX stores reaction forces as vectors R = (Rx, Ry, Rz) for each constrained boundary defined in bcs_list_u_names.
• The computation occurs in calculate_reaction_forces within src/phasefieldx/Reactions/reactions_forces.py, using the negative residual method.
• Results persist in .reaction files (e.g., bottom.reaction) with tab-separated columns and load into pandas DataFrames via result.reaction_files.
• The algorithm assembles the residual, applies Dirichlet lifting, synchronizes parallel data, inverts the sign, and sums components using stride-based indexing for 1D/2D/3D problems.

Frequently Asked Questions

What file format does PhaseFieldX use to store reaction forces?

PhaseFieldX writes reaction forces to plain-text files with the .reaction extension. These files use a tab-separated format with the header #step\tRx\tRy\tRz, allowing direct import into pandas DataFrames or other analysis tools.

How do I access reaction forces for a specific boundary after a simulation?

Access reaction forces through the Result object's reaction_files dictionary using the boundary name with the .reaction extension, such as result.reaction_files['bottom.reaction']. This returns a pandas DataFrame with columns Rx, Ry, and Rz containing the force history.

Why does the reaction force calculation flip the sign of the residual vector?

The calculation inverts the sign (reaction = –residual) because the finite element residual represents internal forces. At constrained boundaries, the physical reaction force must balance these internal forces, requiring the negative sign to obtain the correct directional sense.

Can I compute reaction forces for custom solvers in PhaseFieldX?

Yes, import calculate_reaction_forces from src/phasefieldx/Reactions/reactions_forces.py and call it within your solver loop. Pass the Jacobian form J_u_form, residual form F_u_form, boundary conditions list, solution vector u, function space V_u, and spatial dimension to receive a NumPy array of reaction components.

{
  "mcpServers": {
    "instagit": {
      "command": "npx",
      "args": ["-y", "instagit@latest"]
    }
  }
}