How to Solve a 4th-Order Partial Differential Equation in MOOSE

[zzzmx-josh]

Check these boxes if you have followed the posting rules.

• Q&A General is the most appropriate section for my question
• I have consulted the posting Guidelines on the Discussions front page
• I have searched the Discussions forum and my question has not been asked before
• I have searched the MOOSE website and the documentation does not answer my question
• I have formatted my post following the posting guidelines (screenshots as a last resort, triple back quotes around pasted text)

Question

Hi MOOSE team and community,

I am trying to solve a 4th-order partial differential equation of the form:

with periodic boundary conditions for 𝑢. I would like to ask if there is a straightforward way to solve such an equation in MOOSE? I understand that typically, this kind of high-order equation might be challenging to solve directly. Is the only viable approach to introduce intermediate variables and reduce the system into lower-order equations (for example, by introducing new variables for the second derivatives)?

Any guidance or suggestions on how to implement this in MOOSE would be greatly appreciated!

Thank you!

[GiudGiud]

Hello

I can't think of a direct way here

Is the only viable approach to introduce intermediate variables and reduce the system into lower-order equations (for example, by introducing new variables for the second derivatives)?

this should work. We have coupledSecond in the coupleable API to let you represent these terms.
start with a second order mesh (Mesh/second_order = true) and have at least 2nd order variables for the intermediate variable(s)

[zzzmx-josh]

Thank you for your previous response.
Based on the previous replies, I wrote two kernels, which are as follows:

D2udi2Uii.h

#ifndef D2udi2Uii_H
#define D2udi2Uii_H

#include "Kernel.h"

class D2udi2Uii : public Kernel
{
public:
  D2udi2Uii(const InputParameters & parameters);
  static InputParameters validParams();

protected:

  virtual Real computeQpResidual() override;
  virtual Real computeQpJacobian() override;

private:
  // The component (0: x, 1: y, 2: z) to apply this kernel to
  const unsigned int _ii;

  // Reference to the second derivative (or gradient coupling) of variable 'w'
  const VariableSecond & _second_grad_couple_w;
};

#endif

D2udi2Uii.C

#include "D2udi2Uii.h"

// Register the kernel with the application
registerMooseObject("babblerApp", D2udi2Uii);

// Method to define valid input parameters for this kernel
InputParameters
D2udi2Uii::validParams()
{
  InputParameters params = Kernel::validParams();
  
  // Add a description for this kernel
  params.addClassDescription("Split formulation D2udi2Uii for polar component Kernel");

  // Add required parameter to specify which component (x, y, z) this kernel applies to
  params.addRequiredParam<unsigned int>("component", "The component (0: x, 1: y, 2: z) to apply this kernel to.");
  params.addRequiredCoupledVar("w", "e.g., vx");

  return params;
}

// Constructor for D2udi2Uii
D2udi2Uii::D2udi2Uii(const InputParameters & parameters)
  : Kernel(parameters),
    _ii(getParam<unsigned int>("component")),  // Get the component (x, y, z)
    _second_grad_couple_w(coupledSecond("w"))  // Get second derivative of coupled variable 'w'
{
}

Real
D2udi2Uii::computeQpResidual()
{
  Real residual = (_u[_qp] - _second_grad_couple_w[_qp](_ii, _ii)) * _test[_i][_qp];
  return residual;
}

// Method to compute the Jacobian at a quadrature point
Real
D2udi2Uii::computeQpJacobian()
{
  Real jacobian = _phi[_j][_qp] * _test[_i][_qp];
  
  return jacobian;
}

SecondDerivative.h

#pragma once

#include "ADKernel.h"

/// Class for computing the second derivative kernel
class SecondDerivative : public ADKernel
{
public:
  // Method to define valid input parameters for this kernel
  static InputParameters validParams();

  SecondDerivative(const InputParameters & parameters);

protected:

  virtual ADReal computeQpResidual() override;
  virtual Real computeQpJacobian() override;

private:
  // Index for variable component (e.g., 0: x, 1: y, etc.)
  const unsigned int _ii;

  // Gradient of the coupled variable 'v'
  const VariableGradient & _grad_v;
};

SecondDerivative.C

#include "SecondDerivative.h"

// Register the SecondDerivative kernel with the application
registerMooseObject("babblerApp", SecondDerivative);

// Method to define valid input parameters for this ADKernel
InputParameters SecondDerivative::validParams()
{
  // Initialize parameters from the base ADKernel
  InputParameters params = ADKernel::validParams();

  // Add a required parameter to specify which component (x, y, z) this kernel applies to
  params.addRequiredParam<unsigned int>("component", "The component (0: x, 1: y, 2: z) to apply this ADKernel to.");

  params.addRequiredCoupledVar("variable_to_take_the_derivative_of", "The coupled velocity variable v.");

  return params;
}

// Constructor for SecondDerivative
SecondDerivative::SecondDerivative(const InputParameters &parameters)
    : ADKernel(parameters),
      _ii(getParam<unsigned int>("component")),   // Get the specified component (x, y, or z)
      _grad_v(coupledGradient("variable_to_take_the_derivative_of")) // Get the gradient of the coupled velocity variable
{
}

// Method to compute the residual at a quadrature point
ADReal SecondDerivative::computeQpResidual()
{
  return -_grad_v[_qp](_ii) * _grad_test[_i][_qp](_ii);
}

Real SecondDerivative::computeQpJacobian()
{
  return 0;
}

And the corresponding input file:

[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 20
  ny = 20
  nz = 2
  xmin = -20
  xmax = 20
  ymin = -20
  ymax = 20
  zmin = -1
  zmax = 1
  second_order = true
  # elem_type = HEX8
[]
[Variables]
  [u]
    order = SECOND
  []
  [d2udx2]
    order = SECOND
  []
  [d2udy2]
    order = SECOND
  []
  [d2udz2]
    order = SECOND
  []
[]
[Kernels]
  [u_xx]
    type = D2udi2Uii
    variable = d2udx2
    w = u
    component = 0
  []
  [u_yy]
    type = D2udi2Uii
    variable = d2udy2
    w = u
    component = 1
  []
  [u_zz]
    type = D2udi2Uii
    variable = d2udz2
    w = u
    component = 2
  []
  [d4udx4]
    type = SecondDerivative
    variable = u
    variable_to_take_the_derivative_of = d2udx2
    component = 0
   []
  [d4udy4]
    type = SecondDerivative
    variable = u
    variable_to_take_the_derivative_of = d2udy2
    component = 1
   []
  [d4udz4]
    type = SecondDerivative
    variable = u
    variable_to_take_the_derivative_of = d2udz2
    component = 2
   []
   [f]
    type = BodyForce
    variable = u 
    function = f
  []
[]

It looks like there might still be some issues in the code. Could you please help me identify and resolve them?
Your expertise would be greatly appreciated!

[GiudGiud]

For ADkernels you do not need to add the computeQpJacobian routine

[zzzmx-josh]

Thanks! I will remove the computeQpJacobian function from the ADKernel.
Aside from that, is everything else okay?

[GiudGiud]

Are you sure about taking the component in second derivative?
Now that you use these intermediate variable then do an integration by parts, you should see something like:

  return _grad_d2udx2[_qp] * _grad_test[_i][_qp];

[zzzmx-josh]

Thank you for pointing that out, I made a mistake there. Many thanks again for your help!