cxx_main.cpp

00001 //@HEADER
00002 
00003 // ***********************************************************************
00004 //
00005 //                     Rythmos Package
00006 //                 Copyright (2006) Sandia Corporation
00007 //
00008 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00009 // license for use of this work by or on behalf of the U.S. Government.
00010 //
00011 // This library is free software; you can redistribute it and/or modify
00012 // it under the terms of the GNU Lesser General Public License as
00013 // published by the Free Software Foundation; either version 2.1 of the
00014 // License, or (at your option) any later version.
00015 //
00016 // This library is distributed in the hope that it will be useful, but
00017 // WITHOUT ANY WARRANTY; without even the implied warranty of
00018 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00019 // Lesser General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License along with this library; if not, write to the Free Software
00023 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00024 // USA
00025 // Questions? Contact Todd S. Coffey (tscoffe@sandia.gov)
00026 //
00027 // ***********************************************************************
00028 //@HEADER 
00029 
00030 #include "Epetra_Map.h"
00031 #include "Epetra_Vector.h"
00032 #include "Epetra_Version.h"
00033 #ifdef HAVE_MPI
00034 #include "Epetra_MpiComm.h"
00035 #include "mpi.h"
00036 #else
00037 #include "Epetra_SerialComm.h"
00038 #endif // HAVE_MPI
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 // 10/9/06 tscoffe:  IntegratorDefault includes: 
00050 #include "Rythmos_InterpolationBuffer.hpp"
00051 #include "Rythmos_LinearInterpolator.hpp"
00052 #include "Rythmos_HermiteInterpolator.hpp"
00053 #include "Rythmos_IntegratorDefault.hpp"
00054 
00055 // Includes for Thyra:
00056 #include "Thyra_EpetraThyraWrappers.hpp"
00057 #include "Thyra_EpetraLinearOp.hpp"
00058 #include "Thyra_EpetraModelEvaluator.hpp"
00059 #include "Thyra_LinearNonlinearSolver.hpp"
00060 #include "Rythmos_TimeStepNonlinearSolver.hpp"
00061 #include "Thyra_TestingTools.hpp"
00062 
00063 // Includes for Stratimikos:
00064 #ifdef HAVE_RYTHMOS_STRATIMIKOS
00065 #  include "Thyra_DefaultRealLinearSolverBuilder.hpp"
00066 #endif
00067 
00068 #include <string>
00069 
00070 // Includes for Teuchos:
00071 #include "Teuchos_as.hpp"
00072 #include "Teuchos_RCP.hpp"
00073 #include "Teuchos_ParameterList.hpp"
00074 #include "Teuchos_CommandLineProcessor.hpp"
00075 #include "Teuchos_VerbosityLevelCommandLineProcessorHelpers.hpp"
00076 #include "Teuchos_FancyOStream.hpp"
00077 #include "Teuchos_GlobalMPISession.hpp"
00078 #include "Teuchos_VerboseObject.hpp"
00079 #include "Teuchos_StandardParameterEntryValidators.hpp"
00080 #include "Teuchos_StandardCatchMacros.hpp"
00081 
00082 enum EMethod { METHOD_FE, METHOD_BE, METHOD_ERK, METHOD_BDF };
00083 enum EStepMethod { STEP_TYPE_FIXED, STEP_TYPE_VARIABLE };
00084 
00085 int main(int argc, char *argv[])
00086 {
00087 
00088   using Teuchos::as;
00089   using Teuchos::RCP;
00090   using Teuchos::rcp;
00091   using Teuchos::Array;
00092 
00093   bool verbose = false; // verbosity level.
00094   bool result, success = true; // determine if the run was successfull
00095 
00096   Teuchos::GlobalMPISession mpiSession(&argc,&argv);
00097 
00098   RCP<Teuchos::FancyOStream>
00099     out = Teuchos::VerboseObjectBase::getDefaultOStream();
00100 
00101 #ifdef HAVE_MPI
00102   MPI_Comm mpiComm = MPI_COMM_WORLD;
00103 #endif // HAVE_MPI
00104 
00105   try { // catch exceptions
00106 
00107     double lambda_min = -0.9;   // min ODE coefficient
00108     double lambda_max = -0.01;  // max ODE coefficient
00109     double coeff_s = 0.0;  // Default is no forcing term
00110     std::string lambda_fit = "linear"; // Lambda model
00111     int numElements = 1; // number of elements in vector
00112     double x0 = 10.0; // ODE initial condition
00113     double finalTime = 1.0; // ODE final time
00114     int N = 10;  // number of steps to take
00115     const int num_methods = 4;
00116     const EMethod method_values[] = { METHOD_FE, METHOD_BE, METHOD_ERK, METHOD_BDF };
00117     const char * method_names[] = { "FE", "BE", "ERK", "BDF" };
00118     EMethod method_val = METHOD_ERK;
00119     const int num_step_methods = 2;
00120     const EStepMethod step_method_values[] = { STEP_TYPE_FIXED, STEP_TYPE_VARIABLE };
00121     const char * step_method_names[] = { "fixed", "variable" };
00122     EStepMethod step_method_val = STEP_TYPE_FIXED;
00123     double maxError = 1e-6;
00124     double maxRestepError = 1.0e4*Teuchos::ScalarTraits<double>::prec();
00125     bool version = false;  // display version information 
00126     double reltol = 1.0e-2;
00127     double abstol = 1.0e-4;
00128     int maxOrder = 5;
00129     bool useIntegrator = false;
00130     int buffersize = 100;
00131     Teuchos::EVerbosityLevel verbLevel = Teuchos::VERB_LOW;
00132 
00133     // Parse the command-line options:
00134     Teuchos::CommandLineProcessor  clp(false); // Don't throw exceptions
00135     clp.addOutputSetupOptions(true);
00136 #ifdef HAVE_RYTHMOS_STRATIMIKOS
00137     Thyra::DefaultRealLinearSolverBuilder lowsfCreator;
00138     lowsfCreator.setupCLP(&clp);
00139 #endif
00140     clp.setOption( "x0", &x0, "Constant ODE initial condition." );
00141     clp.setOption( "T", &finalTime, "Final time for simulation." );
00142     clp.setOption( "lambda_min", &lambda_min, "Lower bound for ODE coefficient");
00143     clp.setOption( "lambda_max", &lambda_max, "Upper bound for ODE coefficient");
00144     clp.setOption( "lambda_fit", &lambda_fit, "Lambda model:  random, linear");
00145     clp.setOption( "force_coeff", &coeff_s, "Forcing term coefficient");
00146     clp.setOption( "numelements", &numElements, "Problem size");
00147     clp.setOption( "method", &method_val, num_methods, method_values, method_names, "Integration method" );
00148     clp.setOption( "stepmethod", &step_method_val, num_step_methods, step_method_values, step_method_names, "Stepping method" );
00149     clp.setOption( "numsteps", &N, "Number of integration steps to take" );
00150     clp.setOption( "maxerror", &maxError, "Maximum error" );
00151     clp.setOption( "max-restep-error", &maxRestepError, "Maximum error between solutions of two steppers" );
00152     clp.setOption( "verbose", "quiet", &verbose, "Set if output is printed or not" );
00153     clp.setOption( "version", "run", &version, "Version of this code" );
00154     clp.setOption( "reltol", &reltol, "Relative Error Tolerance" );
00155     clp.setOption( "abstol", &abstol, "Absolute Error Tolerance" );
00156     clp.setOption( "maxorder", &maxOrder, "Maximum Implicit BDF order" );
00157     clp.setOption( "useintegrator", "normal", &useIntegrator, "Use IntegratorDefault as integrator" );
00158     clp.setOption( "buffersize", &buffersize, "Number of solutions to store in InterpolationBuffer" );
00159     setVerbosityLevelOption( "verb-level", &verbLevel, "Overall verbosity level.", &clp );
00160 
00161     Teuchos::CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse(argc,argv);
00162     if( parse_return != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL ) return parse_return;
00163 
00164     // RAB: 2007/05/14: ToDo: In all of the below code that is called change
00165     // from the "outputLevel" interger parameter to the "Verbose Object"
00166     // sublist with its "Verbosity Level" string parameter.
00167 
00168 #ifdef HAVE_RYTHMOS_STRATIMIKOS
00169     lowsfCreator.readParameters(out.get());
00170     *out << "\nThe parameter list after being read in:\n";
00171     lowsfCreator.getParameterList()->print(*out,2,true,false);
00172 #endif
00173 
00174     if (version) // Display version information and exit.
00175     {
00176       *out << Rythmos::Rythmos_Version() << std::endl; 
00177       *out << "basicExample Version 0.1 - 06/23/05" << std::endl;
00178       return(0);
00179     }
00180 
00181     if (lambda_min > lambda_max)
00182     {
00183       *out << "lamba_min must be less than lambda_max" << std::endl;
00184       return(1);
00185     }
00186 
00187     if (finalTime <= 0.0)
00188     {
00189       *out << "Final simulation time must be > 0.0." << std::endl;
00190       return(1);
00191     }
00192 
00193     *out << std::setprecision(15);
00194     
00195     // Set up the parameter list for the application:
00196     Teuchos::ParameterList params;
00197     bool implicitFlag = ((method_val==METHOD_BE) | (method_val==METHOD_BDF));
00198     //*out << "implicitFlag = " << implicitFlag << std::endl;
00199     params.set( "Implicit", implicitFlag );
00200     params.set( "Lambda_min", lambda_min );
00201     params.set( "Lambda_max", lambda_max );
00202     params.set( "Lambda_fit", lambda_fit );
00203     params.set( "NumElements", numElements );
00204     params.set( "x0", x0 );
00205     params.set( "Coeff_s", coeff_s );
00206 #ifdef HAVE_MPI
00207     RCP<Epetra_Comm> epetra_comm_ptr_ = rcp( new Epetra_MpiComm(mpiComm) );
00208 #else
00209     RCP<Epetra_Comm> epetra_comm_ptr_ = rcp( new Epetra_SerialComm  );
00210 #endif // HAVE_MPI
00211 
00212     // Create the factory for the LinearOpWithSolveBase object
00213     RCP<Thyra::LinearOpWithSolveFactoryBase<double> >
00214       W_factory;
00215     if((method_val == METHOD_BE) | (method_val == METHOD_BDF)) {
00216 #ifdef HAVE_RYTHMOS_STRATIMIKOS
00217       W_factory = lowsfCreator.createLinearSolveStrategy("");
00218       *out
00219         << "\nCreated a LinearOpWithSolveFactory described as:\n"
00220         << Teuchos::describe(*W_factory,Teuchos::VERB_MEDIUM);
00221 #endif
00222     }
00223 
00224     // create interface to problem
00225     RCP<ExampleApplication>
00226       epetraModel = rcp(new ExampleApplication(epetra_comm_ptr_, params));
00227     RCP<Thyra::ModelEvaluator<double> >
00228       model = rcp(new Thyra::EpetraModelEvaluator(epetraModel,W_factory));
00229     RCP<ExampleApplication>
00230       epetraModelSlave = rcp(new ExampleApplication(epetra_comm_ptr_, params));
00231     RCP<Thyra::ModelEvaluator<double> >
00232       modelSlave = rcp(new Thyra::EpetraModelEvaluator(epetraModelSlave,W_factory));
00233 
00234     // Create Stepper object depending on command-line input
00235     std::string method;
00236     RCP<Rythmos::StepperBase<double> > stepper_ptr;
00237     RCP<Rythmos::StepperBase<double> > stepperSlave_ptr;
00238     if ( method_val == METHOD_ERK ) {
00239       stepper_ptr = rcp(new Rythmos::ExplicitRKStepper<double>(model));
00240       RCP<Teuchos::ParameterList> ERKparams = rcp(new Teuchos::ParameterList);
00241       ERKparams->set( "outputLevel", as<int>(verbLevel) );
00242       stepper_ptr->setParameterList(ERKparams);
00243       method = "Explicit Runge-Kutta of order 4";
00244       step_method_val = STEP_TYPE_FIXED;
00245     }
00246     else if (method_val == METHOD_FE) {
00247       stepper_ptr = rcp(new Rythmos::ForwardEulerStepper<double>(model));
00248       RCP<Teuchos::ParameterList> FEparams = rcp(new Teuchos::ParameterList);
00249       FEparams->set( "outputLevel", as<int>(verbLevel));
00250       stepper_ptr->setParameterList(FEparams);
00251       method = "Forward Euler";
00252       step_method_val = STEP_TYPE_FIXED;
00253     }
00254     else if ((method_val == METHOD_BE) | (method_val == METHOD_BDF)) {
00255       RCP<Thyra::NonlinearSolverBase<double> >
00256         nonlinearSolver;
00257       RCP<Thyra::NonlinearSolverBase<double> >
00258         nonlinearSolverSlave;
00259       RCP<Rythmos::TimeStepNonlinearSolver<double> >
00260         _nonlinearSolver = rcp(new Rythmos::TimeStepNonlinearSolver<double>());
00261       RCP<Rythmos::TimeStepNonlinearSolver<double> >
00262         _nonlinearSolverSlave = rcp(new Rythmos::TimeStepNonlinearSolver<double>());
00263       RCP<Teuchos::ParameterList>
00264         nonlinearSolverPL = Teuchos::parameterList();
00265       nonlinearSolverPL->set("Default Tol",double(1e-3*maxError));
00266       _nonlinearSolver->setParameterList(nonlinearSolverPL);
00267       _nonlinearSolverSlave->setParameterList(nonlinearSolverPL);
00268       nonlinearSolver = _nonlinearSolver;
00269       nonlinearSolverSlave = _nonlinearSolverSlave;
00270       if (method_val == METHOD_BE) {
00271         stepper_ptr = rcp(
00272           new Rythmos::BackwardEulerStepper<double>(model,nonlinearSolver));
00273         RCP<Teuchos::ParameterList>
00274           BEparams = rcp(new Teuchos::ParameterList);
00275         BEparams->sublist("VerboseObject").set(
00276           "Verbosity Level",
00277           Teuchos::getVerbosityLevelParameterValueName(verbLevel)
00278           );
00279         stepper_ptr->setParameterList(BEparams);
00280         method = "Backward Euler";
00281         step_method_val = STEP_TYPE_FIXED;
00282       } 
00283       else {
00284         RCP<Teuchos::ParameterList>
00285           BDFparams = rcp(new Teuchos::ParameterList);
00286         RCP<Teuchos::ParameterList> BDFStepControlPL =
00287           Teuchos::sublist(BDFparams,RythmosStepControlSettings_name);
00288 
00289         BDFparams->sublist("VerboseObject").set(
00290           "Verbosity Level",
00291           Teuchos::getVerbosityLevelParameterValueName(verbLevel)
00292           );
00293 
00294         BDFStepControlPL->sublist("VerboseObject").set(
00295           "Verbosity Level",
00296           Teuchos::getVerbosityLevelParameterValueName(verbLevel)
00297           );
00298         BDFStepControlPL->set( "stopTime", finalTime );
00299         BDFStepControlPL->set( "maxOrder", maxOrder );
00300         BDFStepControlPL->set( "relErrTol", reltol );
00301         BDFStepControlPL->set( "absErrTol", abstol );
00302         stepper_ptr = rcp(
00303           new Rythmos::ImplicitBDFStepper<double>(model,nonlinearSolver,BDFparams));
00304         stepperSlave_ptr = rcp(
00305           new Rythmos::ImplicitBDFStepper<double>(modelSlave,nonlinearSolverSlave,BDFparams));
00306         method = "Implicit BDF";
00307         // step_method_val setting is left alone in this case
00308       }
00309     }
00310     else {
00311       TEST_FOR_EXCEPT(true);
00312     }
00313     Rythmos::StepperBase<double> &stepper = *stepper_ptr;
00314     if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) )
00315       stepper.describe(*out,verbLevel);
00316 
00317     int numSteps = 0;
00318     double t0 = 0.0;
00319     double dt = (finalTime-t0)/N;
00320     double time = t0;
00321 
00322     RCP<const Thyra::VectorBase<double> > x_computed_thyra_ptr;
00323     if (step_method_val == STEP_TYPE_FIXED)
00324     {
00325       if (useIntegrator)
00326       {
00327         // Set up fixed-step-size integration:
00328         RCP<Teuchos::ParameterList> 
00329           integratorParams = rcp(new Teuchos::ParameterList);
00330         integratorParams->set( "Take Variable Steps", false );
00331         integratorParams->set( "fixed_dt", dt );
00332         integratorParams->sublist("VerboseObject").set(
00333           "Verbosity Level",
00334           Teuchos::getVerbosityLevelParameterValueName(verbLevel)
00335           );
00336         // Create integrator using stepper and linear interpolation buffer:
00337         RCP<Rythmos::InterpolatorBase<double> > 
00338           linearInterpolator = rcp(new Rythmos::LinearInterpolator<double>());
00339         RCP<Rythmos::InterpolationBuffer<double> > 
00340           IB = rcp(new Rythmos::InterpolationBuffer<double>(linearInterpolator,buffersize));
00341         IB->setParameterList(integratorParams);
00342         Rythmos::IntegratorDefault<double> integrator(stepper_ptr,IB,integratorParams);
00343         // Ask for desired time value:
00344         Array<double> time_vals;
00345         for (int i=0 ; i<=N ; ++i)
00346           time_vals.push_back(i*dt);
00347         Array<RCP<const Thyra::VectorBase<double> > > x_vec;
00348         Array<RCP<const Thyra::VectorBase<double> > > xdot_vec;
00349         Array<double> accuracy_vec;
00350         integrator.getPoints(time_vals,&x_vec,&xdot_vec,&accuracy_vec);
00351         // Get solution out of stepper:
00352         x_computed_thyra_ptr = x_vec.back();
00353       }
00354       else
00355       {
00356         // Integrate forward with fixed step sizes:
00357         for (int i=1 ; i<=N ; ++i)
00358         {
00359           double dt_taken = stepper.takeStep(dt,Rythmos::STEP_TYPE_FIXED);
00360           time += dt_taken;
00361           numSteps++;
00362           if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) )
00363             stepper.describe(*out,verbLevel);
00364           if (dt_taken != dt)
00365           {
00366             cerr << "Error, stepper took step of dt = " << dt_taken 
00367               << " when asked to take step of dt = " << dt << std::endl;
00368             break;
00369           }
00370         }
00371         // Get solution out of stepper:
00372         Rythmos::StepStatus<double> stepStatus = stepper.getStepStatus();
00373         x_computed_thyra_ptr = stepStatus.solution;
00374       }
00375     }
00376     else // (step_method_val == STEP_TYPE_VARIABLE)
00377     {
00378       if (useIntegrator)
00379       {
00380         // Set up fixed-step-size integration:
00381         RCP<Teuchos::ParameterList> 
00382           integratorParams = rcp(new Teuchos::ParameterList);
00383         integratorParams->set( "Take Variable Steps", true );
00384         integratorParams->sublist("VerboseObject").set(
00385           "Verbosity Level",
00386           Teuchos::getVerbosityLevelParameterValueName(verbLevel)
00387           );
00388         // Create integrator using stepper and interpolation buffer:
00389         //RCP<Rythmos::InterpolatorBase<double> > 
00390         //  linearInterpolator = rcp(new Rythmos::LinearInterpolator<double>());
00391         //RCP<Rythmos::InterpolationBuffer<double> > 
00392         //  IB = rcp(new Rythmos::InterpolationBuffer<double>(linearInterpolator,buffersize));
00393         RCP<Rythmos::InterpolatorBase<double> > 
00394           hermiteInterpolator = rcp(new Rythmos::HermiteInterpolator<double>());
00395         RCP<Rythmos::InterpolationBuffer<double> > 
00396           IB = rcp(new Rythmos::InterpolationBuffer<double>(hermiteInterpolator,buffersize));
00397         IB->setParameterList(integratorParams);
00398         Rythmos::IntegratorDefault<double> integrator(stepper_ptr,IB,integratorParams);
00399         // Ask for desired time value:
00400         Array<double> time_vals;
00401         for (int i=0 ; i<=N ; ++i)
00402           time_vals.push_back(i*dt);
00403         Array<RCP<const Thyra::VectorBase<double> > > x_vec;
00404         Array<RCP<const Thyra::VectorBase<double> > > xdot_vec;
00405         Array<double> accuracy_vec;
00406         integrator.getPoints(time_vals,&x_vec,&xdot_vec,&accuracy_vec);
00407         // Get solution out of stepper:
00408         x_computed_thyra_ptr = x_vec.back();
00409       }
00410       else
00411       {
00412         // HIGH output data structures:
00413         // Create a place to store the computed solutions
00414         Rythmos::StepStatus<double> stepStatus = stepper.getStepStatus();
00415         x_computed_thyra_ptr = stepStatus.solution;
00416         // Convert Thyra::VectorBase to Epetra_Vector
00417         // Create a place to store the exact numerical solution
00418         RCP<Epetra_Vector> x_numerical_exact_ptr = rcp(new Epetra_Vector(*epetraModel->get_x_map()));
00419         Epetra_Vector& x_numerical_exact = *x_numerical_exact_ptr;
00420         // Create a place to store the relative difference:
00421         RCP<Epetra_Vector> x_rel_diff_ptr = rcp(new Epetra_Vector(*epetraModel->get_x_map()));
00422         Epetra_Vector& x_rel_diff = *x_rel_diff_ptr;
00423         // get lambda from the problem:
00424         RCP<const Epetra_Vector> lambda_ptr = epetraModel->get_coeff();
00425         const Epetra_Vector &lambda = *lambda_ptr;
00426 
00427         while (time < finalTime)
00428         {
00429           double dt_taken = stepper.takeStep(0.0,Rythmos::STEP_TYPE_VARIABLE);
00430           if (method_val == METHOD_BDF) {
00431             stepperSlave_ptr->setStepControlData(stepper);
00432             double slave_dt_taken = stepperSlave_ptr->takeStep(dt_taken,Rythmos::STEP_TYPE_FIXED);
00433             // Check that returned dt matches exactly
00434             TEST_FOR_EXCEPT(dt_taken != slave_dt_taken);
00435             Rythmos::StepStatus<double> stepStatusMaster = stepper.getStepStatus();
00436             Rythmos::StepStatus<double> stepStatusSlave = stepperSlave_ptr->getStepStatus();
00437             // Check that the stepStatus matches exactly:
00438             TEST_FOR_EXCEPT(stepStatusMaster.stepLETStatus != stepStatusSlave.stepLETStatus);
00439             TEST_FOR_EXCEPT(stepStatusMaster.stepSize != stepStatusSlave.stepSize);
00440             TEST_FOR_EXCEPT(stepStatusMaster.time != stepStatusSlave.time);
00441             if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00442               *out << "Master order = " << stepStatusMaster.order << endl;
00443               *out << " Slave order = " << stepStatusSlave.order << endl;
00444               *out << "Master LET Value = " << stepStatusMaster.stepLETValue << endl;
00445               *out << " Slave LET Value = " << stepStatusSlave.stepLETValue << endl;
00446             }
00447             TEST_FOR_EXCEPT(stepStatusMaster.order != stepStatusSlave.order);
00448             // We will allow a difference of some multiplier of machine epsilon:
00449             const double
00450               normLETDiff = std::abs(stepStatusMaster.stepLETValue - stepStatusSlave.stepLETValue);
00451             TEST_FOR_EXCEPTION(
00452               normLETDiff > maxRestepError, std::logic_error,
00453               "Error, normLETDiff = " << normLETDiff << " > maxRestepError = " << maxRestepError << "!" );
00454             if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00455               *out << "normLETDiff = " << normLETDiff << std::endl;
00456             }
00457             // Create a non-const Thyra VectorBase to use as a temp vector
00458             RCP<Thyra::VectorBase<double> > vec_temp = stepStatusSlave.solution->clone_v();
00459             // Check that the solution matches exactly
00460             Thyra::V_StVpStV<double>(&*vec_temp,1.0,*stepStatusMaster.solution,-1.0,*stepStatusSlave.solution);
00461             double normSolutionDiff = Thyra::norm_inf<double>(*vec_temp);
00462             if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00463               *out << "normSolutionDiff = " << normSolutionDiff << endl;
00464             }
00465             const double eps = 1.0e4*Teuchos::ScalarTraits<double>::prec();
00466             TEST_FOR_EXCEPT(normSolutionDiff > eps);
00467             // Check that solution dot matches exactly
00468             {
00469 
00470               RCP<const Thyra::VectorBase<double> >
00471                 master_x_dot = get_xdot(stepper,stepStatusMaster.time),
00472                 slave_x_dot = get_xdot(*stepperSlave_ptr,stepStatusSlave.time);
00473               Thyra::V_VmV<double>(&*vec_temp,*master_x_dot,*slave_x_dot);
00474               double normSolutionDotDiff = Thyra::norm_inf<double>(*vec_temp);
00475               if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00476                 *out << "normSolutionDotDiff = " << normSolutionDotDiff << endl;
00477               }
00478               TEST_FOR_EXCEPT(normSolutionDotDiff > eps);
00479             }
00480             // Do not check that the residual matches because the residual isn't stored in ImplicitBDFStepper.
00481           }
00482           numSteps++;
00483           if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) )
00484             stepper.describe(*out,verbLevel);
00485           if (dt_taken < 0)
00486           {
00487             *out << "Error, stepper failed for some reason with step taken = " << dt_taken << endl;
00488             break;
00489           }
00490           if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) )
00491           {
00492             // Get solution out of stepper:
00493             stepStatus = stepper.getStepStatus();
00494             x_computed_thyra_ptr = stepStatus.solution;
00495             // Convert Thyra::VectorBase to Epetra_Vector
00496             RCP<const Epetra_Vector>
00497               x_computed_ptr = Thyra::get_Epetra_Vector(*(epetraModel->get_x_map()),x_computed_thyra_ptr);
00498             if ((method_val == METHOD_BDF) && (maxOrder == 1))
00499             {
00500               int myN = x_numerical_exact.MyLength();
00501               if (numSteps == 1) // First step
00502               {
00503                 for (int i=0 ; i<myN ; ++i)
00504                   x_numerical_exact[i] = x0;
00505               }
00506               for (int i=0 ; i<myN ; ++i)
00507                 x_numerical_exact[i] = ( x_numerical_exact[i]
00508                     +dt_taken*epetraModel->evalR(time+dt_taken,lambda[i],coeff_s))
00509                     /(1-lambda[i]*dt_taken);
00510               for (int i=0 ; i<myN ; ++i)
00511                 x_rel_diff[i] = (x_numerical_exact[i]-(*x_computed_ptr)[i])/x_numerical_exact[i];
00512               if (myN == 1)
00513                 *out << "Computed x(" << time+dt_taken << ") = " << (*x_computed_ptr)[0] 
00514                     << "  Numerical Exact = " << x_numerical_exact[0] 
00515                     << "  Rel Diff = " << x_rel_diff[0] << std::endl;
00516               else
00517               {
00518                 for (int i=0 ; i<myN ; ++i)
00519                   *out << "Computed x_" << i << "(" << time+dt_taken << ") = " << (*x_computed_ptr)[i] 
00520                       << "  Numerical Exact = " << x_numerical_exact[i] 
00521                       << "  Rel Diff = " << x_rel_diff[i] <<  std::endl;
00522               }
00523             }
00524             else
00525             {
00526               // compute exact answer
00527               RCP<const Epetra_Vector>
00528                 x_star_ptr = epetraModel->getExactSolution(time);
00529               const Epetra_Vector& x_star = *x_star_ptr;
00530               int myN = x_computed_ptr->MyLength();
00531               for (int i=0 ; i<myN ; ++i)
00532                 x_rel_diff[i] = (x_star[i]-(*x_computed_ptr)[i])/x_star[i];
00533               if (myN == 1)
00534                 *out << "Computed x(" << time+dt_taken << ") = " << (*x_computed_ptr)[0] 
00535                     << "  Exact = " << x_star[0] 
00536                     << "  Rel Diff = " << x_rel_diff[0] << std::endl;
00537               else
00538               {
00539                 for (int i=0 ; i<myN ; ++i)
00540                   *out << "Computed x_" << i << "(" << time+dt_taken << ") = " << (*x_computed_ptr)[i] 
00541                       << "  Exact = " << x_star[i] 
00542                       << "  Rel Diff = " << x_rel_diff[i] << std::endl;
00543               }
00544             }
00545           }
00546           time += dt_taken;
00547           *out << "Took stepsize of: " << dt_taken << " time = " << time << endl;
00548         }
00549         // Get solution out of stepper:
00550         stepStatus = stepper.getStepStatus();
00551         x_computed_thyra_ptr = stepStatus.solution;
00552       }
00553     }
00554     *out << "Integrated to time = " << time << endl;
00555 
00556     // Convert solution from Thyra::VectorBase to Epetra_Vector
00557     RCP<const Epetra_Vector>
00558       x_computed_ptr = Thyra::get_Epetra_Vector(*(epetraModel->get_x_map()),x_computed_thyra_ptr);
00559     const Epetra_Vector &x_computed = *x_computed_ptr;
00560 
00561     // compute exact answer
00562     RCP<const Epetra_Vector>
00563       x_star_ptr = epetraModel->getExactSolution(finalTime);
00564     const Epetra_Vector& x_star = *x_star_ptr;
00565     
00566     // get lambda from the problem:
00567     RCP<const Epetra_Vector> lambda_ptr = epetraModel->get_coeff();
00568     const Epetra_Vector &lambda = *lambda_ptr;
00569 
00570     // compute numerical exact answer (for FE and BE)
00571     RCP<const Epetra_Vector> x_numerical_exact_ptr; 
00572     if (method_val == METHOD_FE) 
00573     {
00574       RCP<Epetra_Vector> x_exact_ptr = rcp(new Epetra_Vector(x_star.Map()));
00575       Epetra_Vector& x_exact = *x_exact_ptr;
00576       int myN = x_exact.MyLength();
00577       for ( int i=0 ; i<myN ; ++i)
00578       {
00579         x_exact[i] = x0;
00580         for (int j=1 ; j<=N ; ++j)
00581         {
00582           x_exact[i] = (1+lambda[i]*dt)*x_exact[i]+dt*epetraModel->evalR(0+j*dt,lambda[i],coeff_s);
00583         }
00584         //x_exact[i] = x0*pow(1+lambda[i]*dt,N);
00585       }
00586       x_numerical_exact_ptr = x_exact_ptr;
00587     } 
00588     else if (method_val == METHOD_BE) 
00589     {
00590       RCP<Epetra_Vector> x_exact_ptr = rcp(new Epetra_Vector(x_star.Map()));
00591       Epetra_Vector& x_exact = *x_exact_ptr;
00592       int myN = x_exact.MyLength();
00593       for ( int i=0 ; i<myN ; ++i)
00594       {
00595         x_exact[i] = x0;
00596         for (int j=1 ; j<=N ; ++j)
00597         {
00598           x_exact[i] = (x_exact[i]+dt*epetraModel->evalR(0+j*dt,lambda[i],coeff_s))/(1-lambda[i]*dt);
00599         }
00600         //x_exact[i] = x0*pow(1/(1-lambda[i]*dt),N);
00601       }
00602       x_numerical_exact_ptr = x_exact_ptr;
00603     }
00604     else if (method_val == METHOD_BDF)
00605     {
00606       // exact bdf solution here?
00607     }
00608 
00609     // 06/03/05 tscoffe to get an Epetra_Map associated with an Epetra_Vector:
00610     // x.Map()
00611     // to get an Epetra_Comm associated with an Epetra_Vector:
00612     // x.Comm()
00613     
00614     int MyPID = x_computed.Comm().MyPID();
00615     if (MyPID == 0)
00616     {
00617       *out << "Integrating \\dot{x}=\\lambda x from t = " << t0 
00618                 << " to t = " << finalTime << std::endl;
00619       *out << "using " << method << std::endl;
00620       *out << "with initial x_0 = " << x0
00621                 << ", \\Delta t = " << dt  << "." << std::endl;
00622       *out << "Took " << numSteps << " steps." << std::endl;
00623     }
00624     int MyLength = x_computed.MyLength();
00625     if (verbose)
00626     {
00627       for (int i=0 ; i<MyLength ; ++i)
00628       {
00629         *out << std::setprecision(15);
00630         *out << "lambda[" << MyPID*MyLength+i << "] = " << lambda[i] << std::endl;
00631       }
00632       // Print out computed and exact solutions:
00633       for (int i=0 ; i<MyLength ; ++i)
00634       {
00635         *out << std::setprecision(15);
00636         *out << "Computed: x[" << MyPID*MyLength+i << "] = ";
00637         *out << std::setw(20) << x_computed[i] << "\t";
00638         *out << "Exact: x[" << MyPID*MyLength+i << "] = ";
00639         *out << std::setw(20) << x_star[i] << std::endl;
00640       }
00641     }
00642     
00643     // Check numerics against exact numerical method for FE and BE case:
00644     double numerical_error = 0;
00645     double numerical_error_mag = 0;
00646     if (x_numerical_exact_ptr.get())
00647     {
00648       const Epetra_Vector& x_numerical_exact = *x_numerical_exact_ptr;
00649       for ( int i=0 ; i<MyLength ; ++i)
00650       {
00651         if (verbose) 
00652         {
00653           *out << std::setprecision(15);
00654           *out << "Computed: x[" << MyPID*MyLength+i << "] = ";
00655           *out << std::setw(20) << x_computed[i] << "\t";
00656           *out << "Numerical Exact: x[" << MyPID*MyLength+i << "] = ";
00657           *out << std::setw(20) << x_numerical_exact[i] << std::endl;
00658         }
00659         const double thisError = x_numerical_exact[i]-x_computed[i];
00660         numerical_error += thisError*thisError;
00661         numerical_error_mag += x_numerical_exact[i]*x_numerical_exact[i];
00662       }
00663       numerical_error = sqrt(numerical_error)/sqrt(numerical_error_mag);
00664       result = Thyra::testMaxErr(
00665         "Exact numerical error",numerical_error
00666         ,"maxError",1.0e-10
00667         ,"maxWarning",1.0e-9
00668         ,&*out,""
00669         );
00670       if(!result) success = false;
00671     }
00672 
00673     // Check numerics against exact DE solution:
00674     double error = 0;
00675     double errorMag = 0;
00676     for (int i=0 ; i<MyLength ; ++i)
00677     {
00678       const double thisError = x_computed[i]-x_star[i];
00679       error += thisError*thisError;
00680       errorMag += x_star[i]*x_star[i];
00681     }
00682     error = sqrt(error)/sqrt(errorMag);
00683     result = Thyra::testMaxErr(
00684       "Exact DE solution error",error
00685       ,"maxError",maxError
00686       ,"maxWarning",10.0*maxError
00687       ,&*out,""
00688       );
00689     if(!result) success = false;
00690 
00691 #ifdef HAVE_RYTHMOS_STRATIMIKOS
00692     // Write the final parameters to file
00693     if(W_factory.get())
00694       lowsfCreator.writeParamsFile(*W_factory);
00695 #endif
00696 
00697    } // end try
00698   TEUCHOS_STANDARD_CATCH_STATEMENTS(true,*out,success)
00699 
00700   return success ? 0 : 1;
00701 
00702 } // end main() [Doxygen looks for this!]

Generated on Tue Oct 20 12:46:07 2009 for Rythmos - Transient Integration for Differential Equations by doxygen 1.4.7