epetra/1Dfem/cxx_main.cpp

Go to the documentation of this file.
00001 //@HEADER
00002 // ***********************************************************************
00003 //
00004 //                     Rythmos Package
00005 //                 Copyright (2006) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Todd S. Coffey (tscoffe@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 //@HEADER
00028 
00029 #ifdef HAVE_MPI
00030 #include "Epetra_MpiComm.h"
00031 #include "mpi.h"
00032 #else
00033 #include "Epetra_SerialComm.h"
00034 #endif // HAVE_MPI
00035 
00036 #include "Epetra_Map.h"
00037 #include "Epetra_Vector.h"
00038 #include "Epetra_Version.h"
00039 
00040 #include "ExampleApplication.hpp"
00041 
00042 // Includes for Rythmos:
00043 #include "Rythmos_ConfigDefs.h"
00044 //#include "ExampleApplicationRythmosInterface.hpp"
00045 #include "Rythmos_ForwardEulerStepper.hpp"
00046 #include "Rythmos_BackwardEulerStepper.hpp"
00047 #include "Rythmos_ExplicitRKStepper.hpp"
00048 #include "Rythmos_ImplicitBDFStepper.hpp"
00049 
00050 // Includes for Thyra:
00051 #include "Thyra_EpetraThyraWrappers.hpp"
00052 #include "Thyra_EpetraLinearOp.hpp"
00053 #include "Thyra_EpetraModelEvaluator.hpp"
00054 #include "Thyra_TimeStepNewtonNonlinearSolver.hpp"
00055 #include "Thyra_DiagonalEpetraLinearOpWithSolveFactory.hpp"
00056 #include "Thyra_TestingTools.hpp"
00057 
00058 // Includes for Stratimikos:
00059 #ifdef HAVE_RYTHMOS_STRATIMIKOS
00060 #  include "Thyra_DefaultRealLinearSolverBuilder.hpp"
00061 #endif
00062 
00063 #include <string>
00064 
00065 // Includes for Teuchos:
00066 #include "Teuchos_RefCountPtr.hpp"
00067 #include "Teuchos_ParameterList.hpp"
00068 #include "Teuchos_CommandLineProcessor.hpp"
00069 #include "Teuchos_FancyOStream.hpp"
00070 #include "Teuchos_GlobalMPISession.hpp"
00071 #include "Teuchos_VerboseObject.hpp"
00072 #include "Teuchos_StandardCatchMacros.hpp"
00073 
00074 enum EMethod { METHOD_FE, METHOD_BE, METHOD_ERK, METHOD_BDF };
00075 enum STEP_METHOD { STEP_METHOD_FIXED, STEP_METHOD_VARIABLE };
00076 
00077 int main(int argc, char *argv[])
00078 {
00079   bool verbose = true; // verbosity level.
00080   bool result, success = true; // determine if the run was successfull
00081 
00082   Teuchos::RefCountPtr<Teuchos::FancyOStream>
00083     out = Teuchos::VerboseObjectBase::getDefaultOStream();
00084 
00085   try { // catch exceptions
00086 
00087 #ifdef HAVE_MPI
00088     MPI_Init(&argc,&argv);
00089     MPI_Comm mpiComm = MPI_COMM_WORLD;
00090     int procRank = 0;
00091     int numProc;
00092     MPI_Comm_size( mpiComm, &numProc );
00093     MPI_Comm_rank( mpiComm, &procRank );
00094 #endif // HAVE_MPI
00095 
00096     STEP_METHOD step_method = STEP_METHOD_FIXED;
00097     int numElements = 201; // number of elements in vector
00098     double finalTime = 1.0; // ODE final time
00099     int N = 100;  // number of steps to take
00100     const int num_methods = 4;
00101     const EMethod method_values[] = { METHOD_FE, METHOD_BE, METHOD_ERK, METHOD_BDF };
00102     const char * method_names[] = { "FE", "BE", "ERK", "BDF" };
00103     EMethod method_val = METHOD_BE;
00104     double maxError = 0.01;
00105     bool version = false;  // display version information 
00106     double reltol = 1.0e-2;
00107     double abstol = 1.0e-4;
00108     int maxOrder = 5;
00109 #ifdef Rythmos_DEBUG
00110     int debugLevel = 2; // debugLevel is used when Rythmos_DEBUG ifdef is set.
00111 #endif // Rythmos_DEBUG
00112 
00113     // Parse the command-line options:
00114     Teuchos::CommandLineProcessor  clp(false); // Don't throw exceptions
00115     clp.addOutputSetupOptions(true);
00116 #ifdef HAVE_RYTHMOS_STRATIMIKOS
00117     Thyra::DefaultRealLinearSolverBuilder lowsfCreator;
00118     lowsfCreator.setupCLP(&clp);
00119 #endif
00120 
00121     clp.setOption( "T", &finalTime, "Final time for simulation." );
00122     clp.setOption( "numelements", &numElements, "Problem size");
00123     clp.setOption( "method", &method_val, num_methods, method_values, method_names, "Integration method" );
00124     clp.setOption( "numsteps", &N, "Number of integration steps to take" );
00125     clp.setOption( "maxerror", &maxError, "Maximum error" );
00126     clp.setOption( "reltol", &reltol, "Relative Error Tolerance" );
00127     clp.setOption( "abstol", &abstol, "Absolute Error Tolerance" );
00128     clp.setOption( "maxorder", &maxOrder, "Maximum Implicit BDF order" );
00129     clp.setOption( "verbose", "quiet", &verbose, "Set if output is printed or not" );
00130     clp.setOption( "version", "run", &version, "Version of this code" );
00131 #ifdef Rythmos_DEBUG
00132     clp.setOption( "debuglevel", &debugLevel, "Debug Level for Rythmos" );
00133 #endif // Rythmos_DEBUG
00134 
00135 
00136     Teuchos::CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse(argc,argv);
00137     if( parse_return != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL ) return parse_return;
00138 
00139 #ifdef HAVE_RYTHMOS_STRATIMIKOS
00140     lowsfCreator.readParameters(out.get());
00141     *out << "\nThe parameter list after being read in:\n";
00142     lowsfCreator.getParameterList()->print(*out,2,true,false);
00143 #endif
00144 
00145 
00146     if (version) // Display version information and exit.
00147     {
00148       *out << Rythmos::Rythmos_Version() << std::endl; 
00149       *out << "basicExample Version 0.1 - 06/23/05" << std::endl;
00150       return(0);
00151     }
00152 
00153     if (finalTime <= 0.0)
00154     {
00155       std::cerr << "Final simulation time must be > 0.0." << std::endl;
00156       return(1);
00157     }
00158 
00159     
00160     // Set up the parameter list for the application:
00161     Teuchos::ParameterList params;
00162     params.set( "NumElements", numElements );
00163 #ifdef HAVE_MPI
00164     Teuchos::RefCountPtr<Epetra_Comm> epetra_comm_ptr_ = Teuchos::rcp( new Epetra_MpiComm(mpiComm) );
00165 #else
00166     Teuchos::RefCountPtr<Epetra_Comm> epetra_comm_ptr_ = Teuchos::rcp( new Epetra_SerialComm  );
00167 #endif // HAVE_MPI
00168 
00169     // Create the factory for the LinearOpWithSolveBase object
00170     Teuchos::RefCountPtr<Thyra::LinearOpWithSolveFactoryBase<double> >
00171       W_factory;
00172     if((method_val == METHOD_BE) or (method_val == METHOD_BDF))
00173     {
00174       //W_factory = Teuchos::rcp(new Thyra::DiagonalEpetraLinearOpWithSolveFactory());
00175       //W_factory = Teuchos::rcp(new Thyra::AmesosLinearOpWithSolveFactory());
00176 #ifdef HAVE_RYTHMOS_STRATIMIKOS
00177       W_factory = lowsfCreator.createLinearSolveStrategy("");
00178       *out
00179         << "\nCreated a LinearOpWithSolveFactory described as:\n"
00180         << Teuchos::describe(*W_factory,Teuchos::VERB_MEDIUM);
00181 #endif
00182 
00183     }
00184 
00185     // create interface to problem
00186     Teuchos::RefCountPtr<ExampleApplication1Dfem>
00187       epetraModel = Teuchos::rcp(new ExampleApplication1Dfem(epetra_comm_ptr_,params));
00188     Teuchos::RefCountPtr<Thyra::ModelEvaluator<double> >
00189       model = Teuchos::rcp(new Thyra::EpetraModelEvaluator(epetraModel,W_factory));
00190 
00191     // Create Stepper object depending on command-line input
00192     std::string method;
00193     Teuchos::RefCountPtr<Rythmos::Stepper<double> > stepper_ptr;
00194     if ( method_val == METHOD_ERK ) {
00195       stepper_ptr = Teuchos::rcp(new Rythmos::ExplicitRKStepper<double>(model));
00196       method = "Explicit Runge-Kutta of order 4";
00197     } else if (method_val == METHOD_FE) {
00198       stepper_ptr = Teuchos::rcp(new Rythmos::ForwardEulerStepper<double>(model));
00199       method = "Forward Euler";
00200     } else if (method_val == METHOD_BE) {
00201       Teuchos::RefCountPtr<Thyra::NonlinearSolverBase<double> >
00202         nonlinearSolver;
00203       Teuchos::RefCountPtr<Thyra::TimeStepNewtonNonlinearSolver<double> >
00204         _nonlinearSolver = Teuchos::rcp(new Thyra::TimeStepNewtonNonlinearSolver<double>());
00205       _nonlinearSolver->defaultTol(1e-3*maxError);
00206       nonlinearSolver = _nonlinearSolver;
00207       stepper_ptr = Teuchos::rcp(new Rythmos::BackwardEulerStepper<double>(model,nonlinearSolver));
00208       method = "Backward Euler";
00209     } else if (method_val == METHOD_BDF) {
00210       Teuchos::RefCountPtr<Thyra::NonlinearSolverBase<double> >
00211         nonlinearSolver;
00212       Teuchos::RefCountPtr<Thyra::TimeStepNewtonNonlinearSolver<double> >
00213         _nonlinearSolver = Teuchos::rcp(new Thyra::TimeStepNewtonNonlinearSolver<double>());
00214       _nonlinearSolver->defaultTol(1e-3*maxError);
00215       nonlinearSolver = _nonlinearSolver;
00216       Teuchos::ParameterList BDFparams;
00217       BDFparams.set( "stopTime", finalTime );
00218       BDFparams.set( "maxOrder", maxOrder );
00219       BDFparams.set( "relErrTol", reltol );
00220       BDFparams.set( "absErrTol", abstol );
00221 #ifdef Rythmos_DEBUG
00222       BDFparams.set( "debugLevel", debugLevel );
00223 #endif // Rythmos_DEBUG
00224       stepper_ptr = Teuchos::rcp(new Rythmos::ImplicitBDFStepper<double>(model,nonlinearSolver,BDFparams));
00225       step_method = STEP_METHOD_VARIABLE;
00226       method = "Implicit BDF";
00227     } else {
00228       TEST_FOR_EXCEPT(true);
00229     }
00230     Rythmos::Stepper<double> &stepper = *stepper_ptr;
00231 
00232 #ifdef Rythmos_DEBUG
00233     Teuchos::RefCountPtr<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
00234     if (debugLevel > 1)
00235     {
00236       stepper.describe(*out,Teuchos::VERB_EXTREME);
00237     }
00238 #endif // Rythmos_DEBUG
00239 
00240     double t0 = 0.0;
00241     double t1 = finalTime;
00242     double dt = (t1-t0)/N;
00243     double time = t0;
00244     int numSteps = 0;
00245 
00246     if (step_method == STEP_METHOD_FIXED)
00247     {
00248       // Integrate forward with fixed step sizes:
00249       for (int i=1 ; i<=N ; ++i)
00250       {
00251         double dt_taken = stepper.TakeStep(dt);
00252         numSteps++;
00253         if (dt_taken != dt)
00254         {
00255           cerr << "Error, stepper took step of dt = " << dt_taken << " when asked to take step of dt = " << dt << std::endl;
00256           break;
00257         }
00258         time += dt_taken;
00259       }
00260     }
00261     else // step_method == STEP_METHOD_VARIABLE
00262     {
00263       while (time < finalTime)
00264       {
00265         double dt_taken = stepper.TakeStep();
00266         numSteps++;
00267 #ifdef Rythmos_DEBUG
00268         if (debugLevel > 1)
00269         {
00270           stepper.describe(*out,Teuchos::VERB_EXTREME);
00271         }
00272 #endif // Rythmos_DEBUG
00273         if (dt_taken < 0)
00274         {
00275           cerr << "Error, stepper failed for some reason with step taken = " << dt_taken << endl;
00276           break;
00277         }
00278         time += dt_taken;
00279         *out << "Took stepsize of: " << dt_taken << " time = " << time << endl;
00280       }
00281     }
00282     *out << "Integrated to time = " << time << endl;
00283     // Get solution out of stepper:
00284     Teuchos::RefCountPtr<const Thyra::VectorBase<double> > x_computed_thyra_ptr = stepper.get_solution();
00285     // Convert Thyra::VectorBase to Epetra_Vector
00286     Teuchos::RefCountPtr<const Epetra_Vector>
00287       x_computed_ptr = Thyra::get_Epetra_Vector(*(epetraModel->get_x_map()),x_computed_thyra_ptr);
00288     const Epetra_Vector &x_computed = *x_computed_ptr;
00289 
00290     // compute exact answer
00291     Teuchos::RefCountPtr<Epetra_Vector> x_star_ptr = epetraModel->get_exact_solution(t1);
00292     Epetra_Vector& x_star = *x_star_ptr;
00293 
00294     int MyPID = x_computed.Comm().MyPID();
00295     if (MyPID == 0)
00296     {
00297       *out << "Integrating 1DfemTransient t = " << t0 
00298                 << " to t = " << t1 << std::endl;
00299       *out << "using " << method << "." << std::endl;
00300       *out << "Took " << numSteps << " steps." << std::endl;
00301     }
00302     int MyLength = x_computed.MyLength();
00303     double error = 0;
00304     double errorMag = 0;
00305     for (int i=0 ; i<MyLength ; ++i)
00306     {
00307       if(verbose)
00308       {
00309         *out << std::setprecision(15);
00310         *out << "Computed: x[" << MyPID*MyLength+i << "] = ";
00311         *out << std::setw(20); *out << x_computed[i] << "\t";
00312         *out << "Exact: x[" << MyPID*MyLength+i << "] = ";
00313         *out << std::setw(20); *out << x_star[i] << std::endl;
00314       }
00315       //const double thisError = Thyra::relErr(x_computed[i],x_star[i]);
00316       const double thisError = x_computed[i]-x_star[i];
00317       error += thisError*thisError;
00318       errorMag += x_star[i]*x_star[i];
00319     }
00320     error = sqrt(error)/sqrt(errorMag);
00321     result = Thyra::testMaxErr(
00322       "error",error
00323       ,"maxError",maxError
00324       ,"maxWarning",10.0*maxError
00325       ,&std::cerr,""
00326       );
00327     if(!result) success = false;
00328 
00329 #ifdef HAVE_RYTHMOS_STRATIMIKOS
00330     // Write the final parameters to file
00331     if(W_factory.get())
00332       lowsfCreator.writeParamsFile(*W_factory);
00333 #endif
00334     
00335 #ifdef HAVE_MPI
00336     MPI_Finalize();
00337 #endif // HAVE_MPI
00338 
00339    } // end try
00340    TEUCHOS_STANDARD_CATCH_STATEMENTS(true,*out,success)
00341 
00342   return success ? 0 : 1;
00343 } // end main() [Doxygen looks for this!]
00344 

Generated on Thu Sep 18 12:30:05 2008 for Rythmos - Transient Integration for Differential Equations by doxygen 1.3.9.1