OffDiagonalJacobian for divergence of RealVectorVariable

[mhelel]

Check these boxes if you have followed the posting rules.

• Q&A Navier Stokes is the most appropriate category 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 (avoid screenshots if possible, triple back quotes before/after pasted text, etc)

Question

I'm currently translating parts of the INS Code such that it uses one VectorVariable as velocity $\vec{u}$ instead of three seperate variables. This lead to a fundamental question how to implement the OffDiagJacobian. I know that there is AD, but I wanted to try to write the code myself. So here is the question:

To implement

$\nabla\cdot\vec{u}=0$

I have already written this code:

#include "INSMassVector.h"
#include "Function.h"
registerMooseObject("otterApp", INSMassVector);
InputParameters INSMassVector::validParams() {
    InputParameters params = Kernel::validParams();

    params.addClassDescription("This class computes the mass equation residual and Jacobian "
                             "contributions for the incompressible Navier-Stokes momentum "
                             "equation.");

    params.addParam<bool>("pspg", false, "Whether to perform PSPG stabilization of the mass equation");
    
    params.addRequiredCoupledVar("velocity", "Coupled velocity vector");

    return params;
}
INSMassVector::INSMassVector(const InputParameters & parameters) : 
    Kernel(parameters),
    _pspg(getParam<bool>("pspg")),

    _velocity_id(coupled("velocity")),
    _velocity(coupledVectorValue("velocity")),
    _velocity_var(*getVectorVar("velocity",0)),
    _vector_grad_phi(_assembly.gradPhi(_velocity_var)),

    _grad_velocity(coupledVectorGradient("velocity")),

    _strongPressureTerm(getMaterialProperty<RealVectorValue>("strongPressureTerm")),
    _strongViscousTerm(getMaterialProperty<RealVectorValue>("strongViscousTerm")),
    _forcingTerm(getMaterialProperty<RealVectorValue>("forcingTerm")),
    _tau(getMaterialProperty<Real>("tau"))    
{}
Real INSMassVector::computeQpResidual() {

    Real r = - _grad_velocity[_qp].tr() * _test[_i][_qp];

    if (_pspg) {
        r += computeQpPGResidual();
    }

    return r;
}
Real INSMassVector::computeQpJacobian() {
    Real r = 0;

    if (_pspg) {
        r += computeQpPGJacobian();
    }

    return r;
}

The residual should be implemented correctly and the function 'computeQpPGResidual' as well as 'computeQpPGJacobian' does not matter for now. Now comes the tricky part so I will outline what I understand until now:

The Vector $\vec{u}$ is represented by $\vec{u}=\sum_{j=1}^N u_j \vec{\phi_j}$ with $\vec{\phi_j} = (\phi_j^x, \phi_j^y, \phi_j^x)^T$. So when taking the derivative of the residual on gets

$$\frac{\partial R}{\partial u_j} = \frac{\partial}{\partial u_j} div\vec{u} = \frac{\partial}{\partial u_j}(\partial_x\sum u_j\phi_j^x + \partial_y\sum u_j\phi_j^y + \partial_z\sum u_j\phi_j^z) = \partial_x\phi_j^x + \partial_y\phi_j^y + \partial_z\phi_j^z = div\vec{\phi_j} = \mathbf{1} \cdot\nabla\vec{\phi_j}$$

And now comes the part where I'm not shure how to implement this. I had the understanding, that one treats such equations the same as if $\vec{\phi_j}$ would just be a scalar and the rest is done by the Assembly process. Is this correct? Still then, how would one implement the $\mathbf{1}$?

Additional information

Sorry for the weird formatting after some of the equal signs. I don't have a clue why this happens

[lindsayad]

What's the motivation for making another implementation, e.g. making the velocity variables a vector instead of component scalars? I've worked on both implementation styles and at this point I prefer the latter because it's more likely that you can re-use kernels (e.g. advecting momentum components/passive scalars/energy all with the same kernel, diffusing momentum components/passive scalars/energy all with the same kernel)

[lindsayad]

That only directly addresses part 1 of your question. However, I hope my suggestion to think of vector variables the same as scalars resolves questions 2 and 3. MOOSE does not think about them differently other than that it must use different data members for scalar and vector field shape functions because they are fundamentally different types (e.g. in the EFieldAdvection kernel _phi is of type VariablePhiValue whereas _vector_phi is a pointer to VectorVariablePhiValue)

[lindsayad]

If you still need help with questions 2 and 3 let me know.

In general, yes, I encourage you to leverage what's already in the navier_stokes module. It's less work for you and you'd have help from me (and others) if you wanted to fix or expand capability there

[mhelel]

Thank you a lot. I will just try a few more things and if a question comes up, I will ask you again.

[mhelel]

For Anyone interested: Moose does magic. To implement the OffDiagJacobian for the $div \vec{u}$ one can just use the trace of Matrix as proposed above. So the working code is

#include "INSMassVector.h"
#include "Function.h"

registerMooseObject("otterApp", INSMassVector);

InputParameters INSMassVector::validParams() {
    InputParameters params = Kernel::validParams();

    params.addClassDescription("This class computes the mass equation residual and Jacobian "
                             "contributions for the incompressible Navier-Stokes momentum "
                             "equation.");

    params.addRequiredCoupledVar("velocity", "Coupled velocity vector");

    return params;
}

INSMassVector::INSMassVector(const InputParameters & parameters) : 
    Kernel(parameters),

    _velocity_id(coupled("velocity")),
    _velocity(coupledVectorValue("velocity")),
    _velocity_var(*getVectorVar("velocity",0)),
    _vector_grad_phi(_assembly.gradPhi(_velocity_var)),

    _grad_velocity(coupledVectorGradient("velocity"))
    
{}

Real INSMassVector::computeQpResidual() {
    return _grad_velocity[_qp].tr() * _test[_i][_qp];
}

Real INSMassVector::computeQpJacobian() {
    return 0.0;
}

// This needs a little more review
Real INSMassVector::computeQpOffDiagJacobian(const unsigned int jvar) {
    if (jvar == _velocity_id) {
        return _vector_grad_phi[_j][_qp].tr() * _test[_i][_qp];
    }
    else return 0.0;
}