Convergence issue with KKSMultiACBulk: 2-phase system fails with independent Order Parameters

[js-jixu]

Hi MOOSE Developers and Users,

I am working on reproducing the KKS phase-field model for Zirconium hydrides described in Simon et al. (J. Nucl. Mater., 2021). I am currently simulating a 2D single-variant case (Matrix + 1 Hydride Variant).

I have encountered a numerical stability issue when using the general multi-phase formulation compared to a reduced binary formulation.

The Problem

The model uses the interpolation function proposed by Moelans (2011):

$$ h_{\eta} = \frac{\eta_i^2}{\sum_j \eta_j^2} $$

I created two input files to test this:

1. Input 1 (General Formulation):

   • Variables: eta1 and etam are independent non-linear variables.
   • Kernels: Used KKSMultiACBulkF and KKSMultiACBulkC.
   • Constraint: No explicit Lagrange multiplier enforcing $\sum \eta = 1$, because the paper by Moelans states that it is precisely designed to address the situation where $\sum \eta \neq 1$.
   • Result: Divergence. The solver fails with Linear solve did not converge. I suspect the solver explores a state where both $\eta_1 \approx 0$ and $\eta_m \approx 0$, causing a singularity ($0/0$) in the Jacobian derivative of $h$.

2. Input 2 (Reduced Formulation):

   • Variables: Only eta1 is solved.
   • Material: etam is explicitly defined as 1.0 - eta1:

     # h(eta1)
     [h1_mat]
       type = DerivativeParsedMaterial
       property_name = h1
       coupled_variables = 'eta1'
       expression = 'eta1^2 / ( (1.0 - eta1)^2 + eta1^2)'
     []

   • Result: Converges perfectly.

My Goal

While Input 2 works for this simple case, I intend to extend this to multiple hydride variants (e.g., Matrix + Variant 1 + Variant 2...), where the relation $\eta_m = 1 - \eta_1$ will no longer hold. Therefore, I need to make the General Formulation (Input 1) work robustly.

Questions

1. Handling the Singularity: Is there a standard "best practice" in MOOSE to handle the potential division-by-zero in the KKS interpolation function when using KKSMulti kernels?

   • Should I manually add a small regularization term to the denominator (e.g., + 1.0e-9) in my material derivative expression?
   • Or does MOOSE handle this internally if configured correctly?

2. Solver Strategy: For this type of unconstrained multi-order-parameter system, is PJFNK preferred over NEWTON?

Input Snippet (From the failing Input 1)

Here is how I defined the material causing the suspected issue:

[Variables]
  [eta1]
  []
  [etam]
  []
[]

[Materials]
  # The interpolation function h(eta1)
  # I suspect the denominator causes singularity when solver probes eta ~ 0
  [h1_mat]
    type = DerivativeParsedMaterial
    property_name = h1
    coupled_variables = 'eta1 etam'
    expression = 'eta1^2 / (eta1^2 + etam^2)' 
  []
  
  # ... (similar definition for hm) ...
[]

[Executioner]
  type = Transient
  # I have tried both NEWTON and PJFNK, both struggle without regularization.
  solve_type = 'NEWTON' 
  # ...
[]

---

These are my input files:
input1.txt
input2.txt

Thank you very much for your time and discussion. I really appreciate any insights or suggestions.

[GiudGiud]

Hello

Should I manually add a small regularization term to the denominator (e.g., + 1.0e-9) in my material derivative expression?

Yes you should make sure denominators do not go to 0. Not necessarily with a summed term, you can also do something like:
if (eta > 1e-9, real expression involving the fraction, (else) 0 for the whole term)

Solver Strategy: For this type of unconstrained multi-order-parameter system, is PJFNK preferred over NEWTON?

PJFNK is better suited if the Jacobian is not exact. While DerivativeParsedMaterial are capable of handling derivatives, it could come from somewhere else.
In this page:
https://mooseframework.inl.gov/moose/application_usage/failed_solves.html
we give advice on figuring out if the Jacobian is exact

[js-jixu]

Hi Guillaume,

Thank you so much for the clear explanation! I will immediately try implementing the regularization using the if statement logic you suggested, and I will also switch to PJFNK to handle the inexact Jacobian better.

I noticed that the paper I am trying to reproduce (Simon et al., 2021) was actually authored by the team at INL, and Larry Aagesen is one of the co-authors. It is truly an inspiring piece of work on hydride morphology, which is why I am so keen on reproducing it faithfully.

Since this specific implementation of the Moelans interpolation function is central to that paper, if it is not too much trouble, would it be possible to bring this discussion to @laagesen's attention? His specific insight into how the team handled the unconstrained variables in that work would be incredibly valuable to me.

I will implement these changes immediately and update you with the results as soon as possible.

Thanks again for your help!

[js-jixu]

Hi Guillaume,

I am writing to update you that the convergence issue with the general multi-phase formulation (Input 1) has been successfully resolved!

You were spot on regarding the denominator. As you suggested, the solver was indeed probing states where the order parameters were close to zero, causing a singularity. I implemented a regularization term by adding a small constant (eps = 1.0e-16) to the denominator of the interpolation function.

Combined with the PJFNK , the simulation now runs well.

Thank you very much for your time and helpful suggestions!

Here is the corrected material definition that works:

  [h1_mat]
    type = DerivativeParsedMaterial
    property_name = h1
    coupled_variables = 'eta1 etam' 
    constant_names = 'eps'
    constant_expressions = '1.0e-16' 
    expression = 'eta1^2 / (eta1^2 + etam^2 + eps)'
  []
  [hm_mat]
    type = DerivativeParsedMaterial
    property_name = hm
    coupled_variables = 'eta1 etam'
    constant_names = 'eps'
    constant_expressions = '1.0e-16'
    expression = 'etam^2 / (eta1^2 + etam^2 + eps)'
  []

[Pineapple-blowing-snow]

@js-jixu Hi，Dr.Xu，I recently came across your work on GitHub regarding phase-field modeling of hydride evolution / KKS multiphase simulations], and I found it highly relevant and insightful for my current research. The structure of your model and the implementation details are particularly inspiring.could you kindly share your preferred contact information (e.g., email or other communication channel), or let me know the best way to reach you?Thank you very much for sharing your valuable work with the community. I truly appreciate your time and consideration.
Best regards

[js-jixu]

Thank you for your kind words and for your interest in my work. I am glad to hear that my research on phase-field modeling has been helpful to your project.

I would be more than happy to discuss the model and its implementation details with you. You can reach me via email at [email protected] or [email protected]