SimulatorTester.hpp

Go to the documentation of this file.

//===========================================================================
//
// File: SimulatorTester.hpp
//
// Created: Fri Aug  7 09:21:55 2009
//
// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
//            Bård Skaflestad     <bard.skaflestad@sintef.no>
//
// $Date$
//
// $Revision$
//
//===========================================================================

/*
  Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
  Copyright 2009, 2010 Statoil ASA.

  This file is part of The Open Reservoir Simulator Project (OpenRS).

  OpenRS is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation, either version 3 of the License, or
  (at your option) any later version.

  OpenRS is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
*/

#ifndef OPENRS_SIMULATORTESTER_HEADER
#define OPENRS_SIMULATORTESTER_HEADER


#include <opm/porsol/common/SimulatorBase.hpp>
#include <dune/grid/io/file/vtk/vtkwriter.hh>

#include <iostream>

namespace Opm
{



    template <class SimTraits>
    class SimulatorTester : public SimulatorBase<SimTraits>
    {
    public:
    typedef SimulatorBase<SimTraits> Super;
    void run()
    {
        // Initial saturation.
        std::vector<double> saturation(this->init_saturation_);
        std::vector<double> saturation_old(saturation);
        // Gravity.
        // Dune::FieldVector<double, 3> gravity(0.0);
        // gravity[2] = -Dune::unit::gravity;
        // Compute flow field.
        if (this->gravity_.two_norm() > 0.0) {
        OPM_MESSAGE("Warning: Gravity not handled by flow solver.");
        }

        // Solve some steps.
        for (int i = 0; i < this->simulation_steps_; ++i) {
        std::cout << "\n\n================    Simulation step number " << i
                          << "    ===============" << std::endl;
        // Flow.
        this->flow_solver_.solve(this->res_prop_, saturation, this->bcond_, this->injection_rates_psolver_,
                                         this->residual_tolerance_, this->linsolver_verbosity_, this->linsolver_type_);
//      if (i == 0) {
//          flow_solver_.printSystem("linsys_dump_mimetic");
//      }
        // Transport.
        this->transport_solver_.transportSolve(saturation, this->stepsize_, this->gravity_,
                            this->flow_solver_.getSolution(),
                            this->injection_rates_);
        // Output.
                writeVtkOutput(this->ginterf_,
                               this->res_prop_,
                               this->flow_solver_.getSolution(),
                               saturation,
                               "testsolution-" + boost::lexical_cast<std::string>(i));

                writeField(saturation, "saturation-" + boost::lexical_cast<std::string>(i));

                // Comparing old to new.
                int num_cells = saturation.size();
                double maxdiff = 0.0;
                for (int i = 0; i < num_cells; ++i) {
                    maxdiff = std::max(maxdiff, std::fabs(saturation[i] - saturation_old[i]));
                }
                std::cout << "Maximum saturation change: " << maxdiff << std::endl;

                // Copy to old.
                saturation_old = saturation;
        }
    }

    };



} // namespace Opm


#endif // OPENRS_SIMULATORTESTER_HEADER