How to solve ODEs with spatial variables coupled

[keassyguang]

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 everyone,
I am new to the MOOSE framework, and I am dealing with a set of ODEs with spatial variables. As shown in the example below, x is the spatial position, y is the coordinate in Y direction. m, n, and q are the constants. I want to solve B(x,t) and C(x,t) with given initializations as these two variables will be further used in my material models. I have tried to read the tutorials on the MOOSE website, but I cannot find any relevant examples. Can anyone tell me which Kernels should I use? It will be a great help if someone can show me some examples of how to solve it. Many thanks.

A(x,t) = A(x,t=0)exp(-y);
dB(x,t) / dt = mA(x,t)B(x,t) + nB(x,t)C(x,t);
dC(x,t) / dt = qB(x,t)C(x,t)

[GiudGiud]

Hello

ODEs in moose are solved by:

• scalar kernels for a single ODE with a single degree of freedom for the entire mesh
• nodal kernels for ODEs that live on every node of the mesh
  I think you are in the second case here.
  Use
  https://mooseframework.inl.gov/source/nodalkernels/TimeDerivativeNodalKernel.html for dB/dt
  then modify https://mooseframework.inl.gov/source/nodalkernels/CoupledForceNodalKernel.html to create mAB and nBC terms. Right now this kernel only supports something like mA

Do you want to solve for A? It seems like it's entirely predetermined from its initial condition?

[keassyguang]

if you subsistute in equation4, why do you need adCoupledDot?

Sorry for some typo. The equation 4 should be:
Equation 4: α(x,t) = 1 - C(x,t) / C(x,t=0)

[GiudGiud]

Ah ok.

So you do need to use adCoupledDot but you can use it on C imo

[keassyguang]

Thanks a lot. It works.

[keassyguang]

Hi @GiudGiud,
Thank you for your kind help. Could you please help me get the code below to display properly?
I hope this discussion will be helpful to others.

Here is the header files, source files and input file for the problem discussed:

A(x,t) = A0exp(-y);
dB(x,t) / dt = mA(x,t)B(x,t) + nB(x,t)C(x,t);
dC(x,t) / dt = qB(x,t)C(x,t)

// .h file for A(x,t)

#pragma once

#include "ADMaterial.h"

class AMaterial : public ADMaterial
{
public:
  static InputParameters validParams();
  AMaterial(const InputParameters & parameters);

  virtual void computeQpProperties() override;

protected:
  const Real _A0;

  ADMaterialProperty<Real> & _A;                  
};

// .C file for A(x,t)

#include "AMaterial.h"

registerMooseObject("testApp", AMaterial);

InputParameters
AMaterial::validParams()
{
  InputParameters params = ADMaterial::validParams();
  params.addRequiredParam<Real>("A0", "A0 constant");
  params.addClassDescription("Kernel to compute A(x,t)");

  return params;
}

AMaterial::AMaterial(const InputParameters & parameters)
  : ADMaterial(parameters),
    _A0(getParam<Real>("A0")),
    _A(declareADProperty<Real>("A"))
{
}

void AMaterial::computeQpProperties()
{
  Real y = _q_point[_qp](1);
  _A[_qp] =  _A0 * std::exp(-y);
}

// .h file for B(x,t)

#pragma once

#include "ADKernelValue.h"

class BKernel : public ADKernelValue
{
public:
  static InputParameters validParams();
  BKernel(const InputParameters & parameters);

protected:
  virtual ADReal precomputeQpResidual() override;

  const ADVariableValue & _C;

  const ADMaterialProperty<Real> & _A;    

  const Real _m;
  const Real _n;
};

//

.C file for B(x,t)

#include "BKernel.h"

registerMooseObject("testApp", BKernel);

InputParameters
BKernel::validParams()
{
  InputParameters params = ADKernelValue::validParams();
  params.addRequiredCoupledVar("C", "C variable");
  params.addRequiredParam<MaterialPropertyName>("A", "A property");
  params.addRequiredParam<Real>("m", "m parameter");
  params.addRequiredParam<Real>("n", "n parameter");
  params.addClassDescription("Kernel to compute right hand side of B(x,t)");
  return params;
}

BKernel::BKernel(const InputParameters & parameters)
  : ADKernelValue(parameters),
    _C(adCoupledValue("C")),
    _A(getADMaterialProperty<Real>("A")),
    _m(getParam<Real>("m")),
    _n(getParam<Real>("n"))
{
}

ADReal
BKernel::precomputeQpResidual()
{
  return - (_m * _A[_qp] * _u[_qp] + _n * _u[_qp] * _C[_qp]);
}

// .h file for C(x,t)

#pragma once

#include "ADKernelValue.h"

class CKernel : public ADKernelValue
{
public:
  static InputParameters validParams();
  CKernel(const InputParameters & parameters);

protected:
  virtual ADReal precomputeQpResidual() override;

  const ADVariableValue & _B;

  const Real _q;
};

// .C file for C(x,t)

#include "CKernel.h"

registerMooseObject("testApp", CKernel);

InputParameters
BKernel::validParams()
{
  InputParameters params = ADKernelValue::validParams();
  params.addRequiredCoupledVar("B", "B variable");
  params.addRequiredParam<Real>("q", "q parameter");
  params.addClassDescription("Kernel to compute right hand side of C(x,t)");
  return params;
}

CKernel::CKernel(const InputParameters & parameters)
  : ADKernelValue(parameters),
    _B(adCoupledValue("B")),
    _q(getParam<Real>("q"))
{
}

ADReal
BKernel::precomputeQpResidual()
{
  return -_q * _B[_qp] * _u[_qp];
}

input file example

[Mesh]
  # define your mesh
[]

[Variables]
  [B]
    order = FIRST
    family = LAGRANGE
    initial_condition = B_initial_value
  []

  [C]
    order = FIRST
    family = LAGRANGE
    initial_condition = C_initial_value
  []
[]

[AuxVariables]
  [A_aux] # for visulization
    order = FIRST
    family = MONOMIAL
  []

  [bounds_dummy] # for bounds control
    order = FIRST
    family = LAGRANGE
  []
[]


[Kernels]
  [B_timeDerivative]
    type = ADTimeDerivative
    variable = B
  []
  [B_right_hand_side]
    type = Bkernel
    variable = B
    A = A
    C = C
    m = m
    n = n
  []

  [C_timeDerivative]
    type = ADTimeDerivative
    variable = C
  []
  [C_right_hand_side]
    type = Ckernel
    variable = C
    B = B
    q = q
  []
[]


[AuxKernels]
  [A_aux]
    type = ADMaterialRealAux
    variable = A_aux
    property = A
    execute_on = 'initial timestep_end'
  []
[]


[Materials]
  [A_Material]
    type = AMaterial
    A0 = xxx
  []
[]


[Bounds] # Bounds are quite important at some times
  [B_upper_bound]
    type = ConstantBounds
    variable = bounds_dummy
    bounded_variable = B
    bound_type = upper
    bound_value = B_upper
  []
  [B_lower_bound]
    type = ConstantBounds
    variable = bounds_dummy
    bounded_variable = B
    bound_type = lower
    bound_value = B_lower
  []

  [C_upper_bound]
    type = ConstantBounds
    variable = bounds_dummy
    bounded_variable = C
    bound_type = upper
    bound_value = C_upper
  []
  [C_lower_bound]
    type = ConstantBounds
    variable = bounds_dummy
    bounded_variable = C
    bound_type = lower
    bound_value = C_lower
  []
[]


[Executioner]
# define your executioner
[]

[Outputs]
  exodus = true
[]

[GiudGiud]

you have to use triple back-quotes before and after the code to make it display nicely.

[keassyguang]

Thank you very much.