new Kernel: Porous Flow' style QpJac (Hardly converged) Vs. General hand-coding QpJac (easily converged) #27907

Bearunner · 2024-06-15T19:49:28Z

Bearunner
Jun 15, 2024

Dear all,

Recently, I am learning the dirackernels/PorousFlowPeacemanBorehole and try to implement it in the kernels.

Initially, I finish it by using general hand-coding QpJac and get good converged results.

Recently, I try to use the Porous Flow' style QpJac (e.g. PorousFlowMassTimeDerivative.C and PorousFlowFullySaturatedMassTimeDerivative.C ) to speed up the simulation.

1 General hand-coding QpJac (easily convergence).txt
2 Porous Flow' style QpJac (Hardly convergence).txt
1_log.txt
2_log.txt

Unlucky, I find it is hard to converge even at first step.

Could you please help me on this?

Thanks a lot.

B

General hand-coding QpJac (easily converged)

Real
PorousFlowPeaceman::computeQpResidual()
{
    Real r = 0.0;

    for (unsigned int ph = 0; ph < _num_phases; ++ph)
    {

     r += _fluid_density_node[_i][_ph] / _fluid_viscosity[_i][_ph] * (_u[_i] - 100000));

    }

  return r * _test[_i][_qp];
}


Real
PorousFlowPeaceman::computeQpJacobian()
{
  return computeQpJac(_var.number());
}



Real
PorousFlowPeaceman::computeQpJac(unsigned int jvar) const
{
  const unsigned int pvar = _dictator.porousFlowVariableNum(jvar);

  Real j = 0.0;

  if (jvar == _p_var_number)

  {
    for (unsigned int ph = 0; ph < _num_phases; ++ph)
  {

  j += -_dfluid_viscosity_dvar[_i][ph][pvar] * _phi[_j][_qp] * _fluid_density_node[_i][ph] / std::pow(_fluid_viscosity[_i][ph], 2) + _dfluid_density_node_dvar[_i][ph][pvar] * _phi[_j][_qp] / _fluid_viscosity[_i][ph];
  j *= (_u[_i] - 100000));
  j += _fluid_density_node[_i][ph] / _fluid_viscosity[_i][ph] * _phi[_j][_qp];

  }

  }

  return j * _test[_i][_qp];
}

Real
PorousFlowPeaceman::computeQpOffDiagJacobian(unsigned int jvar)
{
  const unsigned int pvar = _dictator.porousFlowVariableNum(jvar);

  Real j = 0.0;


  if (jvar == _T_var_number)

  {

  for (unsigned int ph = 0; ph < _num_phases; ++ph)

  {

  j += -_dfluid_viscosity_dvar[_i][ph][pvar] * _phi[_j][_qp] * _fluid_density_node[_i][ph] / std::pow(_fluid_viscosity[_i][ph], 2) + _dfluid_density_node_dvar[_i][ph][pvar] * _phi[_j][_qp] / _fluid_viscosity[_i][ph];
  j *= (_u[_i] - 100000));

  }

  }

  return j * _test[_i][_qp];
}

log:easily converged

Time Step 1, time = 0.1, dt = 0.1
    |residual|_2 of individual variables:
               porepressure: 1.76996
 0 Nonlinear |R| = 1.769956e+00
Currently Executing
    Computing Jacobian..........                                                         [ 56.95 s] [  547 MB]
      0 Linear |R| = 1.769956e+00
      1 Linear |R| = 5.037185e-03
      2 Linear |R| = 2.008511e-05
      3 Linear |R| = 5.927711e-08
      4 Linear |R| = 1.849117e-10
      5 Linear |R| = 3.081674e-13
    |residual|_2 of individual variables:
               porepressure: 0.000207102
 1 Nonlinear |R| = 2.071023e-04
    Computing Jacobian.........                                                          [ 54.19 s] [  896 MB]
      0 Linear |R| = 2.071023e-04
      1 Linear |R| = 7.806158e-08
      2 Linear |R| = 1.620506e-10
      3 Linear |R| = 2.768883e-13
      4 Linear |R| = 4.960755e-16
      5 Linear |R| = 1.366996e-18
    |residual|_2 of individual variables:
               porepressure: 2.08483e-05
 2 Nonlinear |R| = 2.084832e-05
    Computing Jacobian.........                                                          [ 53.74 s] [  897 MB]
      0 Linear |R| = 2.084832e-05
      1 Linear |R| = 7.685754e-09
      2 Linear |R| = 1.321019e-11
      3 Linear |R| = 2.502334e-14
      4 Linear |R| = 4.155741e-17
      5 Linear |R| = 8.306119e-20
    |residual|_2 of individual variables:
               porepressure: 1.8536e-06
 3 Nonlinear |R| = 1.853598e-06
    Computing Jacobian.........                                                          [ 54.96 s] [  898 MB]
      0 Linear |R| = 1.853598e-06
      1 Linear |R| = 6.889821e-10
      2 Linear |R| = 1.292019e-12
      3 Linear |R| = 2.389580e-15
      4 Linear |R| = 5.196694e-18
      5 Linear |R| = 2.334065e-20
    |residual|_2 of individual variables:
               porepressure: 1.7901e-07
 4 Nonlinear |R| = 1.790101e-07
 Solve Converged!
  Finished Solving                                                                       [228.06 s] [  559 MB]

Porous Flow' style QpJac (Hardly converged)


Real
PorousFlowPeaceman::computeQpResidual()
{
    Real r = 0.0;


     r += _fluid_density_node[_i][_ph] / _fluid_viscosity[_i][_ph] * (_pp[_i][_ph] - 100000) * _mass_frac[_i][_ph][_fluid_component];

  return r * _test[_i][_qp];
}


Real
PorousFlowPeaceman::computeQpJacobian()
{
  if (!_var_is_porflow_var)
  return 0.0;

  return computeQpJac(_dictator.porousFlowVariableNum(_var.number()));
}


Real
PorousFlowPeaceman::computeQpOffDiagJacobian(unsigned int jvar)
{
  if (_dictator.notPorousFlowVariable(jvar))
    return 0.0;
  return computeQpJac(_dictator.porousFlowVariableNum(jvar));
}


Real
PorousFlowPeaceman::computeQpJac(unsigned int pvar)
{
  Real j = 0.0;

  j += -_dfluid_viscosity_dvar[_i][_ph][pvar] * _fluid_density_node[_i][_ph] / std::pow(_fluid_viscosity[_i][_ph], 2) + _dfluid_density_node_dvar[_i][_ph][pvar] / _fluid_viscosity[_i][_ph];
  j *= (_pp[_i][_ph] - 100000) * _mass_frac[_i][_ph][_fluid_component];
  j += _fluid_density_node[_i][_ph] / _fluid_viscosity[_i][_ph] * (_dpp_dvar[_i][_ph][pvar] * _mass_frac[_i][_ph][_fluid_component] + (_pp[_i][_ph] - 100000) * _dmass_frac_dvar[_i][_ph][_fluid_component][pvar]);

  return j * _test[_i][_qp];
}

log: Hardly converged

Time Step 1, time = 0.1, dt = 0.1

|residual|_2 of individual variables:
           porepressure: 1.1471
0 Nonlinear |R| = 1.147099e+00
Currently Executing
   Computing Jacobian...........                                                        [ 64.59 s] [  467 MB]
     0 Linear |R| = 1.147099e+00
     1 Linear |R| = 1.444729e-02
     2 Linear |R| = 3.728843e-04
     3 Linear |R| = 5.364687e-06
     4 Linear |R| = 1.140004e-07
     5 Linear |R| = 7.141013e-10
     6 Linear |R| = 3.131963e-12
     7 Linear |R| = 7.391971e-15
   |residual|_2 of individual variables:
              porepressure: 0.166439
1 Nonlinear |R| = 1.664389e-01
   Computing Jacobian...........                                                        [ 63.57 s] [  761 MB]
     0 Linear |R| = 1.664389e-01
     1 Linear |R| = 2.626463e-03
     2 Linear |R| = 3.870046e-05
     3 Linear |R| = 8.534643e-07
     4 Linear |R| = 1.141531e-08
     5 Linear |R| = 9.185115e-11
     6 Linear |R| = 4.429794e-13
     7 Linear |R| = 1.007129e-15
   |residual|_2 of individual variables:
              porepressure: 0.0262875
2 Nonlinear |R| = 2.628752e-02
   Computing Jacobian...........                                                        [ 62.79 s] [  762 MB]
     0 Linear |R| = 2.628752e-02
     1 Linear |R| = 4.084056e-04
     2 Linear |R| = 1.022684e-05
     3 Linear |R| = 3.412311e-07
     4 Linear |R| = 4.277030e-09
     5 Linear |R| = 1.550830e-11
     6 Linear |R| = 6.803209e-14
     7 Linear |R| = 1.427683e-16
   |residual|_2 of individual variables:
              porepressure: 0.00664934
3 Nonlinear |R| = 6.649337e-03
   Computing Jacobian...........                                                        [ 62.28 s] [  765 MB]
     0 Linear |R| = 6.649337e-03
     1 Linear |R| = 1.908789e-04
     2 Linear |R| = 5.361577e-06
     3 Linear |R| = 1.166733e-07
     4 Linear |R| = 7.627529e-10
     5 Linear |R| = 2.290047e-12
     6 Linear |R| = 9.519533e-15
     7 Linear |R| = 2.726405e-17
   |residual|_2 of individual variables:
              porepressure: 0.0030829
4 Nonlinear |R| = 3.082901e-03
   Computing Jacobian...........                                                        [ 62.97 s] [  763 MB]
     0 Linear |R| = 3.082901e-03
     1 Linear |R| = 6.883939e-05
     2 Linear |R| = 1.860284e-06
     3 Linear |R| = 5.889309e-08
     4 Linear |R| = 2.524422e-10
     5 Linear |R| = 6.137773e-13
     6 Linear |R| = 2.665876e-15
   |residual|_2 of individual variables:
              porepressure: 0.00163936
5 Nonlinear |R| = 1.639360e-03
   Computing Jacobian...........                                                        [ 62.40 s] [  765 MB]
     0 Linear |R| = 1.639360e-03
     1 Linear |R| = 3.855174e-05
     2 Linear |R| = 1.124493e-06
     3 Linear |R| = 3.161580e-08
     4 Linear |R| = 1.246819e-10
     5 Linear |R| = 2.108129e-13
     6 Linear |R| = 4.080236e-16
   |residual|_2 of individual variables:
              porepressure: 0.000882297
6 Nonlinear |R| = 8.822973e-04
   Computing Jacobian...........                                                        [ 62.45 s] [  766 MB]
     0 Linear |R| = 8.822973e-04
     1 Linear |R| = 2.033103e-05
     2 Linear |R| = 5.844302e-07
     3 Linear |R| = 1.696500e-08
     4 Linear |R| = 6.698822e-11
     5 Linear |R| = 1.217664e-13
     6 Linear |R| = 3.566407e-16
   |residual|_2 of individual variables:
              porepressure: 0.000475276

7 Nonlinear |R| = 4.752755e-04
   Computing Jacobian...........                                                        [ 62.40 s] [  766 MB]
     0 Linear |R| = 4.752755e-04
     1 Linear |R| = 1.100907e-05
     2 Linear |R| = 3.179114e-07
     3 Linear |R| = 9.145917e-09
     4 Linear |R| = 3.603855e-11
     5 Linear |R| = 6.383199e-14
     6 Linear |R| = 1.672315e-16
   |residual|_2 of individual variables:
              porepressure: 0.000256038
8 Nonlinear |R| = 2.560380e-04
   Computing Jacobian...........                                                        [ 62.47 s] [  766 MB]
     0 Linear |R| = 2.560380e-04
     1 Linear |R| = 5.922062e-06
     2 Linear |R| = 1.708009e-07
     3 Linear |R| = 4.925906e-09
     4 Linear |R| = 1.941879e-11
     5 Linear |R| = 3.462604e-14
     6 Linear |R| = 9.367789e-17
   |residual|_2 of individual variables:
              porepressure: 0.000137932
9 Nonlinear |R| = 1.379321e-04
   Computing Jacobian...........                                                        [ 62.07 s] [  767 MB]
     0 Linear |R| = 1.379321e-04
     1 Linear |R| = 3.191600e-06
     2 Linear |R| = 9.208255e-08
     3 Linear |R| = 2.653845e-09
     4 Linear |R| = 1.046051e-11
     5 Linear |R| = 1.861778e-14
     6 Linear |R| = 4.993557e-17
   |residual|_2 of individual variables:
              porepressure: 7.43065e-05
10 Nonlinear |R| = 7.430646e-05
   Computing Jacobian...........                                                        [ 62.06 s] [  765 MB]
     0 Linear |R| = 7.430646e-05
     1 Linear |R| = 1.719180e-06
     2 Linear |R| = 4.959614e-08
     3 Linear |R| = 1.429646e-09
     4 Linear |R| = 5.635367e-12
     5 Linear |R| = 1.003496e-14
     6 Linear |R| = 2.697880e-17
   |residual|_2 of individual variables:
              porepressure: 4.00302e-05
11 Nonlinear |R| = 4.003021e-05
   Computing Jacobian...........                                                        [ 62.03 s] [  767 MB]
     0 Linear |R| = 4.003021e-05
     1 Linear |R| = 9.261817e-07
     2 Linear |R| = 2.671992e-08
     3 Linear |R| = 7.701803e-10
     4 Linear |R| = 3.035851e-12
     5 Linear |R| = 5.405229e-15
     6 Linear |R| = 1.452236e-17
   |residual|_2 of individual variables:
              porepressure: 2.1565e-05
12 Nonlinear |R| = 2.156496e-05
   Computing Jacobian...........                                                        [ 62.31 s] [  767 MB]
     0 Linear |R| = 2.156496e-05
     1 Linear |R| = 4.989457e-07
     2 Linear |R| = 1.439422e-08
     3 Linear |R| = 4.149085e-10
     4 Linear |R| = 1.635471e-12
     5 Linear |R| = 2.912018e-15
     6 Linear |R| = 7.825182e-18
   |residual|_2 of individual variables:
              porepressure: 1.16174e-05
13 Nonlinear |R| = 1.161744e-05
   Computing Jacobian...........                                                        [ 62.77 s] [  767 MB]
     0 Linear |R| = 1.161744e-05
     1 Linear |R| = 2.687917e-07
     2 Linear |R| = 7.754469e-09
     3 Linear |R| = 2.235190e-10
     4 Linear |R| = 8.810599e-13
     5 Linear |R| = 1.568759e-15
     6 Linear |R| = 4.215468e-18
   |residual|_2 of individual variables:
              porepressure: 6.25852e-06
14 Nonlinear |R| = 6.258517e-06
   Computing Jacobian...........                                                        [ 62.96 s] [  767 MB]
     0 Linear |R| = 6.258517e-06
     1 Linear |R| = 1.448026e-07
     2 Linear |R| = 4.177452e-09
     3 Linear |R| = 1.204134e-10
     4 Linear |R| = 4.746418e-13
     5 Linear |R| = 8.450914e-16
     6 Linear |R| = 2.270897e-18
15 Nonlinear |R| = 3.371568e-06
   Computing Jacobian...........                                                        [ 63.14 s] [  768 MB]
     0 Linear |R| = 3.371568e-06
     1 Linear |R| = 7.800762e-08
     2 Linear |R| = 2.250466e-09
     3 Linear |R| = 6.486885e-11
     4 Linear |R| = 2.556965e-13
     5 Linear |R| = 4.552742e-16
     6 Linear |R| = 1.223370e-18
   |residual|_2 of individual variables:
              porepressure: 1.81632e-06
16 Nonlinear |R| = 1.816321e-06
   Computing Jacobian...........                                                        [ 62.88 s] [  768 MB]
     0 Linear |R| = 1.816321e-06
     1 Linear |R| = 4.202415e-08
     2 Linear |R| = 1.212373e-09
     3 Linear |R| = 3.494600e-11
     4 Linear |R| = 1.377448e-13
     5 Linear |R| = 2.452573e-16
     6 Linear |R| = 6.590093e-19
   |residual|_2 of individual variables:
              porepressure: 9.78507e-07

17 Nonlinear |R| = 9.785069e-07
Solve Converged!
 Finished Solving                                                                       [1102.69 s] [  428 MB]

Answered by cpgr

Jun 17, 2024

OK, thinking some more and I think your derivatives are correct for the case where you are using nodal properties. I think you just need to add something like

  if (_i != _j)
    return 0.0;

at the start of QpJac() for the same reason that PorousFlowMassTimeDerivative does (where it has a comment to remind us As the fluid mass is lumped to the nodes, only non-zero terms are for _i==_j). So the derivative of density, mass fraction, should be 0 if _i != _j

View full answer

GiudGiud · 2024-06-15T20:41:46Z

GiudGiud
Jun 15, 2024
Collaborator

Hello

You may want to follow the advice on convergence issues here:
https://mooseframework.inl.gov/moose/application_usage/failed_solves.html

If you suspect the Jacobian is not correct, then there are techniques on that page to check and suggestions on how to remedy this

7 replies

cpgr Jun 17, 2024
Collaborator

Is it just simply a typo somewhere in the derivative when you include mass fractions? That seems to be the only difference?

cpgr Jun 17, 2024
Collaborator

Also, doing it the second way won't speed up the simulation I wouldn't think? If the first way is working, maybe use that (but it looks like you probably want to include mass fractions in that one?)

cpgr Jun 17, 2024
Collaborator

But just realised the first one assumes that the variable if pressure, where it might not be, which is why we used the second form. The second one is missing _phi in there isn't it?

cpgr Jun 17, 2024
Collaborator

OK, thinking some more and I think your derivatives are correct for the case where you are using nodal properties. I think you just need to add something like

  if (_i != _j)
    return 0.0;

at the start of QpJac() for the same reason that PorousFlowMassTimeDerivative does (where it has a comment to remind us As the fluid mass is lumped to the nodes, only non-zero terms are for _i==_j). So the derivative of density, mass fraction, should be 0 if _i != _j

Answer selected by Bearunner

Bearunner Jun 17, 2024
Author

OK, thinking some more and I think your derivatives are correct for the case where you are using nodal properties. I think you just need to add something like
  if (_i != _j)
    return 0.0;
at the start of QpJac() for the same reason that PorousFlowMassTimeDerivative does (where it has a comment to remind us As the fluid mass is lumped to the nodes, only non-zero terms are for _i==_j). So the derivative of density, mass fraction, should be 0 if _i != _j

Thanks for your tips, Chris.

I add it as follows.

Real
PorousFlowPeaceman::computeQpJac(unsigned int pvar)
{
Real j = 0.0;

if (_i != _j)
return 0.0;

j += -_dfluid_viscosity_dvar[_i][_ph][pvar] * _fluid_density_node[_i][_ph] / std::pow(_fluid_viscosity[_i][_ph], 2) + _dfluid_density_node_dvar[_i][_ph][pvar] / _fluid_viscosity[_i][_ph];
j *= (_pp[_i][_ph] - 100000) * _mass_frac[_i][_ph][_fluid_component];
j += _fluid_density_node[_i][_ph] / _fluid_viscosity[_i][_ph] * (_dpp_dvar[_i][_ph][pvar] * _mass_frac[_i][_ph][_fluid_component] + (_pp[_i][_ph] - 100000) * _dmass_frac_dvar[_i][_ph][_fluid_component][pvar]);

return j * _test[_i][_qp];
}

But I got a negative mass_fraction of one component at 5000 seconds around the injection position.

It looks like a [_qp] or [_i] error problem?

Could you please help me on this?

Thanks a lot.

Bearunner Jun 17, 2024
Author

Hello

You may want to follow the advice on convergence issues here: https://mooseframework.inl.gov/moose/application_usage/failed_solves.html

If you suspect the Jacobian is not correct, then there are techniques on that page to check and suggestions on how to remedy this

Thanks for your suggestions, Guillaume.

Bearunner Jun 18, 2024
Author

Hello Chris,

I learned the lines 293-296 in PorousFlowPeacemanBorehole.C and Lines 341-350 in PorousFlowLineSink.C

So I changed the '_pp[_i][_ph]' into _pp[_qp][ph]. Until now, it's working well without negative mass_fraction as above picture.

Could you please help me confirm whether this change is correct or not?

Thanks a lot.

Real
PorousFlowPeaceman::computeQpResidual()
{
    Real r = 0.0;

    for (unsigned ph = 0; ph < _num_phases; ++ph)
    {
     r += _fluid_density_node[_i][ph] / _fluid_viscosity[_i][ph] * (_pp[_qp][ph] - 100000) * _mass_frac[_i][ph][_fluid_component];
    }

  return r * _test[_i][_qp];
}


Real
PorousFlowPeaceman::computeQpJacobian()
{
  if (!_var_is_porflow_var)
  return 0.0;

  return computeQpJac(_dictator.porousFlowVariableNum(_var.number()));
}


Real
PorousFlowPeaceman::computeQpOffDiagJacobian(unsigned int jvar)
{
  // If the variable is not a PorousFlow variable, the OffDiag Jacobian terms are 0
  if (_dictator.notPorousFlowVariable(jvar))
    return 0.0;
  return computeQpJac(_dictator.porousFlowVariableNum(jvar));
}


Real
PorousFlowPeaceman::computeQpJac(unsigned int pvar)
{
Real j = 0.0;

for (unsigned ph = 0; ph < _num_phases; ++ph)
  {
  j += _fluid_density_node[_i][ph] / _fluid_viscosity[_i][ph] * _dpp_dvar[_qp][ph][pvar] * _phi[_j][_qp] * _mass_frac[_i][ph][_fluid_component];
  }

if (_i != _j)
  return j * _test[_i][_qp];

for (unsigned ph = 0; ph < _num_phases; ++ph)
  {
  j += (_pp[_qp][ph] - 100000) * _mass_frac[_i][ph][_fluid_component] * (_dfluid_density_node_dvar[_i][ph][pvar] / _fluid_viscosity[_i][ph] - _dfluid_viscosity_dvar[_i][ph][pvar] * _fluid_density_node[_i][ph] / std::pow(_fluid_viscosity[_i][ph], 2));
  j += (_pp[_qp][ph] - 100000) * _dmass_frac_dvar[_i][ph][_fluid_component][pvar] * _fluid_density_node[_i][ph] / _fluid_viscosity[_i][ph];
  }

  return j * _test[_i][_qp];
}

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

new Kernel: Porous Flow' style QpJac (Hardly converged) Vs. General hand-coding QpJac (easily converged) #27907

{{title}}

{{editor}}'s edit

{{editor}}'s edit

Replies: 1 comment 7 replies

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

Select a reply

new Kernel: Porous Flow' style QpJac (Hardly converged) Vs. General hand-coding QpJac (easily converged) #27907

Bearunner Jun 15, 2024

Replies: 1 comment · 7 replies

GiudGiud Jun 15, 2024 Collaborator

cpgr Jun 17, 2024 Collaborator

cpgr Jun 17, 2024 Collaborator

cpgr Jun 17, 2024 Collaborator

cpgr Jun 17, 2024 Collaborator

Bearunner Jun 17, 2024 Author

Bearunner Jun 17, 2024 Author

Bearunner Jun 18, 2024 Author

Bearunner
Jun 15, 2024

Replies: 1 comment 7 replies

GiudGiud
Jun 15, 2024
Collaborator

cpgr Jun 17, 2024
Collaborator

cpgr Jun 17, 2024
Collaborator

cpgr Jun 17, 2024
Collaborator

cpgr Jun 17, 2024
Collaborator

Bearunner Jun 17, 2024
Author

Bearunner Jun 17, 2024
Author

Bearunner Jun 18, 2024
Author