Conservative Advection with coupled variable - upwind scheme

[ngrilli]

Dear MOOSE Users,

I would like to solve this equation where \rho is the variable and q and p are coupled variables.
v is a velocity that I can import in the kernel as a MaterialProperty.

I would like a kernel for the terms on the right hand-side that has the upwind scheme.
I mean, I want the integral over the whole volume of the terms on the right side to be zero (in case of Neumann BC = zero flux).
The same as ConservativeAdvection in the framework, but with a coupled variable instead.

Is there anything ready for this in the moose modules?

If not, I would like to create my own kernel starting from ConservativeAdvection
and I will need to access the nodal values of the coupled variable to make the upwind scheme.
In the ConservativeAdvection kernel, this is done by:

_u_nodal(_var.dofValues()),

How do I do the same for a coupled variable? Is there anything similar in moose?

Thank you very much in advance for your help.

Best Regards,
Nicolò Grilli
National University of Singapore

[aeslaughter]

@lindsayad Do we have an example of this sort of upwind scheme?

[lindsayad]

Not with coupled variables. @ngrilli you could do:

_q_nodal(getVar("q", 0)->dofValues()),
_p_nodal(getVar("p", 0)->dofValues())

getVar is a method inherited from the Coupleable interface. The arguments to getVar are the name of the coupled variable (whatever the name is in your params.addCoupledVar("name", "doc_string")) and then the variable component, which is usually 0. But we do allow users to pass a vector of variables to a single coupled var parameter (this is not well advertised or known).

[ngrilli]

Dear @aeslaughter and @lindsayad

Thank you very much, this is very helpful.

I will develop my own kernel and will made it available in my github.

Best Regards,
Nicolò Grilli
National University of Singapore

[ngrilli]

Dear MOOSE Users,

I have developed the following kernel for the conservative advection of a coupled variable:
https://github.com/ngrilli/c_pfor_am/blob/main/src/kernels/ConservativeAdvectionCoupled.C

It seems to work when I use the full upwind scheme for my simulations.
However, I am not sure if the out of diagonal Jacobian is correct in the case of full upwind scheme.

Specifically, in:
ConservativeAdvectionCoupled::fullUpwind(JacRes res_or_jac)
I have used the nodal values of the coupled variable instead of the nodal value of the variable used in ConservativeAdvection

However, I am not sure if the following lines are storing the out of diagonal Jacobian properly
(and I have the same doubt for the residual):

if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
{
accumulateTaggedLocalMatrix();

if (_has_diag_save_in)
{
  unsigned int rows = _local_ke.m();
  DenseVector<Number> diag(rows);
  for (unsigned int i = 0; i < rows; i++)
    diag(i) = _local_ke(i, i);

  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
  for (const auto & var : _diag_save_in)
    var->sys().solution().add_vector(diag, var->dofIndices());
}

}

Best Regards,
Nicolò Grilli
National University of Singapore

[lindsayad]

The pasted lines of code, whether right or wrong, will not actually have any impact on the accuracy of your matrix being fed to PETSc. These lines are usually only for visualization purposes; perhaps you already knew that?

You can test the accuracy of your Jacobian using the command line options: -snes_test_jacobian -snes_test_jacobian_view. I would only run this test for smallish problems, e.g. 1000 dofs or less, ideally 30 dofs or less in order to truly understand the output from -snes_test_jacobian_view.

Are you worried about the speed of your simulation? If you do not care about a factor of 1.5 or so, I would highly recommend doing this with automatic differentiation, and then you will have no concern about the accuracy of your Jacobian.