[parthiv1991]

Hello Everyone,

I am trying to implement an Implicit Runge Kutta time integration scheme in OpenFOAM for incompressible flows. I am stuck with the implementation of the Jacobian in OpenFOAM.
I am actually implementing the Rosenbrock-Wanner method, which needs the (numerical) Jacobian at every time-step.

For example, in MATLAB I would do something like this to compute the Jacobian of a vector numerically:

epsilon = 1e-6; % delta
l_x0=length(x0); % length of x0;
f0=feval(f,x0,varargin{:}); % caclulate f0
l_f=size(f0,1); % check the size of f
% z = zeros(l_f,1);


for i=1:l_x0
    dx = [ zeros(i-1,1); epsilon; zeros(l_x0-i,1)];;
    df(:,i) = ( feval(f,x0+dx,varargin{:}) - f0)/epsilon;
end

Could anyone please suggest me how do I compute the Jacobian numerically in OpenFOAM using the volVectorField U ?? I am stuck with this problem for couple of months now and I couldn't find a solution. Any suggestion/example would be a great help for me as my project deadline is fast approaching.

Thank you,

Parthiv.

[ano]

Hello Parthiv,

I understand that you want to calculate the Jacobian/gradient of the velocity field:

It should be

fvc::grad(U)

Openfoam defines the gradient as the transpose of the Jacobian. So consider whether you have to multiply the matrix with a vector from the left or the right (I guess you have to multiply from the right).
You can find the definition in eqn 2.3 in the OF Programmers guide
(You can also find transpose there. However, I would recommend to transform your equations to OF style instead of using the Jacobian.)

[parthiv1991]

Thank you very much Ano,

I will implement it the way you have suggested and then let you know how it works out.

Regards,
Parthiv

[parthiv1991]

Hello Ano,

There is a second problem I am facing and that is with the Identity matrix in Openfoam. I have seen the posts related to it as implemented the code as suggested there:

How to access identity matrix in OF?

but I am getting lots of errors.

I have also converted the equations in openfoam format. But when I multiply the Jacobian, it shows errors like : no match for ‘operator*’.

I have implemented another solver based on ESDIRK method and that works fine. This ROSW is causing a lot of difficulties somehow.

I am new to openfoam and hence I am facing a lot of trouble while coding.

[ano]

Dear Parthiv,

at the first glance it doesn't seem to make sense to write fvc::grad(U)*IdentityMatrix which would be fvc::grad(U) again. perhaps you wanted to write fvc::grad(U)&IdentityMatrix or IdentityMatrix&fvc::grad(U)?

However, some more information is required. Could you please write the exact equation you want to implement, your implemented version of the code and also give the errors (including the variable types) you get?

PS: "no match for ‘operator*’" just tells you that a multiplication between the variable type on the left side of the operator and the one on the right side is not implemented in OF, probably because it is not required.

[parthiv1991]

Dear Ano,

The source of the problem is while calculating the slope "k" of the function.

I'm solving the NS Equations so the function(U) looks like this:

 volVectorField fu = (-fvc::div(phi,U) + nu*fvc::laplacian(U));

Now the ROSW method is written as:

(Identity - gamma*h*Jacobian) * k[i] = h* (function(U+alpha[i,j]*k[j]) + Jacobian*(g[i,j]*k[j]))

where gamma is a scalar, alpha[i,j] and g[i,j] are coefficients from butcher table.

So I need to solve the above equation, which is of the form Ax = b, to obtain the slope "k". Hence the coefficient matrix "A" and rhs vector "b" need to be assembled.
A = (Identity - gamma*h*Jacobian); b = h* (function(U+alpha[i,j]*k[j]) + Jacobian*(g[i,j]*k[j]))

So for the first step, the equation is: (Identity - gamma*h*Jacobian) * k[1] = h* function(U).

I use k as a "volVectorField", Jacobian as "volTensorField"

[ano]

Dear Parthiv,

you want an inner product before the k[i], use a & instead of * (see programmers guide).
And I strongly guess that you want to write "k[i]&(Identity - gamma*h*Jacobian)" (but I could be wrong).

Your grad(U) gives you a tensor, your k[1] has to consist of three values, i.e. be a vector to give you a vector (since your U is also a vector).

Why did you decide to implement this algorithm?

[parthiv1991]

Thank you so much Ano for your suggestions.
Yes, the slope k holds vector values, and causes no problems.

The problems are arising for the Jacobian and Identity matrix implementations.

Even declaring an Identity tensor is generating so many errors as I said in the last post too.

and well, my project is about implementing DIRK and ROSW solvers for incompressible flows in OpenFOAM, and I was not the one to decide the topic, hence I have to somehow implement it and finish it, my degree depends on it.

[parthiv1991]

Hello Ano,

I have attached below the log file for the compilation of the solver.

And also please check your private message for the solver.

I'm implementing the traditional Rosenbrock-Wanner method.