Implementing a set of equations

[gj19866]

Hi all,

I would just like to check that my understanding of how MOOSE works is correct before I get really stuck into debugging my code.

I am trying to solve a rather nasty set of equations, for variables $\psi_R$, $\psi_I$, $\Phi$, and $\underline{j}$:

1. $$0=-u \sqrt{1+\gamma^2\left(\psi_R^2+\psi_I^2\right)}\left(\frac{\partial}{\partial t} \psi_R-\varphi \psi_I\right)+\nabla^2 \psi_R-\gamma^2 \psi_I \nabla^2 \varphi+\left(1-\psi_R^2-\psi_I^2\right) \psi_R$$

2. $$0=-u \sqrt{1+\gamma^2\left(\psi_R^2+\psi_I^2\right)}\left(\frac{\partial}{\partial t} \psi_I+\varphi \psi_R\right)+\nabla^2 \psi_I+\gamma^2 \psi_R \nabla^2 \varphi+\left(1-\psi_R^2-\psi_I^2\right) \psi_I$$

3. $$0=\nabla^2 \varphi-\psi_R \nabla^2 \psi_I+\psi_I \nabla^2 \psi_R$$

4. $$0=\underline{j}+\nabla \varphi-\psi_R \nabla \psi_I+\psi_I \nabla \psi_R$$

To do this I have written one kernel for each equation to solve them.

As I have 4 variables, I have to pick a variable for each equation to solve for, but it is arbitrary which I choose, as long as each variable has at least one active kernel. (is this correct?)
The other variables in the equation need to be coupled in as adCoupledValue etc. (is this correct?)

For example, the kernel for equation two is:

#include "eq66Im.h"
#include <cmath>

registerMooseObject("SCResApp", eq66Im);

InputParameters
eq66Im::validParams()
{
  InputParameters params = ADKernel::validParams();
  params.addClassDescription(
      "Kernel representing Im component of eq66 of bible. Variable is Psi_Re.");
  params.addCoupledVar("Psi_Im",
                       "Imaginary component of Psi");
  params.addCoupledVar("Phi",
                       "Phi");
  params.addParam<MaterialPropertyName>(
      "ucon", "ucon", "u from eq66");
  params.addParam<MaterialPropertyName>(
      "gamma", "gamma", "gamma from eq66");
  return params;
}

eq66Im::eq66Im(const InputParameters & parameters)
  : ADKernel(parameters),
    _Psi_Im(adCoupledValue("Psi_Im")),
    _grad_Psi_Im(adCoupledGradient("Psi_Im")),
    _Psi_Im_dot(adCoupledDot("Psi_Im")),
    _Phi(adCoupledValue("Phi")),
    _grad_Phi(adCoupledGradient("Phi")),
    _ucon(getADMaterialProperty<Real>("ucon")),
    _gamma(getADMaterialProperty<Real>("gamma"))
{
}

ADReal
eq66Im::computeQpResidual()
{

//    part1 = (_test[_i][_qp] * (-_ucon[_qp] * sqrt(1+ _gamma[_qp] * _gamma[_qp] * (_u[_qp] * _u[_qp] + _Psi_Im[_qp] * _Psi_Im[_qp])) ) * (_Psi_Im_dot[_qp] + _Phi[_qp] * _u[_qp]));

//    part2 = (_grad_test[_i][_qp] * (- _grad_Psi_Im[_qp]));

//    part3 = (_grad_test[_i][_qp] * (- _gamma[_qp] * _gamma[_qp] *  _u[_qp] * _grad_Phi[_qp]));

//    part4 = (_test[_i][_qp] * (1 - _u[_qp] * _u[_qp] - _Psi_Im[_qp] * _Psi_Im[_qp]) * _Psi_Im[_qp]);

  return  (_test[_i][_qp] * (-_ucon[_qp] * sqrt(1+ _gamma[_qp] * _gamma[_qp] * (_u[_qp] * _u[_qp] + _Psi_Im[_qp] * _Psi_Im[_qp])) ) * (_Psi_Im_dot[_qp] + _Phi[_qp] * _u[_qp])) + (_grad_test[_i][_qp] * (- _grad_Psi_Im[_qp])) + (_grad_test[_i][_qp] * (- _gamma[_qp] * _gamma[_qp] *  _u[_qp] * _grad_Phi[_qp])) +  (_test[_i][_qp] * (1 - _u[_qp] * _u[_qp] - _Psi_Im[_qp] * _Psi_Im[_qp]) * _Psi_Im[_qp]);

}

Does this look like I'm doing the right thing?
Is this the correct way to approach the problem?

Thank you

[GiudGiud]

Hello

As I have 4 variables, I have to pick a variable for each equation to solve for, but it is arbitrary which I choose, as long as each variable has at least one active kernel. (is this correct?)

Yes, but this choice is consequential. This will determine what is on-diagonal and off-diagonal in the linearized system matrices that you solve during the linear steps of the nonlinear solve.
A good idea is to make sure to have big diagonals for every variable.

The other variables in the equation need to be coupled in as adCoupledValue etc. (is this correct?)

correct.

Does this look like I'm doing the right thing?

I think the idea is right, though:

• i think the integration by parts for term 3 is not correct
• we usually split each term on a per-kernel basis. It just improves the code re-use.

To verify your implementation, I would set 3 of the 4 variables to known functions, and solve for the 4th by hand. Then make sure the solver can solve each of these single equations.
To set a variable, use an auxiliary variable and a FunctionIC

[gj19866]

Amazing thank you GiudGiud! What is wrong with the 3rd term in the equation, please?
Thank you for all the help!

[GiudGiud]

integration by parts of

psi_r div_grad_phi test
leads to
grad_phi * (grad_psi_r test + grad_test * psi_r)

and you only have
_grad_test[_i][_qp] * _u[_qp] * _grad_Phi[_qp])

so missing one term?

[gj19866]

Ah yes, I understand now.
I had tried implementing the weak form as in https://mooseframework.inl.gov/getting_started/examples_and_tutorials/tutorial01_app_development/step04_weak_form.html#:~:text=Procedure%20for%20Generating%20Weak%20Forms&text=Write%20down%20the%20strong%20form,rule%20to%20separate%20the%20terms but I see that I was performing the product rule incorrectly!

Thank you very much GiudGiud, you have made a really big difference for me! I really appreciate all your help!

[gj19866]

I have now corrected and implemented all of these, but it's really struggling to converge.
What is the best plan of attack to fix the convergence issues?

[GiudGiud]

This plan
https://mooseframework.inl.gov/moose/application_usage/failed_solves.html