
[Sponsors] 
April 25, 2012, 08:51 
inverse of a matrix or solution of implicit equation

#1 
Senior Member
Roman Thiele
Join Date: Aug 2009
Location: Stockholm, Sweden
Posts: 366
Rep Power: 12 
Hej,
I am looking for a way to solve the following equation Code:
a_i = c1 * a_j * du_i / dx_j + c2 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 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 OpenFOAM version used is 2.1.
__________________
~roman 

April 25, 2012, 10:35 

#2 
Senior Member
Roman Thiele
Join Date: Aug 2009
Location: Stockholm, Sweden
Posts: 366
Rep Power: 12 
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 (IA_ij) * a_j = S_i > a_j = (IA_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(IcoeffMatrix); // (I_ij  M_ij)^1 const volTensorField invNegCoeffMatrix(inv(negCoeffMatrix)); thetau = invNegCoeffMatrix & (gSource+gradTSource);
__________________
~roman 

April 27, 2012, 05:21 
solution found

#3 
Senior Member
Roman Thiele
Join Date: Aug 2009
Location: Stockholm, Sweden
Posts: 366
Rep Power: 12 
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; 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
__________________
~roman Last edited by romant; April 27, 2012 at 05:32. Reason: addition 

August 6, 2013, 08:02 

#4 
Member

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, Is there any possible to do that by openFOAM? Regards, 

August 6, 2013, 09:02 

#5  
Senior Member
Roman Thiele
Join Date: Aug 2009
Location: Stockholm, Sweden
Posts: 366
Rep Power: 12 
Quote:
__________________
~roman 

Thread Tools  
Display Modes  


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 troubleplease help me understand what my matrix solution is telling me  bzz77  Main CFD Forum  0  March 25, 2010 17:31 