CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Programming & Development

inverse of a matrix or solution of implicit equation

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree2Likes
  • 2 Post By romant

Reply
 
LinkBack Thread Tools Display Modes
Old   April 25, 2012, 08:51
Default inverse of a matrix or solution of implicit equation
  #1
Senior Member
 
romant's Avatar
 
Roman Thiele
Join Date: Aug 2009
Location: Stockholm, Sweden
Posts: 359
Rep Power: 11
romant is on a distinguished road
Hej,

I am looking for a way to solve the following equation

Code:
a_i = c1 * a_j * du_i / dx_j + c2
for i=1 this looks like

Code:
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.

Code:
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.
__________________
~roman
romant is offline   Reply With Quote

Old   April 25, 2012, 10:35
Default
  #2
Senior Member
 
romant's Avatar
 
Roman Thiele
Join Date: Aug 2009
Location: Stockholm, Sweden
Posts: 359
Rep Power: 11
romant is on a distinguished road
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

Code:
        // 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);
__________________
~roman
romant is offline   Reply With Quote

Old   April 27, 2012, 05:21
Default solution found
  #3
Senior Member
 
romant's Avatar
 
Roman Thiele
Join Date: Aug 2009
Location: Stockholm, Sweden
Posts: 359
Rep Power: 11
romant is on a distinguished road
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

Code:
  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

Code:
 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
Tushar@cfd and huangxianbei like this.
__________________
~roman

Last edited by romant; April 27, 2012 at 05:32. Reason: addition
romant is offline   Reply With Quote

Old   August 6, 2013, 08:02
Default
  #4
Member
 
Hossein
Join Date: Apr 2010
Posts: 62
Rep Power: 7
atoof is on a distinguished road
Send a message via Yahoo to atoof
Dear Roman,

I need to calculate a vector using anothor vector and a matrix by formula:
Code:
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,
atoof is offline   Reply With Quote

Old   August 6, 2013, 09:02
Default
  #5
Senior Member
 
romant's Avatar
 
Roman Thiele
Join Date: Aug 2009
Location: Stockholm, Sweden
Posts: 359
Rep Power: 11
romant is on a distinguished road
Quote:
Originally Posted by atoof View Post
Dear Roman,

I need to calculate a vector using anothor vector and a matrix by formula:
Code:
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()
__________________
~roman
romant is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
error message cuteapathy CFX 14 March 20, 2012 07:45
Force can not converge colopolo CFX 13 October 4, 2011 22:03
Exact solution of Buckley_Leveret equation mcaro Main CFD Forum 0 January 25, 2011 07:52
Iterative solution to Poissons Equation Phiper Main CFD Forum 4 December 10, 2010 02:58
Conceptual trouble--please help me understand what my matrix solution is telling me bzz77 Main CFD Forum 0 March 25, 2010 17:31


All times are GMT -4. The time now is 04:37.