Stratimikos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
MixedOrderPhysicsBasedPreconditioner.cpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Stratimikos: Thyra-based strategies for linear solvers
5// Copyright (2006) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
38//
39// ***********************************************************************
40// @HEADER
41
43#include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
44#include "Thyra_PreconditionerFactoryHelpers.hpp"
45#include "Thyra_DefaultInverseLinearOp.hpp"
46#include "Thyra_DefaultMultipliedLinearOp.hpp"
47#include "Thyra_EpetraThyraWrappers.hpp"
48#include "Thyra_EpetraLinearOp.hpp"
49#include "Thyra_get_Epetra_Operator.hpp"
50#include "Thyra_TestingTools.hpp"
51#include "EpetraExt_CrsMatrixIn.h"
52#include "Epetra_CrsMatrix.h"
53#include "Teuchos_GlobalMPISession.hpp"
54#include "Teuchos_VerboseObject.hpp"
55#include "Teuchos_XMLParameterListHelpers.hpp"
56#include "Teuchos_CommandLineProcessor.hpp"
57#include "Teuchos_StandardCatchMacros.hpp"
58#include "Teuchos_TimeMonitor.hpp"
59#ifdef HAVE_MPI
60# include "Epetra_MpiComm.h"
61#else
62# include "Epetra_SerialComm.h"
63#endif
64
65
75namespace {
76
77
78// Read a Epetra_CrsMatrix from a Matrix market file
79Teuchos::RCP<Epetra_CrsMatrix>
80readEpetraCrsMatrixFromMatrixMarket(
81 const std::string fileName, const Epetra_Comm &comm
82 )
83{
84 Epetra_CrsMatrix *A = 0;
85 TEUCHOS_TEST_FOR_EXCEPTION(
86 0!=EpetraExt::MatrixMarketFileToCrsMatrix( fileName.c_str(), comm, A ),
87 std::runtime_error,
88 "Error, could not read matrix file \""<<fileName<<"\"!"
89 );
90 return Teuchos::rcp(A);
91}
92
93
94// Read an Epetra_CrsMatrix in as a wrapped Thyra::EpetraLinearOp object
95Teuchos::RCP<const Thyra::LinearOpBase<double> >
96readEpetraCrsMatrixFromMatrixMarketAsLinearOp(
97 const std::string fileName, const Epetra_Comm &comm,
98 const std::string &label
99 )
100{
101 return Thyra::epetraLinearOp(
102 readEpetraCrsMatrixFromMatrixMarket(fileName,comm),label
103 );
104}
105
106
107} // namespace
108
109
110int main(int argc, char* argv[])
111{
112
113 using Teuchos::describe;
114 using Teuchos::rcp;
115 using Teuchos::rcp_dynamic_cast;
116 using Teuchos::rcp_const_cast;
117 using Teuchos::RCP;
118 using Teuchos::CommandLineProcessor;
119 using Teuchos::ParameterList;
120 using Teuchos::sublist;
121 typedef ParameterList::PrintOptions PLPrintOptions;
122 using Thyra::inverse;
123 using Thyra::initializePreconditionedOp;
124 using Thyra::initializeOp;
125 using Thyra::unspecifiedPrec;
126 using Thyra::solve;
127 typedef RCP<const Thyra::LinearOpBase<double> > LinearOpPtr;
128 typedef RCP<Thyra::VectorBase<double> > VectorPtr;
129
130 bool success = true;
131 bool verbose = true;
132
133 Teuchos::GlobalMPISession mpiSession(&argc,&argv);
134
135 Teuchos::RCP<Teuchos::FancyOStream>
136 out = Teuchos::VerboseObjectBase::getDefaultOStream();
137
138 try {
139
140 //
141 // Read in options from the command line
142 //
143
144 const int numVerbLevels = 6;
145 Teuchos::EVerbosityLevel
146 verbLevelValues[] =
147 {
148 Teuchos::VERB_DEFAULT, Teuchos::VERB_NONE,
149 Teuchos::VERB_LOW, Teuchos::VERB_MEDIUM,
150 Teuchos::VERB_HIGH, Teuchos::VERB_EXTREME
151 };
152 const char*
153 verbLevelNames[] =
154 { "default", "none", "low", "medium", "high", "extreme" };
155
156 Teuchos::EVerbosityLevel verbLevel = Teuchos::VERB_MEDIUM;
157 std::string paramsFile = "";
158 std::string extraParamsFile = "";
159 std::string baseDir = ".";
160 bool useP1Prec = true;
161 bool invertP1 = false;
162 bool showParams = false;
163 double solveTol = 1e-8;
164
165 CommandLineProcessor clp(false); // Don't throw exceptions
166
167 clp.setOption( "verb-level", &verbLevel,
168 numVerbLevels, verbLevelValues, verbLevelNames,
169 "Verbosity level used for all objects."
170 );
171 clp.setOption( "params-file", &paramsFile,
172 "File containing parameters in XML format.",
173 true // Required
174 );
175 clp.setOption( "extra-params-file", &extraParamsFile,
176 "File containing extra parameters in XML format."
177 );
178 clp.setOption( "base-dir", &baseDir,
179 "Base directory for all of the files."
180 );
181 clp.setOption( "use-P1-prec", "use-algebraic-prec", &useP1Prec,
182 "Use the physics-based preconditioner for P2 based on P1 or just an algebraic preconditioner"
183 );
184 clp.setOption( "invert-P1", "prec-P1-only", &invertP1,
185 "In the physics-based preconditioner, use P1's preconditioner or fully invert P1."
186 );
187 clp.setOption( "solve-tol", &solveTol,
188 "Tolerance for the solution to determine success or failure!"
189 );
190 clp.setOption( "show-params", "no-show-params", &showParams,
191 "Show the parameter list or not."
192 );
193 clp.setDocString(
194 "This example program shows an attempt to create a physics-based\n"
195 "preconditioner using the building blocks of Stratimikos and Thyra\n"
196 "implicit linear operators. The idea is to use the linear operator\n"
197 "for first-order Lagrange finite elements as the preconditioner for the\n"
198 "operator using second-order Lagrange finite elements. The details of the\n"
199 "PDE being represented, the mesh being used, or the boundary conditions are\n"
200 "beyond the scope of this example.\n"
201 "\n"
202 "The matrices used in this example are:\n"
203 "\n"
204 " P2: The discretized PDE matrix using second-order Lagrange finite elements.\n"
205 " P1: The discretized PDE matrix using first-order Lagrange finite elements.\n"
206 " M22: The mass matrix for the second-order Lagrange finite-element basis functions.\n"
207 " M11: The mass matrix for the first-order Lagrange finite-element basis functions.\n"
208 " M21: A rectangular matrix that uses mixed first- and second-order basis functions\n"
209 " to map from P2 space to P1 space.\n"
210 " M12: A rectangular matrix that uses mixed first- and second-order basis functions\n"
211 " to map from P1 space to P2 space.\n"
212 "\n"
213 "The above matrices are read from Matrix Market *.mtx files in the directory given by\n"
214 "the --base-dir command-line option.\n"
215 "\n"
216 "The preconditioner operator created in this example program is:\n"
217 "\n"
218 " precP2Op = (inv(M22) * M12) * prec(P1) * (inv(M11) * M21)\n"
219 "\n"
220 "where prec(P1) is either the algebraic preconditioner for P1 (--prec-P1-only)\n"
221 "or is a full solve for P1 (--invert-P1).\n"
222 "\n"
223 "We use Stratimikos to specify the linear solvers and/or algebraic\n"
224 "preconditioners and we use the Thyra implicit operators to build the\n"
225 "implicitly multiplied linear operators associated with the preconditioner.\n"
226 "\n"
227 "Warning! This physics-based preconditioner is singular and can not\n"
228 "be used to solve the linear system given a random right-hand side. However,\n"
229 "singular or not, this example shows how one can use Thyra/Stratimikos to quickly\n"
230 "try out these types of preconditioning ideas (even if they do not work).\n"
231 );
232
233 // Note: Use --help on the command line to see the above documentation
234
235 CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse(argc,argv);
236 if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL ) return parse_return;
237
238
239 //
240 *out << "\nA) Reading in Epetra_CrsMatrix objects for P1, P2, M11, M12, M21, and M22 ...\n";
241 //
242
243#ifdef HAVE_MPI
244 Epetra_MpiComm comm(MPI_COMM_WORLD);
245#else
246 Epetra_SerialComm comm;
247#endif
248
249 LinearOpPtr P1=readEpetraCrsMatrixFromMatrixMarketAsLinearOp(
250 baseDir+"/P1.mtx",comm,"P1");
251 *out << "\nP1 = " << describe(*P1,verbLevel) << "\n";
252 LinearOpPtr P2= readEpetraCrsMatrixFromMatrixMarketAsLinearOp(
253 baseDir+"/P2.mtx",comm,"P2");
254 *out << "\nP2 = " << describe(*P2,verbLevel) << "\n";
255 LinearOpPtr M11=readEpetraCrsMatrixFromMatrixMarketAsLinearOp(
256 baseDir+"/M11.mtx",comm,"M11");
257 *out << "\nM11 = " << describe(*M11,verbLevel) << "\n";
258 LinearOpPtr M22=readEpetraCrsMatrixFromMatrixMarketAsLinearOp(
259 baseDir+"/M22.mtx",comm,"M22");
260 *out << "\nM22 = " << describe(*M22,verbLevel) << "\n";
261 LinearOpPtr M12=readEpetraCrsMatrixFromMatrixMarketAsLinearOp(
262 baseDir+"/M12.mtx",comm,"M12");
263 *out << "\nM12 = " << describe(*M12,verbLevel) << "\n";
264 LinearOpPtr M21=readEpetraCrsMatrixFromMatrixMarketAsLinearOp(
265 baseDir+"/M21.mtx",comm,"M21");
266 *out << "\nM21 = " << describe(*M21,verbLevel) << "\n";
267
268 // ToDo: Replace the above functions with a general Thyra strategy object
269 // to do the reading
270
271 //
272 *out << "\nB) Get the preconditioner and/or linear solver strategies to invert M11, M22, P1, and P2 ...\n";
273 //
274
275 //
276 // Get separate parameter sublists for each square operator separately
277 // that specify the type of linear solver and/or preconditioner to use.
278 //
279
280 RCP<ParameterList> paramList =
281 Teuchos::getParametersFromXmlFile( baseDir+"/"+paramsFile );
282 if (extraParamsFile.length()) {
283 Teuchos::updateParametersFromXmlFile( baseDir+"/"+extraParamsFile, paramList.ptr() );
284 }
285 if (showParams) {
286 *out << "\nRead in parameter list:\n\n";
287 paramList->print(*out,PLPrintOptions().indent(2).showTypes(true));
288 }
289
290 Stratimikos::DefaultLinearSolverBuilder M11_linsolve_strategy_builder;
291 M11_linsolve_strategy_builder.setParameterList(
292 sublist(paramList,"M11 Solver",true) );
293
294 Stratimikos::DefaultLinearSolverBuilder M22_linsolve_strategy_builder;
295 M22_linsolve_strategy_builder.setParameterList(
296 sublist(paramList,"M22 Solver",true) );
297
298 Stratimikos::DefaultLinearSolverBuilder P1_linsolve_strategy_builder;
299 P1_linsolve_strategy_builder.setParameterList(
300 sublist(paramList,"P1 Solver",true) );
301
302 Stratimikos::DefaultLinearSolverBuilder P2_linsolve_strategy_builder;
303 P2_linsolve_strategy_builder.setParameterList(
304 sublist(paramList,"P2 Solver",true) );
305
306 //
307 // Create the linear solver and/or preconditioner strategies
308 // (i.e. factories)
309 //
310
311 // For M11 and M22, we want full linear solver factories with embedded
312 // algebraic preconditioner factories.
313
314 RCP<const Thyra::LinearOpWithSolveFactoryBase<double> > M11_linsolve_strategy
315 = createLinearSolveStrategy(M11_linsolve_strategy_builder);
316
317 RCP<Thyra::LinearOpWithSolveFactoryBase<double> > M22_linsolve_strategy
318 = createLinearSolveStrategy(M22_linsolve_strategy_builder);
319
320 // For P1, we only want its preconditioner factory.
321
322 RCP<const Thyra::LinearOpWithSolveFactoryBase<double> > P1_linsolve_strategy;
323 RCP<const Thyra::PreconditionerFactoryBase<double> > P1_prec_strategy;
324 if(invertP1)
325 P1_linsolve_strategy
326 = createLinearSolveStrategy(P1_linsolve_strategy_builder);
327 else
328 P1_prec_strategy
329 = createPreconditioningStrategy(P1_linsolve_strategy_builder);
330
331 // For P2, we only want a linear solver factory. We will supply the
332 // preconditioner ourselves (that is the whole point of this example).
333
334 RCP<const Thyra::LinearOpWithSolveFactoryBase<double> > P2_linsolve_strategy
335 = createLinearSolveStrategy(P2_linsolve_strategy_builder);
336
337 //
338 *out << "\nC) Create the physics-based preconditioner! ...\n";
339 //
340
341 *out << "\nCreating inv(M11) ...\n";
342 LinearOpPtr invM11 = inverse(*M11_linsolve_strategy, M11);
343 *out << "\ninvM11 = " << describe(*invM11,verbLevel) << "\n";
344
345 *out << "\nCreating inv(M22) ...\n";
346 LinearOpPtr invM22 = inverse(*M22_linsolve_strategy, M22);
347 *out << "\ninvM22 = " << describe(*invM22,verbLevel) << "\n";
348
349 *out << "\nCreating prec(P1) ...\n";
350 LinearOpPtr invP1;
351 if(invertP1) {
352 *out << "\nCreating prec(P1) as a full solver ...\n";
353 invP1 = inverse(*P1_linsolve_strategy, P1);
354 }
355 else {
356 *out << "\nCreating prec(P1) as just an algebraic preconditioner ...\n";
357 RCP<Thyra::PreconditionerBase<double> >
358 precP1 = prec(*P1_prec_strategy,P1);
359 *out << "\nprecP1 = " << describe(*precP1,verbLevel) << "\n";
360 invP1 = precP1->getUnspecifiedPrecOp();
361 }
362 rcp_const_cast<Thyra::LinearOpBase<double> >(
363 invP1)->setObjectLabel("invP1"); // Cast to set label ...
364 *out << "\ninvP1 = " << describe(*invP1,verbLevel) << "\n";
365
366 LinearOpPtr P2ToP1 = multiply( invM11, M21 );
367 *out << "\nP2ToP1 = " << describe(*P2ToP1,verbLevel) << "\n";
368
369 LinearOpPtr P1ToP2 = multiply( invM22, M12 );
370 *out << "\nP1ToP2 = " << describe(*P1ToP2,verbLevel) << "\n";
371
372 LinearOpPtr precP2Op = multiply( P1ToP2, invP1, P2ToP1 );
373 *out << "\nprecP2Op = " << describe(*precP2Op,verbLevel) << "\n";
374
375 //
376 *out << "\nD) Setup the solver for P2 ...\n";
377 //
378
379 RCP<Thyra::LinearOpWithSolveBase<double> >
380 P2_lows = P2_linsolve_strategy->createOp();
381 if(useP1Prec) {
382 *out << "\nCreating the solver P2 using the specialized precP2Op\n";
383 initializePreconditionedOp<double>( *P2_linsolve_strategy, P2,
384 unspecifiedPrec(precP2Op), P2_lows.ptr());
385 }
386 else {
387 *out << "\nCreating the solver P2 using algebraic preconditioner\n";
388 initializeOp(*P2_linsolve_strategy, P2, P2_lows.ptr());
389 }
390 *out << "\nP2_lows = " << describe(*P2_lows, verbLevel) << "\n";
391
392 //
393 *out << "\nE) Solve P2 for a random RHS ...\n";
394 //
395
396 VectorPtr x = createMember(P2->domain());
397 VectorPtr b = createMember(P2->range());
398 Thyra::randomize(-1.0, +1.0, b.ptr());
399 Thyra::assign(x.ptr(), 0.0); // Must give an initial guess!
400
401 Thyra::SolveStatus<double>
402 solveStatus = solve<double>( *P2_lows, Thyra::NOTRANS, *b, x.ptr() );
403
404 *out << "\nSolve status:\n" << solveStatus;
405
406 *out << "\nSolution ||x|| = " << Thyra::norm(*x) << "\n";
407
408 if(showParams) {
409 *out << "\nParameter list after use:\n\n";
410 paramList->print(*out, PLPrintOptions().indent(2).showTypes(true));
411 }
412
413 //
414 *out << "\nF) Checking the error in the solution of r=b-P2*x ...\n";
415 //
416
417 VectorPtr P2x = Thyra::createMember(b->space());
418 Thyra::apply( *P2, Thyra::NOTRANS, *x, P2x.ptr() );
419 VectorPtr r = Thyra::createMember(b->space());
420 Thyra::V_VmV<double>(r.ptr(), *b, *P2x);
421
422 double
423 P2x_nrm = Thyra::norm(*P2x),
424 r_nrm = Thyra::norm(*r),
425 b_nrm = Thyra::norm(*b),
426 r_nrm_over_b_nrm = r_nrm / b_nrm;
427
428 bool result = r_nrm_over_b_nrm <= solveTol;
429 if(!result) success = false;
430
431 *out
432 << "\n||P2*x|| = " << P2x_nrm << "\n";
433
434 *out
435 << "\n||P2*x-b||/||b|| = " << r_nrm << "/" << b_nrm
436 << " = " << r_nrm_over_b_nrm << " <= " << solveTol
437 << " : " << Thyra::passfail(result) << "\n";
438
439 Teuchos::TimeMonitor::summarize(*out<<"\n");
440 }
441 TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose,std::cerr,success)
442
443 if (verbose) {
444 if(success) *out << "\nCongratulations! All of the tests checked out!\n";
445 else *out << "\nOh no! At least one of the tests failed!\n";
446 }
447
448 return ( success ? EXIT_SUCCESS : EXIT_FAILURE );
449
450}
int main(int argc, char *argv[])
LinearSolverBuilder< double > DefaultLinearSolverBuilder