[romant]

Hej,

I am looking for a way to solve the following equation

a_i = c1 * a_j * du_i / dx_j + c2

for i=1 this looks like

a_1 = c1 * a_1 * du_1/dx_1 + c1 * a_2 * du_1 / dx_2 + c1 * a_3 * du_1 / dx_3 + c2

where a is a vector, u is a vector and i and j run from 1-3. This can be solved by either matrix inversion where A_ij = (I- (c * du_i / dx_j)) is the coefficient matrix in every cell (I is the identity matrix).

The idea is to rewrite this into an explicit form, by taking a_1 to the other side.

I was also looking into a form to solve this with the build in fvMatrix of OF, unfortunately you can basically only insert div, grad (or other derivative forms) and source. However, source (fvm::Sp(scalar, solvableVariable))cannot take a tensor as it would be required from du_i/dx_j as a multiplication factor for vector a

The other idea is to take the coefficient matrix (A_ij) in every cell and divide the right hand side by this.

A_ij*a_j = S_i => a_j = A_ij^-1 * S_i

Does anybody know how to do this in OF?

OpenFOAM version used is 2.1.

[romant]

Further checking this, I found a function inv(const tensorField&), which I believe takes the inverse of a tensor, but I am unsure for what it takes the inverse, either the tensor in every cell or the overall tensor (the field itself).

I now assumed that it takes the inverse of every tensor in the cells and the following code came out for the following equation

(I-A_ij) * a_j = S_i
-> a_j = (I-A_ij)^1 * S_i

        // velocity gradient matrix, M_ij
      const volTensorField coeffMatrix(-sThetauCoeff*Cthetau2*fvc::grad(U));
        // I_ij - M_ij
      const volTensorField negCoeffMatrix(I-coeffMatrix);
        // (I_ij - M_ij)^1
      const volTensorField invNegCoeffMatrix(inv(negCoeffMatrix));

      thetau = invNegCoeffMatrix & (gSource+gradTSource);

[romant]

The solution was a little bit easier than first assumed. OpenFOAM has a great operator definition which lets you divide by a tensor

So in case you have an equation

A_ij x_j = s_i

and you want to solve for x_j, where A is a second order tensor (for example the gradient of velocity, grad(U) ) and x and s are vectors, the one simply divides by A_ij

as in

  volVectorField s(g);  // g is the gravity field
  volTensorField A(fvc::grad(U));
  volVectorField x = s / A;

I guess you could do the same with

 volVectorField s(g); // g is the gravity field
 volTensorField A(fvc::grad(U));
 volTensorField invA(inv(A));  // taking the inverse of A
 volVectorField x = (1/det(A)*invA) & s; // det(A) takes the determinant

the previous solution is of course much simpler and more elegant

[atoof]

Dear Roman,

I need to calculate a vector using anothor vector and a matrix by formula:

a=[X^T X]^-1 X^T y,

where y is a vector and X is a square matrix (T means transpose and -1 means inverse of a matrix), but the matrix and vectors are not volScalarField.

Is there any possible to do that by openFOAM?

Regards,

[romant]

Originally Posted by atoof

Dear Roman,

I need to calculate a vector using anothor vector and a matrix by formula:

a=[X^T X]^-1 X^T y,

where y is a vector and X is a square matrix (T means transpose and -1 means inverse of a matrix), but the matrix and vectors are not volScalarField.

Is there any possible to do that by openFOAM?

Regards,

X should be a volTensorField I think and then y is a volVectorField, there is an operator for transpose, which is X.T() for a tensor X, so an inverse (^1) you get with inv(X) and a transpose with X.T()