Sorry, the right hand side should be A^t*b:
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] = …
Dune::FieldVector<Scalar, numSurfaceSpecies> Atb;
A.mv(b, Atb);
// solution
Dune::FieldVector<Scalar, numSurfaceSpecies> x;
AtA.solve(x, Atb);
________________________________
Von: Dumux <[email protected]> im Auftrag von Flemisch,
Bernd, apl.Prof.Dr. <[email protected]>
Gesendet: Montag, 11. Dezember 2017 10:31:07
An: DuMuX User Forum
Betreff: Re: [DuMuX] Dune::BCRSMatrix
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 <[email protected]> im Auftrag von
[email protected] <[email protected]>
Gesendet: Montag, 11. Dezember 2017 10:10:20
An: [email protected]
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:[email protected]] 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
<[email protected]<mailto:[email protected]>>
im Auftrag von [email protected]<mailto:[email protected]>
<[email protected]<mailto:[email protected]>>
Gesendet: Montag, 11. Dezember 2017 09:29:07
An: [email protected]<mailto:[email protected]>
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<MatrixBlock> 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 matrices somewhere? Any hints on how to set up this matrix are greatly
appreciated!
Best regards
Georg Futter
——————————————————————————
German Aerospace Center (DLR)
Institute of Engineering Thermodynamics | Computational Electrochemistry |
Pfaffenwaldring 38-40 | 70569 Stuttgart
Dipl.-Ing. Georg Futter | Ph.D. student
Telefon 0711/6862-8135 | [email protected]<mailto:[email protected]>
www.DLR.de<http://www.dlr.de/>
_______________________________________________
Dumux mailing list
[email protected]
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
