Re: [DuMuX] Dune::BCRSMatrix

[Flemisch, Bernd, apl.Prof.Dr.]

In that case (numSurfaceReactions > numSurfaceSpecies), you have to go for a 
least-squares solution

min_{x in R^numSurfaceSpecies} || A*x - b ||

which can be achieved by solving

(A^t*A)*x = b

I have not seen a dedicated ISTL solver for this.


My indexing was wrong in the first place. You should modify to

Dune::FieldMatrix<Scalar, numSurfaceReactions, numSurfaceSpecies> A;
Dune::FieldMatrix<Scalar, numSurfaceSpecies, numSurfaceReactions> At;
for(int speciesIdx=0; speciesIdx<numSurfaceSpecies; speciesIdx++)
for(int reactionIdx=0; reactionIdx<numSurfaceReactions; 
reactionIdx++)
{
A[reactionIdx][speciesIdx] = ...
At[speciesIdx][reactionIdx] = ...
}
auto AtA = A.leftmultiply(At);

// right hand side
Dune::FieldVector<Scalar, numSurfaceReactions> b;
for(int reactionIdx=0; reactionIdx<numSurfaceReactions; reactionIdx++)
b[reactionIdx] = …

// solution
Dune::FieldVector<Scalar, numSurfaceSpecies> x;
AtA.solve(x, b);



Von: Dumux <dumux-boun...@listserv.uni-stuttgart.de> im Auftrag von 
georg.fut...@dlr.de <georg.fut...@dlr.de>
Gesendet: Montag, 11. Dezember 2017 10:10:20
An: dumux@listserv.uni-stuttgart.de
Betreff: Re: [DuMuX] Dune::BCRSMatrix

Hi Bernd,

Actually, the system might be overdetermined, so FieldMatrix solve() will not 
necessarily work. Do I have any options to solve such a system with ISTL?

Kind regards
Georg

Von: Dumux [mailto:dumux-boun...@listserv.uni-stuttgart.de] Im Auftrag von 
Flemisch, Bernd, apl.Prof.Dr.
Gesendet: Montag, 11. Dezember 2017 10:03
An: DuMuX User Forum
Betreff: Re: [DuMuX] Dune::BCRSMatrix


Hi Georg,



for the system that you describe, there is no need to wrap things in 
BCRSMatrices and BlockVectors. You can use directly FieldMatrices and -Vectors:

Dune::FieldMatrix<Scalar, numSurfaceReactions, numSurfaceSpecies> A;

for(int speciesIdx=0; speciesIdx<numSurfaceSpecies; speciesIdx++)

for(int reactionIdx=0; reactionIdx<numSurfaceReactions; 
reactionIdx++)

A[speciesIdx][reactionIdx] = ...



// right hand side

Dune::FieldVector<Scalar, numSurfaceReactions> b;
for(int reactionIdx=0; reactionIdx<numSurfaceReactions; reactionIdx++)

b[reactionIdx] = …



// solution

Dune::FieldVector<Scalar, numSurfaceSpecies> x;

A.solve(x, b);


I hope that numSurfaceReactions == numSurfaceSpecies?



You would only need the more complicated setting if your system size was not 
known at compile time.



Kind regards

Bernd


Von: Dumux 
<dumux-boun...@listserv.uni-stuttgart.de<mailto:dumux-boun...@listserv.uni-stuttgart.de>>
 im Auftrag von georg.fut...@dlr.de<mailto:georg.fut...@dlr.de> 
<georg.fut...@dlr.de<mailto:georg.fut...@dlr.de>>
Gesendet: Montag, 11. Dezember 2017 09:29:07
An: dumux@listserv.uni-stuttgart.de<mailto:dumux@listserv.uni-stuttgart.de>
Betreff: [DuMuX] Dune::BCRSMatrix

Hello Dumux,

I would like to solve a linear system of equations (A*x = b) using the Dune 
ISTL. To use the  Dumux implementation (seqsolverbackend.hh), 
Dune::FieldVectors need to be wrapped in Dune::Blockvectors and a 
Dune::FieldMatrix is wrapped in a Dune::BCRSMatrix. I am not familiar with 
these data types and the Dune documentation on the homepage is really 
confusing. Since the system is not too large, I guess I can solve everything in 
a single block? This is the code I use right now:

Dune::FieldMatrix<Scalar, numSurfaceReactions, numSurfaceSpecies> A;
for(int speciesIdx=0; speciesIdx<numSurfaceSpecies; speciesIdx++)
for(int reactionIdx=0; reactionIdx<numSurfaceReactions; 
reactionIdx++)
assign values…

// right hand side
Dune::BlockVector<Dune::FieldVector<Scalar, numSurfaceReactions> > b;
b.resize(1); // does this mean there is only one block containing a 
FieldVector ?
for(int reactionIdx=0; reactionIdx<numSurfaceReactions; reactionIdx++)
b[0][reactionIdx] = …

// solution
Dune::BlockVector<Dune::FieldVector<Scalar, numSurfaceSpecies> > x;
x.resize(1);

typedef typename Dune::FieldMatrix<Scalar, numSurfaceReactions, 
numSurfaceSpecies> MatrixBlock;
typedef typename Dune::BCRSMatrix ISTLMatrix;

ISTLMatrix A2(1, 1, ISTLMatrix::random); // does this mean there is only one 
block containing a FieldMatrix ?

for(int speciesIdx=0; speciesIdx<numSurfaceSpecies; speciesIdx++)
for(int reactionIdx=0; reactionIdx<numSurfaceReactions; 
reactionIdx++)
A2[0][0][reactionIdx][speciesIdx] = … HOW TO ASSIGN VALUES???

I get a segmentation fault from the last line. Is there better information on 
BCRS mat

Re: [DuMuX] Dune::BCRSMatrix

[Flemisch, Bernd, apl.Prof.Dr.]

Hi Georg,


for the system that you describe, there is no need to wrap things in 
BCRSMatrices and BlockVectors. You can use directly FieldMatrices and -Vectors:

Dune::FieldMatrix A;

for(int speciesIdx=0; speciesIdx b;

for(int reactionIdx=0; reactionIdx x;

A.solve(x, b);


I hope that numSurfaceReactions == numSurfaceSpecies?


You would only need the more complicated setting if your system size was not 
known at compile time.


Kind regards

Bernd


Von: Dumux  im Auftrag von 
georg.fut...@dlr.de 
Gesendet: Montag, 11. Dezember 2017 09:29:07
An: dumux@listserv.uni-stuttgart.de
Betreff: [DuMuX] Dune::BCRSMatrix

Hello Dumux,

I would like to solve a linear system of equations (A*x = b) using the Dune 
ISTL. To use the  Dumux implementation (seqsolverbackend.hh), 
Dune::FieldVectors need to be wrapped in Dune::Blockvectors and a 
Dune::FieldMatrix is wrapped in a Dune::BCRSMatrix. I am not familiar with 
these data types and the Dune documentation on the homepage is really 
confusing. Since the system is not too large, I guess I can solve everything in 
a single block? This is the code I use right now:

Dune::FieldMatrix A;
for(int speciesIdx=0; speciesIdx > b;
b.resize(1); // does this mean there is only one block containing a 
FieldVector ?
for(int reactionIdx=0; reactionIdx > x;
x.resize(1);

typedef typename Dune::FieldMatrix MatrixBlock;
typedef typename Dune::BCRSMatrix ISTLMatrix;

ISTLMatrix A2(1, 1, ISTLMatrix::random); // does this mean there is only one 
block containing a FieldMatrix ?

for(int speciesIdx=0; speciesIdx
www.DLR.de


___
Dumux mailing list
Dumux@listserv.uni-stuttgart.de
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux