Go to the documentation of this file.00001 
00002 
00003 
00004 
00005 
00006 
00007 
00008 
00009 
00010 
00011 
00012 
00013 
00014 
00015 
00016 
00017 
00018 
00019 
00020 
00021 
00022 
00023 
00024 
00025 
00026 #include "ode.H"
00027 #include <chemistryModel/ODEChemistryModel.H>
00028 
00029 
00030 
00031 template<class CompType, class ThermoType>
00032 Foam::ode<CompType, ThermoType>::ode
00033 (
00034     ODEChemistryModel<CompType, ThermoType>& model,
00035     const word& modelName
00036 )
00037 :
00038     chemistrySolver<CompType, ThermoType>(model, modelName),
00039     coeffsDict_(model.subDict(modelName + "Coeffs")),
00040     solverName_(coeffsDict_.lookup("ODESolver")),
00041     odeSolver_(ODESolver::New(solverName_, model)),
00042     eps_(readScalar(coeffsDict_.lookup("eps"))),
00043     scale_(readScalar(coeffsDict_.lookup("scale")))
00044 {}
00045 
00046 
00047 
00048 
00049 template<class CompType, class ThermoType>
00050 Foam::ode<CompType, ThermoType>::~ode()
00051 {}
00052 
00053 
00054 
00055 
00056 template<class CompType, class ThermoType>
00057 Foam::scalar Foam::ode<CompType, ThermoType>::solve
00058 (
00059     scalarField& c,
00060     const scalar T,
00061     const scalar p,
00062     const scalar t0,
00063     const scalar dt
00064 ) const
00065 {
00066     label nSpecie = this->model_.nSpecie();
00067     scalarField c1(this->model_.nEqns(), 0.0);
00068 
00069     
00070     for (label i=0; i<nSpecie; i++)
00071     {
00072         c1[i] = c[i];
00073     }
00074     c1[nSpecie] = T;
00075     c1[nSpecie+1] = p;
00076 
00077     scalar dtEst = dt;
00078 
00079     odeSolver_->solve
00080     (
00081         this->model_,
00082         t0,
00083         t0 + dt,
00084         c1,
00085         eps_,
00086         dtEst
00087     );
00088 
00089     for (label i=0; i<c.size(); i++)
00090     {
00091         c[i] = max(0.0, c1[i]);
00092     }
00093 
00094     return dtEst;
00095 }
00096 
00097 
00098