solver.h

/*

freesteam - IAPWS-IF97 steam tables library
Copyright (C) 2004-2005  John Pye

This program 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 2
of the License, or (at your option) any later version.

This program 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 this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

*/

#ifndef SOLVER_H
#define SOLVER_H

#include "steamcalculator.h"
#include "steamproperty.h"
#include "zeroin.h"
#include "exception.h"
#include <string>
using namespace std;


template<class MainProp, class OtherProp, class VaryProp, int MainPropAlt,int OtherPropAlt,int VaryPropAlt>
class SolverBase : public DesignByContract
{

        friend class ZeroIn<SolverBase,VaryProp,OtherProp>;

        protected:
                SolverBase(const MainProp &mp,const OtherProp &op){
                        //cerr << "Running SolverBase creator" << endl;
                        this->mp=mp;
                        this->op=op;
                }

                virtual ~SolverBase(){}

        public:

                virtual void setLowerBound(const VaryProp &lowerbound){
                        this->lowerbound = lowerbound;
                }
                virtual void setUpperBound(const VaryProp &upperbound){
                        this->upperbound = upperbound;
                }


                void setTarget(const OtherProp &op){
                        this->op=op;
                }

                void setMainProp(const MainProp &mp){
                        this->mp=mp;
                }


                SteamCalculator
                solve(
                        OtherProp   maxerror
                        ,VaryProp   tol
                ){
                        ZeroIn<SolverBase,OtherProp,VaryProp> z;

                        //cerr << endl << "Solving for "<< SteamProperty<OtherProp,OtherPropAlt>::name() << " = " << op << " by varying " << SteamProperty<VaryProp,VaryPropAlt>::name() << ", and with " << SteamProperty<MainProp,MainPropAlt>::name() << " fixed at " << mp << endl;

                        try{

                                //S.set_pT(P_CRIT,T);

                                z.setLowerBound(lowerbound);
                                z.setUpperBound(upperbound);
                                z.setTolerance(tol);

                                z.setMethod(&SolverBase::getError_vp);

                                z.visit(this);

                                if(!z.isSolved(maxerror)){
                                        stringstream s;
                                        s.flags(ios_base::showbase);
                                        s << "Failed solution: target "<< SteamProperty<OtherProp,OtherPropAlt>::name() << " = " << op << ", with " << SteamProperty<MainProp,MainPropAlt>::name() << " fixed at " << mp << endl;

                                        s << " (error was " << z.getError() << ", max allowed is " << maxerror << ")";
                                        throw new Exception(s.str());
                                }

                                setVaryProp(z.getSolution());

                                return S;
                        }catch(Exception *e){
                                stringstream s;
                                s << "Solver<" << SteamProperty<MainProp,MainPropAlt>::name() << "," << SteamProperty<OtherProp,OtherPropAlt>::name() << "," << SteamProperty<VaryProp,VaryPropAlt>::name() << ">::solve: " << e->what();
                                delete e;
                                throw new Exception(s.str());
                        }
                }
        protected:

                virtual void setVaryProp(const VaryProp &vp) = 0;

                OtherProp getError_vp(const VaryProp &vp){
                        //cerr << endl << "Solver::getError_vp: " << SteamProperty<VaryProp,VaryPropAlt>::name() << " = " << vp;

                        setVaryProp(vp);

                        //cerr << " -> " << SteamProperty<OtherProp,OtherPropAlt>::name() << " = " << SteamProperty<OtherProp,OtherPropAlt>::get(S) ;

                        return SteamProperty<OtherProp,OtherPropAlt>::get(S) - op;
                }

        protected:

                OtherProp op;
                MainProp mp;
                SteamCalculator S;
                VaryProp lowerbound;
                VaryProp upperbound;
};


template<class MainProp, class OtherProp, class VaryProp, int MainPropAlt=0,int OtherPropAlt=0,int VaryPropAlt=0>
class Solver
        : public SolverBase<MainProp,OtherProp,VaryProp,MainPropAlt,OtherPropAlt,VaryPropAlt>
{
        public:

                Solver(const MainProp &mp,const OtherProp &op)
                        : SolverBase<MainProp,OtherProp,VaryProp,MainPropAlt,OtherPropAlt,VaryPropAlt>(mp,op){
                        // not implemented
                        throw new Exception("Not implemented");
                }


};

/*

// Solving for specific volume given temperature
class Solver<Temperature,SpecificVolume,Pressure,0,0,0>
        : public SolverBase<Temperature,SpecificVolume,Pressure,0,0,0>
{
        public:

                Solver(Temperature T,SpecificVolume v, bool isSuperheated=false)
                        : SolverBase<Temperature,SpecificVolume,Pressure,0,0,0>(T,v){


                        if(T > T_CRIT){
                                lowerbound = P_MIN;
                                upperbound = P_MAX;
                        }else{
                                if(isSuperheated){
                                        lowerbound = Boundaries::getSatPres_T(T);
                                        upperbound = P_MAX;
                                }else{
                                        lowerbound = P_MIN;
                                        upperbound = Boundaries::getSatPres_T(T);
                                }
                        }
                }

        protected:

                virtual void setVaryProp(const Pressure &p){
                        S.set_pT(p,mp);
                }

};

*/


template<class OtherProp, int OtherPropAlt>
class Solver<Temperature,OtherProp,Pressure,0,OtherPropAlt,0>
        : public SolverBase<Temperature,OtherProp,Pressure,0,OtherPropAlt,0>
{
        public:

                Solver(Temperature T,OtherProp op)
                        : SolverBase<Temperature,OtherProp,Pressure,0,OtherPropAlt,0>(T,op){
                        SolverBase<Temperature,OtherProp,Pressure,0,OtherPropAlt,0>::lowerbound = P_TRIPLE;
                        SolverBase<Temperature,OtherProp,Pressure,0,OtherPropAlt,0>::upperbound = P_MAX;
                }

        protected:

                virtual void setVaryProp(const Pressure &p){
                        SolverBase<Temperature,OtherProp,Pressure,0,OtherPropAlt,0>::S.set_pT(p,SolverBase<Temperature,OtherProp,Pressure,0,OtherPropAlt,0>::mp);
                }
};


template<class OtherProp, int OtherPropAlt>
class Solver<Pressure,OtherProp,Temperature,0,OtherPropAlt,0>
        : public SolverBase<Pressure,OtherProp,Temperature,0,OtherPropAlt,0>
{
        public:

                Solver(Pressure p, OtherProp op)
                        : SolverBase<Pressure,OtherProp,Temperature,0,OtherPropAlt,0>(p,op){
                        SolverBase<Pressure,OtherProp,Temperature,0,OtherPropAlt,0>::lowerbound = T_MIN;
                        SolverBase<Pressure,OtherProp,Temperature,0,OtherPropAlt,0>::upperbound = T_MAX;
                }

        protected:

                virtual void setVaryProp(const Temperature &T){
                        SolverBase<Pressure,OtherProp,Temperature,0,OtherPropAlt,0>::S.set_pT(SolverBase<Pressure,OtherProp,Temperature,0,OtherPropAlt,0>::mp,T);
                }
};

#endif

---