From c07bf8ca71ece94c93834ff3f145270df412e3b2 Mon Sep 17 00:00:00 2001 From: michael Date: Fri, 19 Nov 2021 13:28:15 +0100 Subject: [PATCH] Implemented ICoCo v2 --- CoreFlows/Models/inc/DiffusionEquation.hxx | 39 +- CoreFlows/Models/inc/ProblemCoreFlows.hxx | 179 ++- CoreFlows/Models/inc/ProblemFluid.hxx | 9 +- CoreFlows/Models/inc/SinglePhase.hxx | 2 +- .../inc/StationaryDiffusionEquation.hxx | 44 +- CoreFlows/Models/inc/TransportEquation.hxx | 62 +- CoreFlows/Models/src/DiffusionEquation.cxx | 974 ++++++++-------- CoreFlows/Models/src/ProblemCoreFlows.cxx | 99 +- CoreFlows/Models/src/ProblemFluid.cxx | 45 +- .../src/StationaryDiffusionEquation.cxx | 1021 +++++++++-------- CoreFlows/Models/src/TransportEquation.cxx | 638 +++++----- 11 files changed, 1790 insertions(+), 1322 deletions(-) diff --git a/CoreFlows/Models/inc/DiffusionEquation.hxx b/CoreFlows/Models/inc/DiffusionEquation.hxx index 2a3110b..614d749 100755 --- a/CoreFlows/Models/inc/DiffusionEquation.hxx +++ b/CoreFlows/Models/inc/DiffusionEquation.hxx @@ -54,7 +54,7 @@ public : DiffusionEquation( int dim,bool FECalculation=true,double rho=10000,double cp=300,double lambda=5, MPI_Comm comm = MPI_COMM_WORLD); - //Gestion du calcul + //Gestion du calcul (ICoCo) void initialize(); void terminate();//vide la mémoire et enregistre le résultat final bool initTimeStep(double dt); @@ -97,29 +97,32 @@ public : void setConductivity(double conductivite){ _conductivity=conductivite; }; - void setFluidTemperatureField(Field coupledTemperatureField){ - _fluidTemperatureField=coupledTemperatureField; - _fluidTemperatureFieldSet=true; - }; - void setDiffusiontensor(Matrix DiffusionTensor){ _DiffusionTensor=DiffusionTensor; }; - void setFluidTemperature(double fluidTemperature){ - _fluidTemperature=fluidTemperature; - } - //get output fields for postprocessing or coupling - vector getOutputFieldsNames() ;//liste tous les champs que peut fournir le code pour le postraitement - Field& getOutputField(const string& nameField );//Renvoie un champs pour le postraitement + /** Set input fields to prepare the simulation or coupling **/ + vector getInputFieldsNames(); + void setInputField(const string& nameField, Field& inputField );//supply of a required input field - Field& getRodTemperatureField(){ - return _VV; + void setFluidTemperatureField(Field coupledTemperatureField){ + _fluidTemperatureField=coupledTemperatureField; + _fluidTemperatureFieldSet=true; + }; + void setFluidTemperature(double fluidTemperature){ + _fluidTemperature=fluidTemperature; } Field& getFluidTemperatureField(){ return _fluidTemperatureField; } + + /*** get output fields names for postprocessing or coupling ***/ + vector getOutputFieldsNames() ;//liste tous les champs que peut fournir le code pour le postraitement + Field& getOutputField(const string& nameField );//Renvoie un champs pour le postraitement + + Field& getOutputTemperatureField();//Return the main unknown if present (initialize() should be called first) + Field& getRodTemperatureField();//Return the main unknown if present (initialize() should be called first) /*********** Generic functions for finite element method ***********/ static Vector gradientNodal(Matrix M, vector< double > v);//gradient of nodal shape functions @@ -141,7 +144,12 @@ protected : Vector _normale; Matrix _DiffusionTensor; Vec _Tn, _deltaT, _Tk, _Tkm1, _b0; + Vec _Tn_seq; // Local sequential copy of the parallel vector _Tn, used for saving result files + double _dt_diffusion, _dt_src; + TimeScheme _timeScheme; + map _limitField; + /************ Data for FE calculation *************/ int _NunknownNodes;/* number of unknown nodes for FE calculation */ @@ -150,9 +158,6 @@ protected : std::vector< int > _boundaryNodeIds;/* List of boundary nodes */ std::vector< int > _dirichletNodeIds;/* List of boundary nodes with Dirichlet BC */ - TimeScheme _timeScheme; - map _limitField; - /********* Possibility to set a boundary field as DirichletNeumann boundary condition *********/ bool _dirichletValuesSet; bool _neumannValuesSet; diff --git a/CoreFlows/Models/inc/ProblemCoreFlows.hxx b/CoreFlows/Models/inc/ProblemCoreFlows.hxx index ac9e73c..6ec3c60 100755 --- a/CoreFlows/Models/inc/ProblemCoreFlows.hxx +++ b/CoreFlows/Models/inc/ProblemCoreFlows.hxx @@ -170,7 +170,7 @@ public : * \param [in] bool, passage par reférence. * \param [out] bool * */ - virtual bool iterateTimeStep(bool &converged) = 0; //?? + virtual bool iterateTimeStep(bool &converged) = 0; /** \fn isStationary * \brief vérifie la stationnairité du problème . @@ -188,27 +188,121 @@ public : * */ virtual double presentTime() const; + /** \fn setStationaryMode + * \brief Perform the search of a stationary regime + * \details + * \param [in] bool + * \param [out] + * */ + virtual void setStationaryMode(bool stationaryMode){ _stationaryMode=stationaryMode;}; + + /** \fn getStationaryMode + * \brief Tells if we are seeking a stationary regime + * \details + * \param [in] + * \param [out] bool + * */ + virtual bool getStationaryMode(){return _stationaryMode;}; + + /** \fn resetTime + * \brief sets the current time (typically to start a new calculation) + * \details + * \param [in] double + * \param [out] void + * */ + void resetTime (double time); + /* //Coupling interface - virtual vector getInputFieldsNames()=0 ;//Renvoie les noms des champs dont le problème a besoin (données initiales) - virtual Field& getInputFieldTemplate(const string& name)=0;//Renvoie le format de champs attendu (maillage, composantes etc) - virtual void setInputField(const string& name, const Field& afield)=0;//enregistre les valeurs d'une donnée initiale + virtual void getInputMEDDoubleFieldTemplate(const std::string& name, MEDDoubleField& afield) const;//Renvoie le format de champs attendu (maillage, composantes etc) virtual vector getOutputFieldsNames()=0 ;//liste tous les champs que peut fournir le code pour le postraitement - virtual Field& getOutputField(const string& nameField )=0;//Renvoie un champs pour le postraitement + virtual void getOutputMEDDoubleField(const std::string& name, MEDDoubleField& afield)//Renvoie un champs pour le postraitement ou le couplage */ + /** Set input fields to prepare the simulation or coupling **/ + virtual vector getInputFieldsNames()=0; + virtual void setInputField(const string& nameField, Field& inputField )=0;//supply of a required input field + + /*! @brief (Optional) Provide the code with a scalar double data. + * + * See Problem documentation for more details on the time semantic of a scalar value. + * + * @param[in] name name of the scalar value that is given to the code. + * @param[in] val value passed to the code. + * @throws ICoCo::WrongArgument exception if the scalar name ('name' parameter) is invalid. + */ + /* virtual void setInputDoubleValue(const std::string& name, const double& val); */ + + /*! @brief (Optional) Retrieve a scalar double value from the code. + * + * See Problem documentation for more details on the time semantic of a scalar value. + * + * @param[in] name name of the scalar value to be read from the code. + * @return the double value read from the code. + * @throws ICoCo::WrongArgument exception if the scalar name ('name' parameter) is invalid. + */ + /* virtual double getOutputDoubleValue(const std::string& name) const; */ + + /*! @brief (Optional) Similar to setInputDoubleValue() but for an int value. + * @sa setInputDoubleValue() + */ + /* virtual void setInputIntValue(const std::string& name, const int& val); */ + + /*! @brief (Optional) Similar to getOutputDoubleValue() but for an int value. + * @sa getOutputDoubleValue() + */ + /* virtual int getOutputIntValue(const std::string& name) const; */ + + /*! @brief (Optional) Similar to setInputDoubleValue() but for an string value. + * @sa setInputDoubleValue() + */ + /* virtual void setInputStringValue(const std::string& name, const std::string& val); */ + + /*! @brief (Optional) Similar to getOutputDoubleValue() but for an string value. + * @sa getOutputDoubleValue() + */ + /* virtual std::string getOutputStringValue(const std::string& name) const; */ + + /*! @brief Return ICoCo interface major version number. + * @return ICoCo interface major version number (2 at present) + */ + static int GetICoCoMajorVersion() { return 2; } + + /*! @brief (Optional) Get MEDCoupling major version, if the code was built with MEDCoupling support. + * + * This can be used to assess compatibility between codes when coupling them. + * + * @return the MEDCoupling major version number (typically 7, 8, 9, ...) + */ + virtual int getMEDCouplingMajorVersion() const{ return MEDCOUPLING_VERSION_MAJOR; }; + + /*! @brief (Optional) Indicate whether the code was built with a 64-bits version of MEDCoupling. + * + * Implemented if the code was built with MEDCoupling support. + * This can be used to assess compatibility between codes when coupling them. + * + * @return the MEDCoupling major version number + */ + virtual bool isMEDCoupling64Bits() const; + + /*! @brief (Optional) Get the list of input scalars accepted by the code. + * + * @return the list of scalar names that represent inputs of the code + * @throws ICoCo::WrongContext exception if called before initialize() or after terminate(). + */ + /* virtual std::vector getInputValuesNames() const; */ + + /*! @brief (Optional) Get the list of output scalars that can be provided by the code. + * + * @return the list of scalar names that can be returned by the code + * @throws ICoCo::WrongContext exception if called before initialize() or after terminate(). + */ + /* virtual std::vector getOutputValuesNames() const; */ + Field getUnknownField() const; //paramètres du calcul -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* - /** \fn setPresentTime - * \brief met à jour _time (le temps courant du calcul) - * \details - * \param [in] double - * \param [out] void - * */ - void setPresentTime (double time); - /** \fn setMaxNbOfTimeStep * \brief met à jour _maxNbOfTimeStep ( le nombre maximum d'itération du calcul ) * \details @@ -574,6 +668,7 @@ public : void setHeatPowerField(Field heatPower){ _heatPowerField=heatPower; _heatPowerFieldSet=true; + _isStationary=false;//Source term may be changed after previously reaching a stationary state } /** \fn setHeatPowerField @@ -586,20 +681,22 @@ public : void setHeatPowerField(string fileName, string fieldName){ _heatPowerField=Field(fileName, CELLS,fieldName); _heatPowerFieldSet=true; + _isStationary=false;//Source term may be changed after previously reaching a stationary state } /** \fn setHeatSource - * \brief met à jour la puissance thermique ( _phi ) + * \brief sets a constant heat power field * \details * \param [in] double * \param [out] void * */ void setHeatSource(double phi){ _heatSource=phi; + _isStationary=false;//Source term may be changed after previously reaching a stationary state } /** \fn getHeatPowerField - * \brief renvoie le champs ?? ( _heatPowerField ) + * \brief returns the heat power field * \details * \param [in] void * \param [out] Field @@ -608,37 +705,6 @@ public : return _heatPowerField; } - /** \fn setRodTemperatureField ?? - * \brief - * \details - * \param [in] Field - * \param [out] void - * */ - void setRodTemperatureField(Field rodTemperature){ - _rodTemperatureField=rodTemperature; - _rodTemperatureFieldSet=true; - } - - /** \fn setRodTemperature ?? - * \brief - * \details - * \param [in] double - * \param [out] void - * */ - void setRodTemperature(double rodTemp){ - _rodTemperature=rodTemp; - } - - /** \fn getRodTemperatureField - * \brief - * \details - * \param [in] void - * \param [out] Field - * */ - virtual Field& getRodTemperatureField(){ // ?? je ne retrouve pas cet attribut dans le file.cxx - return _rodTemperatureField; - } - /** \fn setHeatTransfertCoeff * \brief set the heat transfert coefficient for heat exchange between fluid and solid * \details @@ -647,6 +713,7 @@ public : * */ void setHeatTransfertCoeff(double heatTransfertCoeff){ _heatTransfertCoeff=heatTransfertCoeff; + _isStationary=false;//Source term may be changed after previously reaching a stationary state } /** \fn setVerbose @@ -697,7 +764,7 @@ protected : int _Nmailles;//number of cells int _Nnodes;//number of nodes int _Nfaces;//number of faces - int _neibMaxNb;//maximum number of neighbours around a cell + int _neibMaxNbCells;//maximum number of neighbours around a cell int _neibMaxNbNodes;/* maximum number of nodes around a node */ Mesh _mesh; Field _perimeters; @@ -735,6 +802,7 @@ protected : bool _isStationary; bool _initialDataSet; bool _initializedMemory; + bool _stationaryMode;//ICoCo V2 bool _restartWithNewTimeScheme; bool _restartWithNewFileName; double _timeMax,_time; @@ -747,10 +815,10 @@ protected : ofstream * _runLogFile;//for creation of a log file to save the history of the simulation //Heat transfert variables - Field _heatPowerField, _rodTemperatureField; - bool _heatPowerFieldSet, _rodTemperatureFieldSet; + Field _heatPowerField; + bool _heatPowerFieldSet; double _heatTransfertCoeff; - double _heatSource, _rodTemperature; + double _heatSource; double _hsatv, _hsatl;//all models appart from DiffusionEquation will need this //Display variables @@ -759,11 +827,12 @@ protected : string _path;//path to execution directory used for saving results saveFormat _saveFormat;//file saving format : MED, VTK or CSV - //MPI variables - PetscMPIInt _size; /* size of communicator */ - PetscMPIInt _rank; /* processor rank */ - - + //MPI related variables + PetscMPIInt _mpi_size; /* size of communicator */ + PetscMPIInt _mpi_rank; /* processor rank */ + VecScatter _scat; /* For the distribution of a local vector */ + int _globalNbUnknowns, _localNbUnknowns; + int _d_nnz, _o_nnz; /* local and "non local" numbers of non zeros corfficients */ }; #endif /* PROBLEMCOREFLOWS_HXX_ */ diff --git a/CoreFlows/Models/inc/ProblemFluid.hxx b/CoreFlows/Models/inc/ProblemFluid.hxx index 9989b61..0858618 100755 --- a/CoreFlows/Models/inc/ProblemFluid.hxx +++ b/CoreFlows/Models/inc/ProblemFluid.hxx @@ -473,15 +473,16 @@ public : * */ void setNumericalScheme(SpaceScheme scheme, TimeScheme method=Explicit); - //données initiales + /* ICoCo code coupling interface */ /* - virtual vector getInputFieldsNames()=0 ;//Renvoie les noms des champs dont le problème a besoin (données initiales) virtual Field& getInputFieldTemplate(const string& name)=0;//Renvoie le format de champs attendu (maillage, composantes etc) - virtual void setInputField(const string& name, const Field& afield)=0;//enregistre les valeurs d'une donnée initiale virtual vector getOutputFieldsNames()=0 ;//liste tous les champs que peut fournir le code pour le postraitement virtual Field& getOutputField(const string& nameField )=0;//Renvoie un champs pour le postraitement */ - + /** Set input fields to prepare the simulation or coupling **/ + vector getInputFieldsNames(); + void setInputField(const string& nameField, Field& inputField );//supply of a required input field + Field getConservativeField() const ; Field getPrimitiveField() const; diff --git a/CoreFlows/Models/inc/SinglePhase.hxx b/CoreFlows/Models/inc/SinglePhase.hxx index f012b59..a1c19fa 100755 --- a/CoreFlows/Models/inc/SinglePhase.hxx +++ b/CoreFlows/Models/inc/SinglePhase.hxx @@ -146,7 +146,7 @@ public : double getReferencePressure() { return _Pref; }; double getReferenceTemperature() { return _Tref; }; - //get output fields for postprocessing or coupling + /* get output fields for postprocessing or coupling */ vector getOutputFieldsNames() ;//liste tous les champs que peut fournir le code pour le postraitement Field& getOutputField(const string& nameField );//Renvoie un champs pour le postraitement Field& getPressureField(); diff --git a/CoreFlows/Models/inc/StationaryDiffusionEquation.hxx b/CoreFlows/Models/inc/StationaryDiffusionEquation.hxx index ab0f1c3..b47aa6b 100755 --- a/CoreFlows/Models/inc/StationaryDiffusionEquation.hxx +++ b/CoreFlows/Models/inc/StationaryDiffusionEquation.hxx @@ -57,7 +57,6 @@ public : _fileName = fileName; } bool solveStationaryProblem(); - Field getOutputTemperatureField(); //Linear system and spectrum void setLinearSolver(linearSolver kspType, preconditioner pcType); @@ -106,6 +105,15 @@ public : void setConductivity(double conductivite){ _conductivity=conductivite; }; + void setDiffusiontensor(Matrix DiffusionTensor){ + _DiffusionTensor=DiffusionTensor; + }; + + + //get input fields to prepare the simulation or coupling + vector getInputFieldsNames(); + void setInputField(const string& nameField, Field& inputField );//supply of a required input field + void setFluidTemperatureField(Field coupledTemperatureField){ _fluidTemperatureField=coupledTemperatureField; _fluidTemperatureFieldSet=true; @@ -113,16 +121,9 @@ public : void setFluidTemperature(double fluidTemperature){ _fluidTemperature=fluidTemperature; } - Field& getRodTemperatureField(){ - return _VV; - } Field& getFluidTemperatureField(){ return _fluidTemperatureField; } - void setDiffusiontensor(Matrix DiffusionTensor){ - _DiffusionTensor=DiffusionTensor; - }; - /** \fn setHeatPowerField * \brief set the heat power field (variable in space) * \details @@ -146,6 +147,22 @@ public : _heatPowerFieldSet=true; } + /** \fn getHeatPowerField + * \brief returns the heat power field + * \details + * \param [in] void + * \param [out] Field + * */ + Field getHeatPowerField(){ + return _heatPowerField; + } + //get output fields names for postprocessing or coupling + vector getOutputFieldsNames() ;//liste tous les champs que peut fournir le code pour le postraitement + Field& getOutputField(const string& nameField );//Renvoie un champs pour le postraitement + + Field& getOutputTemperatureField(); + Field& getRodTemperatureField(); + /** \fn setVerbose * \brief Updates display options * \details @@ -200,6 +217,8 @@ protected : Vector _normale; Matrix _DiffusionTensor; Vec _deltaT, _Tk, _Tkm1, _b0; + Vec _Tk_seq; // Local sequential copy of the parallel vector _Tk, used for saving result files + double _dt_src; //Heat transfert variables @@ -236,9 +255,12 @@ protected : std::map< int, double> _dirichletBoundaryValues; std::map< int, double> _neumannBoundaryValues; - //MPI variables - PetscMPIInt _size; /* size of communicator */ - PetscMPIInt _rank; /* processor rank */ + /**** MPI related variables ***/ + PetscMPIInt _mpi_size; /* size of communicator */ + PetscMPIInt _mpi_rank; /* processor rank */ + VecScatter _scat; /* For the distribution of a local vector */ + int _globalNbUnknowns, _localNbUnknowns; + int _d_nnz, _o_nnz; /* local and "non local" numbers of non zeros corfficients */ }; #endif /* StationaryDiffusionEquation_HXX_ */ diff --git a/CoreFlows/Models/inc/TransportEquation.hxx b/CoreFlows/Models/inc/TransportEquation.hxx index 6f17e7f..94f01c0 100755 --- a/CoreFlows/Models/inc/TransportEquation.hxx +++ b/CoreFlows/Models/inc/TransportEquation.hxx @@ -126,7 +126,44 @@ public : _vitesseTransport=v; }; - //get output fields for postprocessing or coupling + /* set input fields to prepare the simulation */ + vector getInputFieldsNames(); + void setInputField(const string& nameField, Field& inputField );//supply of a required input field + + /** \fn setRodTemperatureField + * \brief Set the rod temperature field + * \details + * \param [in] Field + * \param [out] void + * */ + void setRodTemperatureField(Field rodTemperature){ + _rodTemperatureField=rodTemperature; + _rodTemperatureFieldSet=true; + _isStationary=false;//Source term may be changed after previously reaching a stationary state + } + + /** \fn setRodTemperature + * \brief Set a constant rod temperature field + * \details + * \param [in] double + * \param [out] void + * */ + void setRodTemperature(double rodTemp){ + _rodTemperature=rodTemp; + _isStationary=false;//Source term may be changed after previously reaching a stationary state + } + + /** \fn getRodTemperatureField + * \brief + * \details + * \param [in] void + * \param [out] Field + * */ + Field& getRodTemperatureField(){ // ?? je ne retrouve pas cet attribut dans le file.cxx + return _rodTemperatureField; + } + + /* get output fields for postprocessing or coupling */ vector getOutputFieldsNames() ;//liste tous les champs que peut fournir le code pour le postraitement Field& getOutputField(const string& nameField );//Renvoie un champs pour le postraitement @@ -138,10 +175,24 @@ public : return _VV; } + Field& getVoidFractionField(){ + return _Alpha; + } + + Field& getDensityField(){ + return _Rho; + } + protected : double computeTransportMatrix(); double computeRHS(); void updatePrimitives(); + + /* Postprocessing fields */ + Field _TT, _Alpha, _Rho;//Fields of temperature, void fraction, density. Unknown field is enthalpy (_VV) + double _rhosatv, _rhosatl; + double _Tref, _href, _cpref; + double temperature(double h){ return _Tref+(h-_href)/_cpref; }; @@ -157,15 +208,18 @@ protected : return alpha*_rhosatv+(1-alpha)*_rhosatl; }; - Field _TT, _Alpha, _Rho;//Fields of temperature and coupled temperature - double _rhosatv, _rhosatl; - double _Tref, _href, _cpref; Vector _vitesseTransport, _normale; bool _transportMatrixSet; Vec _Hn, _deltaH, _Hk, _Hkm1, _b0; + Vec _Hn_seq; // Local sequential copy of the parallel vector _Hn, used for saving result files double _dt_transport, _dt_src; map _limitField; + + /* source terms */ + bool _rodTemperatureFieldSet; + Field _rodTemperatureField; + double _rodTemperature; }; #endif /* TransportEquation_HXX_ */ diff --git a/CoreFlows/Models/src/DiffusionEquation.cxx b/CoreFlows/Models/src/DiffusionEquation.cxx index b8993c6..51c78ea 100755 --- a/CoreFlows/Models/src/DiffusionEquation.cxx +++ b/CoreFlows/Models/src/DiffusionEquation.cxx @@ -130,126 +130,146 @@ DiffusionEquation::DiffusionEquation(int dim, bool FECalculation,double rho,doub void DiffusionEquation::initialize() { - _runLogFile->open((_fileName+".log").c_str(), ios::out | ios::trunc);;//for creation of a log file to save the history of the simulation - - if(_Ndim != _mesh.getSpaceDimension() or _Ndim!=_mesh.getMeshDimension())//for the moment we must have space dim=mesh dim + if(_mpi_rank==0) { - PetscPrintf(PETSC_COMM_WORLD,"Problem : dimension defined is %d but mesh dimension= %d, and space dimension is %d",_Ndim,_mesh.getMeshDimension(),_mesh.getSpaceDimension()); - *_runLogFile<< "Problem : dim = "<<_Ndim<< " but mesh dim= "<<_mesh.getMeshDimension()<<", mesh space dim= "<<_mesh.getSpaceDimension()<close(); - throw CdmathException("DiffusionEquation::initialize: mesh has incorrect dimension"); - } - - if(!_initialDataSet) - throw CdmathException("DiffusionEquation::initialize() set initial data first"); - else - { - PetscPrintf(PETSC_COMM_WORLD,"Initialising the diffusion of a solid temperature using "); - *_runLogFile<<"Initialising the diffusion of a solid temperature using "; - if(!_FECalculation) - { - PetscPrintf(PETSC_COMM_WORLD,"Finite volumes method\n\n"); - *_runLogFile<< "Finite volumes method"<open((_fileName+".log").c_str(), ios::out | ios::trunc);;//for creation of a log file to save the history of the simulation + + if(_Ndim != _mesh.getSpaceDimension() or _Ndim!=_mesh.getMeshDimension())//for the moment we must have space dim=mesh dim + { + PetscPrintf(PETSC_COMM_SELF,"Problem : dimension defined is %d but mesh dimension= %d, and space dimension is %d",_Ndim,_mesh.getMeshDimension(),_mesh.getSpaceDimension()); + *_runLogFile<< "Problem : dim = "<<_Ndim<< " but mesh dim= "<<_mesh.getMeshDimension()<<", mesh space dim= "<<_mesh.getSpaceDimension()<close(); + throw CdmathException("DiffusionEquation::initialize: mesh has incorrect dimension"); + } + + if(!_initialDataSet) + throw CdmathException("DiffusionEquation::initialize() set initial data first"); + else + { + PetscPrintf(PETSC_COMM_SELF,"Initialising the diffusion of a solid temperature using "); + *_runLogFile<<"Initialising the diffusion of a solid temperature using "; + if(!_FECalculation) + { + PetscPrintf(PETSC_COMM_SELF,"Finite volumes method\n\n"); + *_runLogFile<< "Finite volumes method"<close(); + throw CdmathException("DiffusionEquation::initialize: mesh has incorrect cell types"); + } + } + else if(_Ndim==3) + { + if( _mesh.isTetrahedral() )//Mesh dim=3 + PetscPrintf(PETSC_COMM_SELF,"3D Finite element method on tetrahedra\n"); + else if (_mesh.getMeshDimension()==2 and _mesh.isTriangular())//Mesh dim=2 + PetscPrintf(PETSC_COMM_SELF,"2D Finite element method on a 3D surface : space dimension is %d, mesh dimension is %d\n",_Ndim,_mesh.getMeshDimension()); + else if (_mesh.getMeshDimension()==1)//Mesh dim=1 + PetscPrintf(PETSC_COMM_SELF,"1D Finite element method on a 3D network : space dimension is %d, mesh dimension is %d\n",_Ndim,_mesh.getMeshDimension()); + else + { + PetscPrintf(PETSC_COMM_SELF,"Error Finite element with space dimension %d, and mesh dimension %d, mesh should be either tetrahedral, either a triangularised surface or 1D network",_Ndim,_mesh.getMeshDimension()); + *_runLogFile<<"DiffusionEquation::initialize mesh has incorrect dimension"<close(); + throw CdmathException("DiffusionEquation::initialize: mesh has incorrect cell types"); + } + } + + _boundaryNodeIds = _mesh.getBoundaryNodeIds(); + _NboundaryNodes=_boundaryNodeIds.size(); + + if(_NboundaryNodes==_Nnodes) + PetscPrintf(PETSC_COMM_SELF,"!!!!! Warning : all nodes are boundary nodes !!!!!"); + + for(int i=0; i<_NboundaryNodes; i++) + if(_limitField[(_mesh.getNode(_boundaryNodeIds[i])).getGroupName()].bcType==DirichletDiffusion) + _dirichletNodeIds.push_back(_boundaryNodeIds[i]); + _NdirichletNodes=_dirichletNodeIds.size(); + _NunknownNodes=_Nnodes - _NdirichletNodes; + PetscPrintf(PETSC_COMM_SELF,"Number of unknown nodes %d, Number of boundary nodes %d, Number of Dirichlet boundary nodes %d\n\n", _NunknownNodes,_NboundaryNodes, _NdirichletNodes); + *_runLogFile<<"Number of unknown nodes " << _NunknownNodes <<", Number of boundary nodes " << _NboundaryNodes<< ", Number of Dirichlet boundary nodes " << _NdirichletNodes <close(); - throw CdmathException("DiffusionEquation::initialize: mesh has incorrect cell types"); - } - } - else if(_Ndim==3) - { - if( _mesh.isTetrahedral() )//Mesh dim=3 - PetscPrintf(PETSC_COMM_WORLD,"3D Finite element method on tetrahedra\n"); - else if (_mesh.getMeshDimension()==2 and _mesh.isTriangular())//Mesh dim=2 - PetscPrintf(PETSC_COMM_WORLD,"2D Finite element method on a 3D surface : space dimension is %d, mesh dimension is %d\n",_Ndim,_mesh.getMeshDimension()); - else if (_mesh.getMeshDimension()==1)//Mesh dim=1 - PetscPrintf(PETSC_COMM_WORLD,"1D Finite element method on a 3D network : space dimension is %d, mesh dimension is %d\n",_Ndim,_mesh.getMeshDimension()); - else - { - PetscPrintf(PETSC_COMM_WORLD,"Error Finite element with space dimension %d, and mesh dimension %d, mesh should be either tetrahedral, either a triangularised surface or 1D network",_Ndim,_mesh.getMeshDimension()); - *_runLogFile<<"DiffusionEquation::initialize mesh has incorrect dimension"<close(); - throw CdmathException("DiffusionEquation::initialize: mesh has incorrect cell types"); - } - } - - _boundaryNodeIds = _mesh.getBoundaryNodeIds(); - _NboundaryNodes=_boundaryNodeIds.size(); - - if(_NboundaryNodes==_Nnodes) - PetscPrintf(PETSC_COMM_WORLD,"!!!!! Warning : all nodes are boundary nodes !!!!!"); - - for(int i=0; i<_NboundaryNodes; i++) - if(_limitField[(_mesh.getNode(_boundaryNodeIds[i])).getGroupName()].bcType==DirichletDiffusion) - _dirichletNodeIds.push_back(_boundaryNodeIds[i]); - _NdirichletNodes=_dirichletNodeIds.size(); - _NunknownNodes=_Nnodes - _NdirichletNodes; - PetscPrintf(PETSC_COMM_WORLD,"Number of unknown nodes %d, Number of boundary nodes %d, Number of Dirichlet boundary nodes %d\n\n", _NunknownNodes,_NboundaryNodes, _NdirichletNodes); - *_runLogFile<<"Number of unknown nodes " << _NunknownNodes <<", Number of boundary nodes " << _NboundaryNodes<< ", Number of Dirichlet boundary nodes " << _NdirichletNodes <1 && _mpi_rank == 0) + VecCreateSeq(PETSC_COMM_SELF, _globalNbUnknowns, &_Tn_seq);//For saving results on proc 0 + VecScatterCreateToZero(_Tn,&_scat,&_Tn_seq); //Linear solver - KSPCreate(PETSC_COMM_SELF, &_ksp); + KSPCreate(PETSC_COMM_WORLD, &_ksp); KSPSetType(_ksp, _ksptype); // if(_ksptype == KSPGMRES) KSPGMRESSetRestart(_ksp,10000); KSPSetTolerances(_ksp,_precision,_precision,PETSC_DEFAULT,_maxPetscIts); @@ -269,6 +289,15 @@ double DiffusionEquation::computeTimeStep(bool & stop){ _dt_src=computeRHS(stop); + VecAssemblyBegin(_b); + VecAssemblyEnd( _b); + + if(_verbose or _system) + { + PetscPrintf(PETSC_COMM_WORLD,"Right hand side of the linear system\n"); + VecView(_b,PETSC_VIEWER_STDOUT_WORLD); + } + stop=false; return min(_dt_diffusion,_dt_src); } @@ -277,131 +306,136 @@ double DiffusionEquation::computeDiffusionMatrix(bool & stop) { double result; + MatZeroEntries(_A); + VecZeroEntries(_b0); + if(_FECalculation) result=computeDiffusionMatrixFE(stop); else result=computeDiffusionMatrixFV(stop); + PetscPrintf(PETSC_COMM_WORLD,"Maximum diffusivity is %.2f, CFL = %.2f, Delta x = %.2f\n",_maxvp,_cfl,_minl); + + MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY); + MatAssemblyEnd( _A, MAT_FINAL_ASSEMBLY); + VecAssemblyBegin(_b0); + VecAssemblyEnd( _b0); + + //Contribution from the solid/fluid heat exchange with assumption of constant heat transfer coefficient //update value here if variable heat transfer coefficient if(_timeScheme == Implicit and _heatTransfertCoeff/(_rho*_cp)>_precision) MatShift(_A,_heatTransfertCoeff/(_rho*_cp));//Contribution from the liquit/solid heat transfer if(_verbose or _system) - MatView(_A,PETSC_VIEWER_STDOUT_SELF); + MatView(_A,PETSC_VIEWER_STDOUT_WORLD); return result; } double DiffusionEquation::computeDiffusionMatrixFE(bool & stop){ - Cell Cj; - string nameOfGroup; - double coeff;//Diffusion coefficients between nodes i and j - MatZeroEntries(_A); - VecZeroEntries(_b); - - Matrix M(_Ndim+1,_Ndim+1);//cell geometry matrix - std::vector< Vector > GradShapeFuncs(_Ndim+1);//shape functions of cell nodes - std::vector< int > nodeIds(_Ndim+1);//cell node Ids - std::vector< Node > nodes(_Ndim+1);//cell nodes - int i_int, j_int; //index of nodes j and k considered as unknown nodes - bool dirichletCell_treated; - - std::vector< vector< double > > values(_Ndim+1,vector< double >(_Ndim+1,0));//values of shape functions on cell node - for (int idim=0; idim<_Ndim+1;idim++) - values[idim][idim]=1; - - /* parameters for boundary treatment */ - vector< double > valuesBorder(_Ndim+1); - Vector GradShapeFuncBorder(_Ndim+1); - - for (int j=0; j<_Nmailles;j++) - { - Cj = _mesh.getCell(j); - - for (int idim=0; idim<_Ndim+1;idim++){ - nodeIds[idim]=Cj.getNodeId(idim); - nodes[idim]=_mesh.getNode(nodeIds[idim]); - for (int jdim=0; jdim<_Ndim;jdim++) - M(idim,jdim)=nodes[idim].getPoint()[jdim]; - M(idim,_Ndim)=1; - } - for (int idim=0; idim<_Ndim+1;idim++) - GradShapeFuncs[idim]=gradientNodal(M,values[idim])/fact(_Ndim); - - /* Loop on the edges of the cell */ - for (int idim=0; idim<_Ndim+1;idim++) - { - if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),nodeIds[idim])==_dirichletNodeIds.end())//!_mesh.isBorderNode(nodeIds[idim]) - {//First node of the edge is not Dirichlet node - i_int=unknownNodeIndex(nodeIds[idim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing - dirichletCell_treated=false; - for (int jdim=0; jdim<_Ndim+1;jdim++) - { - if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),nodeIds[jdim])==_dirichletNodeIds.end())//!_mesh.isBorderNode(nodeIds[jdim]) - {//Second node of the edge is not Dirichlet node - j_int= unknownNodeIndex(nodeIds[jdim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing - MatSetValue(_A,i_int,j_int,(_DiffusionTensor*GradShapeFuncs[idim])*GradShapeFuncs[jdim]/Cj.getMeasure(), ADD_VALUES); - } - else if (!dirichletCell_treated) - {//Second node of the edge is a Dirichlet node - dirichletCell_treated=true; - for (int kdim=0; kdim<_Ndim+1;kdim++) - { - std::map::iterator it=_dirichletBoundaryValues.find(nodeIds[kdim]); - if( it != _dirichletBoundaryValues.end() ) - { - if( _dirichletValuesSet ) - valuesBorder[kdim]=_dirichletBoundaryValues[it->second]; - else - valuesBorder[kdim]=_limitField[_mesh.getNode(nodeIds[kdim]).getGroupName()].T; - } - else - valuesBorder[kdim]=0; - } - GradShapeFuncBorder=gradientNodal(M,valuesBorder)/fact(_Ndim); - coeff =-1.*(_DiffusionTensor*GradShapeFuncBorder)*GradShapeFuncs[idim]/Cj.getMeasure(); - VecSetValue(_b,i_int,coeff, ADD_VALUES); - } - } - } - } + if(_mpi_rank == 0) + { + Cell Cj; + string nameOfGroup; + double coeff;//Diffusion coefficients between nodes i and j + + Matrix M(_Ndim+1,_Ndim+1);//cell geometry matrix + std::vector< Vector > GradShapeFuncs(_Ndim+1);//shape functions of cell nodes + std::vector< int > nodeIds(_Ndim+1);//cell node Ids + std::vector< Node > nodes(_Ndim+1);//cell nodes + int i_int, j_int; //index of nodes j and k considered as unknown nodes + bool dirichletCell_treated; + + std::vector< vector< double > > values(_Ndim+1,vector< double >(_Ndim+1,0));//values of shape functions on cell node + for (int idim=0; idim<_Ndim+1;idim++) + values[idim][idim]=1; + + /* parameters for boundary treatment */ + vector< double > valuesBorder(_Ndim+1); + Vector GradShapeFuncBorder(_Ndim+1); + + for (int j=0; j<_Nmailles;j++) + { + Cj = _mesh.getCell(j); + + for (int idim=0; idim<_Ndim+1;idim++){ + nodeIds[idim]=Cj.getNodeId(idim); + nodes[idim]=_mesh.getNode(nodeIds[idim]); + for (int jdim=0; jdim<_Ndim;jdim++) + M(idim,jdim)=nodes[idim].getPoint()[jdim]; + M(idim,_Ndim)=1; + } + for (int idim=0; idim<_Ndim+1;idim++) + GradShapeFuncs[idim]=gradientNodal(M,values[idim])/fact(_Ndim); + + /* Loop on the edges of the cell */ + for (int idim=0; idim<_Ndim+1;idim++) + { + if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),nodeIds[idim])==_dirichletNodeIds.end())//!_mesh.isBorderNode(nodeIds[idim]) + {//First node of the edge is not Dirichlet node + i_int=unknownNodeIndex(nodeIds[idim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing + dirichletCell_treated=false; + for (int jdim=0; jdim<_Ndim+1;jdim++) + { + if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),nodeIds[jdim])==_dirichletNodeIds.end())//!_mesh.isBorderNode(nodeIds[jdim]) + {//Second node of the edge is not Dirichlet node + j_int= unknownNodeIndex(nodeIds[jdim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing + MatSetValue(_A,i_int,j_int,(_DiffusionTensor*GradShapeFuncs[idim])*GradShapeFuncs[jdim]/Cj.getMeasure(), ADD_VALUES); + } + else if (!dirichletCell_treated) + {//Second node of the edge is a Dirichlet node + dirichletCell_treated=true; + for (int kdim=0; kdim<_Ndim+1;kdim++) + { + std::map::iterator it=_dirichletBoundaryValues.find(nodeIds[kdim]); + if( it != _dirichletBoundaryValues.end() ) + { + if( _dirichletValuesSet ) + valuesBorder[kdim]=_dirichletBoundaryValues[it->second]; + else + valuesBorder[kdim]=_limitField[_mesh.getNode(nodeIds[kdim]).getGroupName()].T; + } + else + valuesBorder[kdim]=0; + } + GradShapeFuncBorder=gradientNodal(M,valuesBorder)/fact(_Ndim); + coeff =-1.*(_DiffusionTensor*GradShapeFuncBorder)*GradShapeFuncs[idim]/Cj.getMeasure(); + VecSetValue(_b0,i_int,coeff, ADD_VALUES); + } + } + } + } + } + + //Calcul de la contribution de la condition limite de Neumann au second membre + if( _NdirichletNodes !=_NboundaryNodes) + { + vector< int > boundaryFaces = _mesh.getBoundaryFaceIds(); + int NboundaryFaces=boundaryFaces.size(); + for(int i = 0; i< NboundaryFaces ; i++)//On parcourt les faces du bord + { + Face Fi = _mesh.getFace(boundaryFaces[i]); + for(int j = 0 ; j<_Ndim ; j++)//On parcourt les noeuds de la face + { + if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),Fi.getNodeId(j))==_dirichletNodeIds.end())//node j is a Neumann BC node (not a Dirichlet BC node) + { + j_int=unknownNodeIndex(Fi.getNodeId(j), _dirichletNodeIds);//indice du noeud j en tant que noeud inconnu + if( _neumannValuesSet ) + coeff =Fi.getMeasure()/_Ndim*_neumannBoundaryValues[Fi.getNodeId(j)]; + else + coeff =Fi.getMeasure()/_Ndim*_limitField[_mesh.getNode(Fi.getNodeId(j)).getGroupName()].normalFlux; + VecSetValue(_b0, j_int, coeff, ADD_VALUES); + } + } + } + } } - - //Calcul de la contribution de la condition limite de Neumann au second membre - if( _NdirichletNodes !=_NboundaryNodes) - { - vector< int > boundaryFaces = _mesh.getBoundaryFaceIds(); - int NboundaryFaces=boundaryFaces.size(); - for(int i = 0; i< NboundaryFaces ; i++)//On parcourt les faces du bord - { - Face Fi = _mesh.getFace(boundaryFaces[i]); - for(int j = 0 ; j<_Ndim ; j++)//On parcourt les noeuds de la face - { - if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),Fi.getNodeId(j))==_dirichletNodeIds.end())//node j is a Neumann BC node (not a Dirichlet BC node) - { - j_int=unknownNodeIndex(Fi.getNodeId(j), _dirichletNodeIds);//indice du noeud j en tant que noeud inconnu - if( _neumannValuesSet ) - coeff =Fi.getMeasure()/_Ndim*_neumannBoundaryValues[Fi.getNodeId(j)]; - else - coeff =Fi.getMeasure()/_Ndim*_limitField[_mesh.getNode(Fi.getNodeId(j)).getGroupName()].normalFlux; - VecSetValue(_b, j_int, coeff, ADD_VALUES); - } - } - } - } - - MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY); - MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY); - VecAssemblyBegin(_b); - VecAssemblyEnd(_b); _diffusionMatrixSet=true; - stop=false ; _maxvp=_diffusivity;//To do : optimise value with the mesh while respecting stability - PetscPrintf(PETSC_COMM_SELF,"Maximum diffusivity is %.2f, CFL = %.2f, Delta x = %.2f\n",_maxvp,_cfl,_minl); if(fabs(_maxvp)<_precision) throw CdmathException("DiffusionEquation::computeDiffusionMatrixFE(): Error computing time step ! Maximum diffusivity is zero => division by zero"); else @@ -409,107 +443,103 @@ double DiffusionEquation::computeDiffusionMatrixFE(bool & stop){ } double DiffusionEquation::computeDiffusionMatrixFV(bool & stop){ - long nbFaces = _mesh.getNumberOfFaces(); - Face Fj; - Cell Cell1,Cell2; - string nameOfGroup; - double inv_dxi, inv_dxj; - double barycenterDistance; - Vector normale(_Ndim); - double dn; - PetscInt idm, idn; - std::vector< int > idCells; - MatZeroEntries(_A); - VecZeroEntries(_b0); - - for (int j=0; j_maxvp) - _maxvp=fabs(dn); - - // compute 1/dxi = volume of Ci/area of Fj - inv_dxi = Fj.getMeasure()/Cell1.getMeasure(); - - // If Fj is on the boundary - if (Fj.getNumberOfCells()==1) { - if(_verbose ) - { - cout << "face numero " << j << " cellule frontiere " << idCells[0] << " ; vecteur normal=("; - for(int p=0; p<_Ndim; p++) - cout << normale[p] << ","; - cout << ") "< idCells; + + for (int j=0; j_maxvp) + _maxvp=fabs(dn); + + // compute 1/dxi = volume of Ci/area of Fj + inv_dxi = Fj.getMeasure()/Cell1.getMeasure(); + + // If Fj is on the boundary + if (Fj.getNumberOfCells()==1) { + if(_verbose ) + { + cout << "face numero " << j << " cellule frontiere " << idCells[0] << " ; vecteur normal=("; + for(int p=0; p<_Ndim; p++) + cout << normale[p] << ","; + cout << ") "< 1) + inv_dxj = Fj.getMeasure()/Cell2.getMeasure(); + else + inv_dxj = 1/Cell2.getMeasure(); + + barycenterDistance=Cell1.getBarryCenter().distance(Cell2.getBarryCenter()); + + MatSetValue(_A,idm,idm, dn*inv_dxi/barycenterDistance, ADD_VALUES); + MatSetValue(_A,idm,idn,-dn*inv_dxi/barycenterDistance, ADD_VALUES); + MatSetValue(_A,idn,idn, dn*inv_dxj/barycenterDistance, ADD_VALUES); + MatSetValue(_A,idn,idm,-dn*inv_dxj/barycenterDistance, ADD_VALUES); } - Cell2 = _mesh.getCell(idCells[1]); - idn = idCells[1]; - if (_Ndim > 1) - inv_dxj = Fj.getMeasure()/Cell2.getMeasure(); else - inv_dxj = 1/Cell2.getMeasure(); - - barycenterDistance=Cell1.getBarryCenter().distance(Cell2.getBarryCenter()); - - MatSetValue(_A,idm,idm, dn*inv_dxi/barycenterDistance, ADD_VALUES); - MatSetValue(_A,idm,idn,-dn*inv_dxi/barycenterDistance, ADD_VALUES); - MatSetValue(_A,idn,idn, dn*inv_dxj/barycenterDistance, ADD_VALUES); - MatSetValue(_A,idn,idm,-dn*inv_dxj/barycenterDistance, ADD_VALUES); + throw CdmathException("DiffusionEquation::computeDiffusionMatrixFV(): incompatible number of cells around a face"); } - else - throw CdmathException("DiffusionEquation::computeDiffusionMatrixFV(): incompatible number of cells around a face"); } - MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY); - MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY); - VecAssemblyBegin(_b0); - VecAssemblyEnd(_b0); - _diffusionMatrixSet=true; - stop=false ; - cout<<"Maximum diffusivity is "<<_maxvp<< " CFL = "<<_cfl<<" Delta x = "<<_minl< division by zero"); else @@ -517,42 +547,38 @@ double DiffusionEquation::computeDiffusionMatrixFV(bool & stop){ } double DiffusionEquation::computeRHS(bool & stop){//Contribution of the PDE RHS to the linear systemm RHS (boundary conditions do contribute to the system RHS via the function computeDiffusionMatrix - VecAssemblyBegin(_b); - double Ti; - if(!_FECalculation) - for (int i=0; i<_Nmailles;i++) - //Contribution due to fluid/solide heat exchange + Contribution of the volumic heat power - if(_timeScheme == Explicit) - { - VecGetValues(_Tn, 1, &i, &Ti); - VecSetValue(_b,i,(_heatTransfertCoeff*(_fluidTemperatureField(i)-Ti)+_heatPowerField(i))/(_rho*_cp),ADD_VALUES); - } - else//Implicit scheme - VecSetValue(_b,i,(_heatTransfertCoeff* _fluidTemperatureField(i) +_heatPowerField(i))/(_rho*_cp) ,ADD_VALUES); - else - { - Cell Ci; - std::vector< int > nodesId; - for (int i=0; i<_Nmailles;i++) - { - Ci=_mesh.getCell(i); - nodesId=Ci.getNodesId(); - for (int j=0; jdirichletNodes.upper_bound() - { - double coeff = (_heatTransfertCoeff*_fluidTemperatureField(nodesId[j]) + _heatPowerField(nodesId[j]))/(_rho*_cp); - VecSetValue(_b,unknownNodeIndex(nodesId[j], _dirichletNodeIds), coeff*Ci.getMeasure()/(_Ndim+1),ADD_VALUES); - } - } - } - VecAssemblyEnd(_b); - if(_verbose or _system) + if(_mpi_rank == 0) { - PetscPrintf(PETSC_COMM_WORLD,"Right hand side of the linear system\n"); - VecView(_b,PETSC_VIEWER_STDOUT_SELF); + double Ti; + if(!_FECalculation) + for (int i=0; i<_Nmailles;i++) + //Contribution due to fluid/solide heat exchange + Contribution of the volumic heat power + if(_timeScheme == Explicit) + { + VecGetValues(_Tn, 1, &i, &Ti); + VecSetValue(_b,i,(_heatTransfertCoeff*(_fluidTemperatureField(i)-Ti)+_heatPowerField(i))/(_rho*_cp),ADD_VALUES); + } + else//Implicit scheme + VecSetValue(_b,i,(_heatTransfertCoeff* _fluidTemperatureField(i) +_heatPowerField(i))/(_rho*_cp) ,ADD_VALUES); + else + { + Cell Ci; + std::vector< int > nodesId; + for (int i=0; i<_Nmailles;i++) + { + Ci=_mesh.getCell(i); + nodesId=Ci.getNodesId(); + for (int j=0; jdirichletNodes.upper_bound() + { + double coeff = (_heatTransfertCoeff*_fluidTemperatureField(nodesId[j]) + _heatPowerField(nodesId[j]))/(_rho*_cp); + VecSetValue(_b,unknownNodeIndex(nodesId[j], _dirichletNodeIds), coeff*Ci.getMeasure()/(_Ndim+1),ADD_VALUES); + } + } + } } - stop=false ; + if(_heatTransfertCoeff>_precision) return _rho*_cp/_heatTransfertCoeff; else @@ -561,19 +587,18 @@ double DiffusionEquation::computeRHS(bool & stop){//Contribution of the PDE RHS bool DiffusionEquation::initTimeStep(double dt){ - if(_dt>0 and dt>0) + /* tricky because of code coupling */ + if(_dt>0 and dt>0)//Previous time step was set and used { - //Remove the contribution from dt to prepare for new initTimeStep. The diffusion matrix is not recomputed + //Remove the contribution from dt to prepare for new time step. The diffusion matrix is not recomputed if(_timeScheme == Implicit) MatShift(_A,-1/_dt+1/dt); //No need to remove the contribution to the right hand side since it is recomputed from scratch at each time step } else if(dt>0)//_dt==0, first time step { - MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY); - MatAssemblyEnd( _A, MAT_FINAL_ASSEMBLY); if(_timeScheme == Implicit) - MatShift(_A,1/_dt); + MatShift(_A,1/dt); } else//dt<=0 { @@ -581,25 +606,23 @@ bool DiffusionEquation::initTimeStep(double dt){ throw CdmathException("Error DiffusionEquation::initTimeStep : cannot set time step to zero"); } //At this stage _b contains _b0 + power + heat exchange - VecAXPY(_b, 1/_dt, _Tn); + VecAXPY(_b, 1/dt, _Tn); _dt = dt; if(_verbose && (_nbTimeStep-1)%_freqSave ==0) { PetscPrintf(PETSC_COMM_WORLD,"Matrix of the linear system\n"); - MatView(_A,PETSC_VIEWER_STDOUT_SELF); + MatView(_A,PETSC_VIEWER_STDOUT_WORLD); } return _dt>0; } void DiffusionEquation::abortTimeStep(){ - //Remove contribution od dt to the RHS + //Remove contribution of dt to the RHS VecAXPY(_b, -1/_dt, _Tn); - MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY); - MatAssemblyEnd( _A, MAT_FINAL_ASSEMBLY); - //Remove contribution od dt to the matrix + //Remove contribution of dt to the matrix if(_timeScheme == Implicit) MatShift(_A,-1/_dt); _dt = 0; @@ -628,8 +651,6 @@ bool DiffusionEquation::iterateTimeStep(bool &converged) } else { - MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY); - MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY); #if PETSC_VERSION_GREATER_3_5 KSPSetOperators(_ksp, _A, _A); @@ -655,16 +676,7 @@ bool DiffusionEquation::iterateTimeStep(bool &converged) { VecCopy(_Tk, _deltaT);//ici on a deltaT=Tk VecAXPY(_deltaT, -1, _Tkm1);//On obtient deltaT=Tk-Tkm1 - _erreur_rel= 0; - double Ti, dTi; - - for(int i=0; i<_VV.getNumberOfElements(); i++) - { - VecGetValues(_deltaT, 1, &i, &dTi); - VecGetValues(_Tk, 1, &i, &Ti); - if(_erreur_rel < fabs(dTi/Ti)) - _erreur_rel = fabs(dTi/Ti); - } + VecNorm(_deltaT,NORM_INFINITY,&_erreur_rel); converged = (_erreur_rel <= _precision) ;//converged=convergence des iterations de Newton } } @@ -677,34 +689,22 @@ void DiffusionEquation::validateTimeStep() { VecCopy(_Tk, _deltaT);//ici Tk=Tnp1 donc on a deltaT=Tnp1 VecAXPY(_deltaT, -1, _Tn);//On obtient deltaT=Tnp1-Tn - - _erreur_rel= 0; - double Ti, dTi; - - for(int i=0; i<_VV.getNumberOfElements(); i++) - { - VecGetValues(_deltaT, 1, &i, &dTi); - VecGetValues(_Tk, 1, &i, &Ti); - if(_erreur_rel < fabs(dTi/Ti)) - _erreur_rel = fabs(dTi/Ti); - } + VecNorm(_deltaT,NORM_INFINITY,&_erreur_rel); _isStationary =(_erreur_rel <_precision); VecCopy(_Tk, _Tn); VecCopy(_Tk, _Tkm1); - if(_verbose && (_nbTimeStep-1)%_freqSave ==0) - PetscPrintf(PETSC_COMM_WORLD,"Valeur propre locale max: %.2f\n", _maxvp); - _time+=_dt; _nbTimeStep++; + if ((_nbTimeStep-1)%_freqSave ==0 || _isStationary || _time>=_timeMax || _nbTimeStep>=_maxNbOfTimeStep) - save(); + save(); } void DiffusionEquation::save(){ - PetscPrintf(PETSC_COMM_SELF,"Saving numerical results\n\n"); + PetscPrintf(PETSC_COMM_WORLD,"Saving numerical results\n\n"); *_runLogFile<< "Saving numerical results"<< endl<1){ + VecScatterBegin(_scat,_Tn,_Tn_seq,INSERT_VALUES,SCATTER_FORWARD); + VecScatterEnd( _scat,_Tn,_Tn_seq,INSERT_VALUES,SCATTER_FORWARD); + } + if(_verbose or _system) { PetscPrintf(PETSC_COMM_WORLD,"Unknown of the linear system :\n"); - VecView(_Tk,PETSC_VIEWER_STDOUT_SELF); + VecView(_Tn,PETSC_VIEWER_STDOUT_WORLD); } - //On remplit le champ - double Ti; - if(!_FECalculation) - for(int i =0; i<_Nmailles;i++) - { - VecGetValues(_Tn, 1, &i, &Ti); - _VV(i)=Ti; - } - else - { - int globalIndex; - for(int i=0; i<_NunknownNodes; i++) - { - VecGetValues(_Tk, 1, &i, &Ti); - globalIndex = globalNodeIndex(i, _dirichletNodeIds); - _VV(globalIndex)=Ti;//Assumes node numbering starts with border nodes - } - - Node Ni; - string nameOfGroup; - for(int i=0; i<_NdirichletNodes; i++)//Assumes node numbering starts with border nodes - { - Ni=_mesh.getNode(_dirichletNodeIds[i]); - nameOfGroup = Ni.getGroupName(); - _VV(_dirichletNodeIds[i])=_limitField[nameOfGroup].T; - } - } - _VV.setTime(_time,_nbTimeStep); - - // create mesh and component info - if (_nbTimeStep ==0 || _restartWithNewFileName){ - if (_restartWithNewFileName) - _restartWithNewFileName=false; - string suppress ="rm -rf "+resultFile+"_*"; - system(suppress.c_str());//Nettoyage des précédents calculs identiques - - _VV.setInfoOnComponent(0,"Temperature_(K)"); - switch(_saveFormat) - { - case VTK : - _VV.writeVTK(resultFile); - break; - case MED : - _VV.writeMED(resultFile); - break; - case CSV : - _VV.writeCSV(resultFile); - break; + + if(_mpi_rank==0){ + //On remplit le champ + double Ti; + if(!_FECalculation) + for(int i =0; i<_Nmailles;i++) + { + if(_mpi_size>1) + VecGetValues(_Tn_seq, 1, &i, &Ti); + else + VecGetValues(_Tn , 1, &i, &Ti); + + _VV(i)=Ti; + } + else + { + int globalIndex; + for(int i=0; i<_NunknownNodes; i++) + { + if(_mpi_size>1) + VecGetValues(_Tn_seq, 1, &i, &Ti); + else + VecGetValues(_Tk , 1, &i, &Ti); + globalIndex = globalNodeIndex(i, _dirichletNodeIds); + _VV(globalIndex)=Ti; + } + + Node Ni; + string nameOfGroup; + for(int i=0; i<_NdirichletNodes; i++)//Assumes node numbering starts with border nodes + { + Ni=_mesh.getNode(_dirichletNodeIds[i]); + nameOfGroup = Ni.getGroupName(); + _VV(_dirichletNodeIds[i])=_limitField[nameOfGroup].T; + } + } + _VV.setTime(_time,_nbTimeStep); + + // create mesh and component info + if (_nbTimeStep ==0 || _restartWithNewFileName){ + if (_restartWithNewFileName) + _restartWithNewFileName=false; + string suppress ="rm -rf "+resultFile+"_*"; + system(suppress.c_str());//Nettoyage des précédents calculs identiques + + _VV.setInfoOnComponent(0,"Temperature_(K)"); + switch(_saveFormat) + { + case VTK : + _VV.writeVTK(resultFile); + break; + case MED : + _VV.writeMED(resultFile); + break; + case CSV : + _VV.writeCSV(resultFile); + break; + } } - } - else{ // do not create mesh - switch(_saveFormat) - { - case VTK : - _VV.writeVTK(resultFile,false); - break; - case MED : - _VV.writeMED(resultFile,false); - break; - case CSV : - _VV.writeCSV(resultFile); - break; + else{ // do not create mesh + switch(_saveFormat) + { + case VTK : + _VV.writeVTK(resultFile,false); + break; + case MED : + _VV.writeMED(resultFile,false); + break; + case CSV : + _VV.writeCSV(resultFile); + break; + } } + + if(_isStationary) + { + resultFile+="_Stat"; + switch(_saveFormat) + { + case VTK : + _VV.writeVTK(resultFile); + break; + case MED : + _VV.writeMED(resultFile); + break; + case CSV : + _VV.writeCSV(resultFile); + break; + } + } } - - if(_isStationary) - { - resultFile+="_Stat"; - switch(_saveFormat) - { - case VTK : - _VV.writeVTK(resultFile); - break; - case MED : - _VV.writeMED(resultFile); - break; - case CSV : - _VV.writeCSV(resultFile); - break; - } - } } void DiffusionEquation::terminate(){ @@ -808,6 +823,8 @@ void DiffusionEquation::terminate(){ VecDestroy(&_b0); VecDestroy(&_b); MatDestroy(&_A); + if(_mpi_size>1 && _mpi_rank == 0) + VecDestroy(&_Tn_seq); } void @@ -825,21 +842,45 @@ DiffusionEquation::setNeumannValues(map< int, double> neumannBoundaryValues) } -vector DiffusionEquation::getOutputFieldsNames() +Field& +DiffusionEquation::getOutputTemperatureField() +{ + if(!_initializedMemory) + throw("Computation not initialized. No temperature field available"); + else + return _VV; +} + +Field& +DiffusionEquation::getRodTemperatureField() +{ + return getOutputTemperatureField(); +} + +vector +DiffusionEquation::getInputFieldsNames() { vector result(2); result[0]="FluidTemperature"; - result[1]="RodTemperature"; + result[1]="HeatPower"; + + return result; +} +vector +DiffusionEquation::getOutputFieldsNames() +{ + vector result(1); + + result[0]="RodTemperature"; return result; } -Field& DiffusionEquation::getOutputField(const string& nameField ) +Field& +DiffusionEquation::getOutputField(const string& nameField ) { - if(nameField=="FluidTemperature" || nameField=="FLUIDTEMPERATURE" || nameField=="TemperatureFluide" || nameField=="TEMPERATUREFLUIDE" ) - return getFluidTemperatureField(); - else if(nameField=="RodTemperature" || nameField=="RODTEMPERATURE" || nameField=="TEMPERATURECOMBUSTIBLE" || nameField=="TemperatureCombustible" ) + if(nameField=="RodTemperature" || nameField=="RODTEMPERATURE" || nameField=="TEMPERATURECOMBUSTIBLE" || nameField=="TemperatureCombustible" ) return getRodTemperatureField(); else { @@ -848,3 +889,16 @@ Field& DiffusionEquation::getOutputField(const string& nameField ) } } +void +DiffusionEquation::setInputField(const string& nameField, Field& inputField ) +{ + if(nameField=="FluidTemperature" || nameField=="FLUIDTEMPERATURE" || nameField=="TemperatureFluide" || nameField=="TEMPERATUREFLUIDE") + return setFluidTemperatureField( inputField) ; + else if(nameField=="HeatPower" || nameField=="HEATPOWER" || nameField=="PuissanceThermique" || nameField=="PUISSANCETHERMIQUE" ) + return setHeatPowerField( inputField ); + else + { + cout<<"Error : Field name "<< nameField << " is not an input field name, call getInputFieldsNames first" << endl; + throw CdmathException("DiffusionEquation::setInputField error : Unknown Field name"); + } +} diff --git a/CoreFlows/Models/src/ProblemCoreFlows.cxx b/CoreFlows/Models/src/ProblemCoreFlows.cxx index 05ec96b..b53cb75 100755 --- a/CoreFlows/Models/src/ProblemCoreFlows.cxx +++ b/CoreFlows/Models/src/ProblemCoreFlows.cxx @@ -32,18 +32,18 @@ ProblemCoreFlows::ProblemCoreFlows(MPI_Comm comm) PETSC_COMM_WORLD = comm; PetscInitialize(NULL,NULL,0,0);//Note this is ok if MPI has been been initialised independently from PETSC } - MPI_Comm_rank(PETSC_COMM_WORLD,&_rank); - MPI_Comm_size(PETSC_COMM_WORLD,&_size); - PetscPrintf(PETSC_COMM_WORLD,"Simulation on %d processors\n",_size);//Prints to standard out, only from the first processor in the communicator. Calls from other processes are ignored. - PetscSynchronizedPrintf(PETSC_COMM_WORLD,"Processor [%d] ready for action\n",_rank);//Prints synchronized output from several processors. Output of the first processor is followed by that of the second, etc. + MPI_Comm_rank(PETSC_COMM_WORLD,&_mpi_rank); + MPI_Comm_size(PETSC_COMM_WORLD,&_mpi_size); + PetscPrintf(PETSC_COMM_WORLD,"Simulation on %d processors\n",_mpi_size);//Prints to standard out, only from the first processor in the communicator. Calls from other processes are ignored. + PetscSynchronizedPrintf(PETSC_COMM_WORLD,"Processor [%d] ready for action\n",_mpi_rank);//Prints synchronized output from several processors. Output of the first processor is followed by that of the second, etc. PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT); - if(_size>1) + if(_mpi_size>1) { PetscPrintf(PETSC_COMM_WORLD,"---- More than one processor detected : running a parallel simulation ----\n"); PetscPrintf(PETSC_COMM_WORLD,"---- Limited parallelism : input and output fields remain sequential ----\n"); - PetscPrintf(PETSC_COMM_WORLD,"---- Only the matrixoperations are done in parallel thanks to PETSc----\n"); - PetscPrintf(PETSC_COMM_WORLD,"---- Processor %d is in charge of building the mesh, saving the results, filling and then distributing the matrix to other processors.\n\n",_rank); + PetscPrintf(PETSC_COMM_WORLD,"---- Only the matrix operations are done in parallel thanks to PETSc----\n"); + PetscPrintf(PETSC_COMM_WORLD,"---- Processor %d is in charge of building the mesh, saving the results, filling and then distributing the matrix to other processors.\n\n",_mpi_rank); } /* Numerical parameter */ @@ -57,6 +57,7 @@ ProblemCoreFlows::ProblemCoreFlows(MPI_Comm comm) _precision_Newton=_precision; _erreur_rel= 0; _isStationary=false; + _stationaryMode=false; _timeScheme=Explicit; _wellBalancedCorrection=false; _FECalculation=false; @@ -76,7 +77,8 @@ ProblemCoreFlows::ProblemCoreFlows(MPI_Comm comm) /* Mesh parameters */ _Ndim=0; _minl=0; - _neibMaxNb=0; + _neibMaxNbCells=0; + _neibMaxNbNodes=0; /* Memory and restart */ _initialDataSet=false; @@ -89,12 +91,14 @@ ProblemCoreFlows::ProblemCoreFlows(MPI_Comm comm) _maxNewtonIts=50; _PetscIts=0;//the number of iterations of the linear solver _ksptype = (char*)&KSPGMRES; - _pctype = (char*)&PCLU; + if( _mpi_size>1) + _pctype = (char*)&PCNONE; + else + _pctype = (char*)&PCLU; - /* Physical parameter */ + /* Physical parameters */ _heatPowerFieldSet=false; _heatTransfertCoeff=0; - _rodTemperatureFieldSet=false; _heatSource=0; //extracting current directory @@ -128,7 +132,7 @@ double ProblemCoreFlows::presentTime() const void ProblemCoreFlows::setTimeMax(double timeMax){ _timeMax = timeMax; } -void ProblemCoreFlows::setPresentTime (double time) +void ProblemCoreFlows::resetTime (double time) { _time=time; } @@ -145,7 +149,6 @@ void ProblemCoreFlows::setPrecision(double precision) } void ProblemCoreFlows::setInitialField(const Field &VV) { - if(_Ndim != VV.getSpaceDimension()){ *_runLogFile<<"ProblemCoreFlows::setInitialField: mesh has incorrect space dimension"<close(); @@ -174,6 +177,7 @@ void ProblemCoreFlows::setInitialField(const Field &VV) _VV.setName("SOLVERLAB results"); _time=_VV.getTime(); _mesh=_VV.getMesh(); + _initialDataSet=true; //Mesh data @@ -189,33 +193,53 @@ void ProblemCoreFlows::setInitialField(const Field &VV) Face Fk; if(_verbose) - PetscPrintf(PETSC_COMM_WORLD,"Computing cell perimeters and mesh minimal diameter\n"); + PetscPrintf(PETSC_COMM_SELF,"Processor %d Computing cell perimeters and mesh minimal diameter\n", _mpi_rank); //Compute the maximum number of neighbours for nodes or cells if(VV.getTypeOfField()==NODES) _neibMaxNbNodes=_mesh.getMaxNbNeighbours(NODES); else - _neibMaxNb=_mesh.getMaxNbNeighbours(CELLS); + { + _neibMaxNbCells=_mesh.getMaxNbNeighbours(CELLS); - //Compute Delta x and the cell perimeters - for (int i=0; i<_mesh.getNumberOfCells(); i++){ - Ci = _mesh.getCell(i); - if (_Ndim > 1){ - _perimeters(i)=0; - for (int k=0 ; k 1){ + _perimeters(i)=0; + for (int k=0 ; kopen((_fileName+".log").c_str(), ios::out | ios::trunc);;//for creation of a log file to save the history of the simulation *_runLogFile<< "Running test case "<< _fileName<2 && _Ndim==1 ){//inner face with more than two neighbours if(_verbose && (_nbTimeStep-1)%_freqSave ==0) - cout<<"lattice mesh junction at face "< +ProblemFluid::getInputFieldsNames() +{ + vector result(1); + + result[0]="HeatPower"; + result[1]="Porosity"; + result[2]="PressureLoss"; + result[3]="Section"; + + return result; +} + +void +ProblemFluid::setInputField(const string& nameField, Field& inputField ) +{ + if(nameField=="Porosity" || nameField=="POROSITY" || nameField=="Porosité" || nameField=="POROSITE") + return setPorosityField( inputField) ; + else if(nameField=="HeatPower" || nameField=="HEATPOWER" || nameField=="PuissanceThermique" || nameField=="PUISSANCETHERMIQUE" ) + return setHeatPowerField( inputField ); + else if(nameField=="PressureLoss" || nameField=="PRESSURELOSS" || nameField=="PerteDeCharge" || nameField=="PPERTEDECHARGE") + return setPressureLossField( inputField) ; + else if(nameField=="Section" || nameField=="SECTION" ) + return setSectionField( inputField ); + else + { + cout<<"Error : Field name "<< nameField << " is not an input field name, call getInputFieldsNames first" << endl; + throw CdmathException("SinglePhase::setInputField error : Unknown Field name"); + } +} diff --git a/CoreFlows/Models/src/StationaryDiffusionEquation.cxx b/CoreFlows/Models/src/StationaryDiffusionEquation.cxx index 9e6fdc0..f813ad5 100755 --- a/CoreFlows/Models/src/StationaryDiffusionEquation.cxx +++ b/CoreFlows/Models/src/StationaryDiffusionEquation.cxx @@ -22,10 +22,10 @@ StationaryDiffusionEquation::StationaryDiffusionEquation(int dim, bool FECalcula PETSC_COMM_WORLD = comm; PetscInitialize(NULL,NULL,0,0);//Note this is ok if MPI has been been initialised independently from PETSC } - MPI_Comm_rank(PETSC_COMM_WORLD,&_rank); - MPI_Comm_size(PETSC_COMM_WORLD,&_size); - PetscPrintf(PETSC_COMM_WORLD,"Simulation on %d processors\n",_size);//Prints to standard out, only from the first processor in the communicator. Calls from other processes are ignored. - PetscSynchronizedPrintf(PETSC_COMM_WORLD,"Processor [%d] ready for action\n",_rank);//Prints synchronized output from several processors. Output of the first processor is followed by that of the second, etc. + MPI_Comm_rank(PETSC_COMM_WORLD,&_mpi_rank); + MPI_Comm_size(PETSC_COMM_WORLD,&_mpi_size); + PetscPrintf(PETSC_COMM_WORLD,"Simulation on %d processors\n",_mpi_size);//Prints to standard out, only from the first processor in the communicator. Calls from other processes are ignored. + PetscSynchronizedPrintf(PETSC_COMM_WORLD,"Processor [%d] ready for action\n",_mpi_rank);//Prints synchronized output from several processors. Output of the first processor is followed by that of the second, etc. PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT); if(lambda < 0.) @@ -71,7 +71,11 @@ StationaryDiffusionEquation::StationaryDiffusionEquation(int dim, bool FECalcula _NEWTON_its=0; int _PetscIts=0;//the number of iterations of the linear solver _ksptype = (char*)&KSPGMRES; - _pctype = (char*)&PCLU; + if( _mpi_size>1) + _pctype = (char*)&PCNONE; + else + _pctype = (char*)&PCLU; + _conditionNumber=false; _erreur_rel= 0; @@ -104,106 +108,118 @@ StationaryDiffusionEquation::StationaryDiffusionEquation(int dim, bool FECalcula void StationaryDiffusionEquation::initialize() { - _runLogFile->open((_fileName+".log").c_str(), ios::out | ios::trunc);;//for creation of a log file to save the history of the simulation - - if(!_meshSet) - throw CdmathException("StationaryDiffusionEquation::initialize() set mesh first"); - else - { - cout<<"!!!! Initialisation of the computation of the temperature diffusion in a solid using "; - *_runLogFile<<"!!!!! Initialisation of the computation of the temperature diffusion in a solid using "; - if(!_FECalculation) - { - cout<< "Finite volumes method"<open((_fileName+".log").c_str(), ios::out | ios::trunc);;//for creation of a log file to save the history of the simulation + + if(!_meshSet) + throw CdmathException("StationaryDiffusionEquation::initialize() set mesh first"); + else + { + cout<<"!!!! Initialisation of the computation of the temperature diffusion in a solid using "; + *_runLogFile<<"!!!!! Initialisation of the computation of the temperature diffusion in a solid using "; + if(!_FECalculation) + { + cout<< "Finite volumes method"<::iterator it=_dirichletBoundaryValues.find(_boundaryNodeIds[i]); + if( it != _dirichletBoundaryValues.end() ) + _dirichletNodeIds.push_back(_boundaryNodeIds[i]); + else if( _mesh.getNode(_boundaryNodeIds[i]).getGroupNames().size()==0 ) + { + cout<<"!!! No boundary group set for boundary node" << _boundaryNodeIds[i]<< endl; + *_runLogFile<< "!!! No boundary group set for boundary node" << _boundaryNodeIds[i]<close(); + throw CdmathException("Missing boundary group"); + } + else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType==NoneBCStationaryDiffusion) + { + cout<<"!!! No boundary condition set for boundary node " << _boundaryNodeIds[i]<< endl; + cout<<"!!! Accepted boundary conditions are DirichletStationaryDiffusion "<< DirichletStationaryDiffusion <<" and NeumannStationaryDiffusion "<< NeumannStationaryDiffusion << endl; + *_runLogFile<< "!!! No boundary condition set for boundary node " << _boundaryNodeIds[i]<close(); + throw CdmathException("Missing boundary condition"); + } + else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType==DirichletStationaryDiffusion) + _dirichletNodeIds.push_back(_boundaryNodeIds[i]); + else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType!=NeumannStationaryDiffusion) + { + cout<<"!!! Wrong boundary condition "<< _limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType<< " set for boundary node " << _boundaryNodeIds[i]<< endl; + cout<<"!!! Accepted boundary conditions are DirichletStationaryDiffusion "<< DirichletStationaryDiffusion <<" and NeumannStationaryDiffusion "<< NeumannStationaryDiffusion << endl; + *_runLogFile<< "!!! Wrong boundary condition "<< _limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType<< " set for boundary node " << _boundaryNodeIds[i]<close(); + throw CdmathException("Wrong boundary condition"); + } + } + _NdirichletNodes=_dirichletNodeIds.size(); + _NunknownNodes=_Nnodes - _NdirichletNodes; + cout<<"Number of unknown nodes " << _NunknownNodes <<", Number of boundary nodes " << _NboundaryNodes<< ", Number of Dirichlet boundary nodes " << _NdirichletNodes <::iterator it=_dirichletBoundaryValues.find(_boundaryNodeIds[i]); - if( it != _dirichletBoundaryValues.end() ) - _dirichletNodeIds.push_back(_boundaryNodeIds[i]); - else if( _mesh.getNode(_boundaryNodeIds[i]).getGroupNames().size()==0 ) - { - cout<<"!!! No boundary group set for boundary node" << _boundaryNodeIds[i]<< endl; - *_runLogFile<< "!!! No boundary group set for boundary node" << _boundaryNodeIds[i]<close(); - throw CdmathException("Missing boundary group"); - } - else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType==NoneBCStationaryDiffusion) - { - cout<<"!!! No boundary condition set for boundary node " << _boundaryNodeIds[i]<< endl; - cout<<"!!! Accepted boundary conditions are DirichletStationaryDiffusion "<< DirichletStationaryDiffusion <<" and NeumannStationaryDiffusion "<< NeumannStationaryDiffusion << endl; - *_runLogFile<< "!!! No boundary condition set for boundary node " << _boundaryNodeIds[i]<close(); - throw CdmathException("Missing boundary condition"); - } - else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType==DirichletStationaryDiffusion) - _dirichletNodeIds.push_back(_boundaryNodeIds[i]); - else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType!=NeumannStationaryDiffusion) - { - cout<<"!!! Wrong boundary condition "<< _limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType<< " set for boundary node " << _boundaryNodeIds[i]<< endl; - cout<<"!!! Accepted boundary conditions are DirichletStationaryDiffusion "<< DirichletStationaryDiffusion <<" and NeumannStationaryDiffusion "<< NeumannStationaryDiffusion << endl; - *_runLogFile<< "!!! Wrong boundary condition "<< _limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType<< " set for boundary node " << _boundaryNodeIds[i]<close(); - throw CdmathException("Wrong boundary condition"); - } - } - _NdirichletNodes=_dirichletNodeIds.size(); - _NunknownNodes=_Nnodes - _NdirichletNodes; - cout<<"Number of unknown nodes " << _NunknownNodes <<", Number of boundary nodes " << _NboundaryNodes<< ", Number of Dirichlet boundary nodes " << _NdirichletNodes <1 && _mpi_rank == 0) + VecCreateSeq(PETSC_COMM_SELF,_globalNbUnknowns,&_Tk_seq);//For saving results on proc 0 + if(_mpi_size==0) + _Tk_seq=_Tk; + VecScatterCreateToZero(_Tk,&_scat,&_Tk_seq); //Linear solver - KSPCreate(PETSC_COMM_SELF, &_ksp); + KSPCreate(PETSC_COMM_WORLD, &_ksp); KSPSetType(_ksp, _ksptype); // if(_ksptype == KSPGMRES) KSPGMRESSetRestart(_ksp,10000); KSPSetTolerances(_ksp,_precision,_precision,PETSC_DEFAULT,_maxPetscIts); @@ -259,124 +275,126 @@ double StationaryDiffusionEquation::computeDiffusionMatrix(bool & stop) { double result; + MatZeroEntries(_A); + VecZeroEntries(_b); + if(_FECalculation) result=computeDiffusionMatrixFE(stop); else result=computeDiffusionMatrixFV(stop); + MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY); + MatAssemblyEnd( _A, MAT_FINAL_ASSEMBLY); + //Contribution from the solid/fluid heat exchange with assumption of constant heat transfer coefficient //update value here if variable heat transfer coefficient if(_heatTransfertCoeff>_precision) MatShift(_A,_heatTransfertCoeff);//Contribution from the liquit/solid heat transfer if(_verbose or _system) - MatView(_A,PETSC_VIEWER_STDOUT_SELF); + MatView(_A,PETSC_VIEWER_STDOUT_WORLD); return result; } double StationaryDiffusionEquation::computeDiffusionMatrixFE(bool & stop){ - Cell Cj; - string nameOfGroup; - double coeff; - MatZeroEntries(_A); - VecZeroEntries(_b); - - Matrix M(_Ndim+1,_Ndim+1);//cell geometry matrix - std::vector< Vector > GradShapeFuncs(_Ndim+1);//shape functions of cell nodes - std::vector< int > nodeIds(_Ndim+1);//cell node Ids - std::vector< Node > nodes(_Ndim+1);//cell nodes - int i_int, j_int; //index of nodes j and k considered as unknown nodes - bool dirichletCell_treated; - - std::vector< vector< double > > values(_Ndim+1,vector< double >(_Ndim+1,0));//values of shape functions on cell node - for (int idim=0; idim<_Ndim+1;idim++) - values[idim][idim]=1; - - /* parameters for boundary treatment */ - vector< double > valuesBorder(_Ndim+1); - Vector GradShapeFuncBorder(_Ndim+1); - - for (int j=0; j<_Nmailles;j++) - { - Cj = _mesh.getCell(j); - - for (int idim=0; idim<_Ndim+1;idim++){ - nodeIds[idim]=Cj.getNodeId(idim); - nodes[idim]=_mesh.getNode(nodeIds[idim]); - for (int jdim=0; jdim<_Ndim;jdim++) - M(idim,jdim)=nodes[idim].getPoint()[jdim]; - M(idim,_Ndim)=1; - } - for (int idim=0; idim<_Ndim+1;idim++) - GradShapeFuncs[idim]=DiffusionEquation::gradientNodal(M,values[idim])/DiffusionEquation::fact(_Ndim); - - /* Loop on the edges of the cell */ - for (int idim=0; idim<_Ndim+1;idim++) - { - if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),nodeIds[idim])==_dirichletNodeIds.end())//!_mesh.isBorderNode(nodeIds[idim]) - {//First node of the edge is not Dirichlet node - i_int=DiffusionEquation::unknownNodeIndex(nodeIds[idim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing - dirichletCell_treated=false; - for (int jdim=0; jdim<_Ndim+1;jdim++) - { - if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),nodeIds[jdim])==_dirichletNodeIds.end())//!_mesh.isBorderNode(nodeIds[jdim]) - {//Second node of the edge is not Dirichlet node - j_int= DiffusionEquation::unknownNodeIndex(nodeIds[jdim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing - MatSetValue(_A,i_int,j_int,(_DiffusionTensor*GradShapeFuncs[idim])*GradShapeFuncs[jdim]/Cj.getMeasure(), ADD_VALUES); - } - else if (!dirichletCell_treated) - {//Second node of the edge is a Dirichlet node - dirichletCell_treated=true; - for (int kdim=0; kdim<_Ndim+1;kdim++) - { - std::map::iterator it=_dirichletBoundaryValues.find(nodeIds[kdim]); - if( it != _dirichletBoundaryValues.end() ) - { - if( _dirichletValuesSet ) - valuesBorder[kdim]=_dirichletBoundaryValues[it->second]; - else - valuesBorder[kdim]=_limitField[_mesh.getNode(nodeIds[kdim]).getGroupName()].T; - } - else - valuesBorder[kdim]=0; - } - GradShapeFuncBorder=DiffusionEquation::gradientNodal(M,valuesBorder)/DiffusionEquation::fact(_Ndim); - coeff =-1.*(_DiffusionTensor*GradShapeFuncBorder)*GradShapeFuncs[idim]/Cj.getMeasure(); - VecSetValue(_b,i_int,coeff, ADD_VALUES); - } - } - } - } + if(_mpi_rank == 0) + { + Cell Cj; + string nameOfGroup; + double coeff; + + Matrix M(_Ndim+1,_Ndim+1);//cell geometry matrix + std::vector< Vector > GradShapeFuncs(_Ndim+1);//shape functions of cell nodes + std::vector< int > nodeIds(_Ndim+1);//cell node Ids + std::vector< Node > nodes(_Ndim+1);//cell nodes + int i_int, j_int; //index of nodes j and k considered as unknown nodes + bool dirichletCell_treated; + + std::vector< vector< double > > values(_Ndim+1,vector< double >(_Ndim+1,0));//values of shape functions on cell node + for (int idim=0; idim<_Ndim+1;idim++) + values[idim][idim]=1; + + /* parameters for boundary treatment */ + vector< double > valuesBorder(_Ndim+1); + Vector GradShapeFuncBorder(_Ndim+1); + + for (int j=0; j<_Nmailles;j++) + { + Cj = _mesh.getCell(j); + + for (int idim=0; idim<_Ndim+1;idim++){ + nodeIds[idim]=Cj.getNodeId(idim); + nodes[idim]=_mesh.getNode(nodeIds[idim]); + for (int jdim=0; jdim<_Ndim;jdim++) + M(idim,jdim)=nodes[idim].getPoint()[jdim]; + M(idim,_Ndim)=1; + } + for (int idim=0; idim<_Ndim+1;idim++) + GradShapeFuncs[idim]=DiffusionEquation::gradientNodal(M,values[idim])/DiffusionEquation::fact(_Ndim); + + /* Loop on the edges of the cell */ + for (int idim=0; idim<_Ndim+1;idim++) + { + if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),nodeIds[idim])==_dirichletNodeIds.end())//!_mesh.isBorderNode(nodeIds[idim]) + {//First node of the edge is not Dirichlet node + i_int=DiffusionEquation::unknownNodeIndex(nodeIds[idim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing + dirichletCell_treated=false; + for (int jdim=0; jdim<_Ndim+1;jdim++) + { + if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),nodeIds[jdim])==_dirichletNodeIds.end())//!_mesh.isBorderNode(nodeIds[jdim]) + {//Second node of the edge is not Dirichlet node + j_int= DiffusionEquation::unknownNodeIndex(nodeIds[jdim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing + MatSetValue(_A,i_int,j_int,(_DiffusionTensor*GradShapeFuncs[idim])*GradShapeFuncs[jdim]/Cj.getMeasure(), ADD_VALUES); + } + else if (!dirichletCell_treated) + {//Second node of the edge is a Dirichlet node + dirichletCell_treated=true; + for (int kdim=0; kdim<_Ndim+1;kdim++) + { + std::map::iterator it=_dirichletBoundaryValues.find(nodeIds[kdim]); + if( it != _dirichletBoundaryValues.end() ) + { + if( _dirichletValuesSet ) + valuesBorder[kdim]=_dirichletBoundaryValues[it->second]; + else + valuesBorder[kdim]=_limitField[_mesh.getNode(nodeIds[kdim]).getGroupName()].T; + } + else + valuesBorder[kdim]=0; + } + GradShapeFuncBorder=DiffusionEquation::gradientNodal(M,valuesBorder)/DiffusionEquation::fact(_Ndim); + coeff =-1.*(_DiffusionTensor*GradShapeFuncBorder)*GradShapeFuncs[idim]/Cj.getMeasure(); + VecSetValue(_b,i_int,coeff, ADD_VALUES); + } + } + } + } + } + + //Calcul de la contribution de la condition limite de Neumann au second membre + if( _NdirichletNodes !=_NboundaryNodes) + { + vector< int > boundaryFaces = _mesh.getBoundaryFaceIds(); + int NboundaryFaces=boundaryFaces.size(); + for(int i = 0; i< NboundaryFaces ; i++)//On parcourt les faces du bord + { + Face Fi = _mesh.getFace(boundaryFaces[i]); + for(int j = 0 ; j<_Ndim ; j++)//On parcourt les noeuds de la face + { + if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),Fi.getNodeId(j))==_dirichletNodeIds.end())//node j is a Neumann BC node (not a Dirichlet BC node) + { + j_int=DiffusionEquation::unknownNodeIndex(Fi.getNodeId(j), _dirichletNodeIds);//indice du noeud j en tant que noeud inconnu + if( _neumannValuesSet ) + coeff =Fi.getMeasure()/_Ndim*_neumannBoundaryValues[Fi.getNodeId(j)]; + else + coeff =Fi.getMeasure()/_Ndim*_limitField[_mesh.getNode(Fi.getNodeId(j)).getGroupName()].normalFlux; + VecSetValue(_b, j_int, coeff, ADD_VALUES); + } + } + } + } } - - //Calcul de la contribution de la condition limite de Neumann au second membre - if( _NdirichletNodes !=_NboundaryNodes) - { - vector< int > boundaryFaces = _mesh.getBoundaryFaceIds(); - int NboundaryFaces=boundaryFaces.size(); - for(int i = 0; i< NboundaryFaces ; i++)//On parcourt les faces du bord - { - Face Fi = _mesh.getFace(boundaryFaces[i]); - for(int j = 0 ; j<_Ndim ; j++)//On parcourt les noeuds de la face - { - if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),Fi.getNodeId(j))==_dirichletNodeIds.end())//node j is a Neumann BC node (not a Dirichlet BC node) - { - j_int=DiffusionEquation::unknownNodeIndex(Fi.getNodeId(j), _dirichletNodeIds);//indice du noeud j en tant que noeud inconnu - if( _neumannValuesSet ) - coeff =Fi.getMeasure()/_Ndim*_neumannBoundaryValues[Fi.getNodeId(j)]; - else - coeff =Fi.getMeasure()/_Ndim*_limitField[_mesh.getNode(Fi.getNodeId(j)).getGroupName()].normalFlux; - VecSetValue(_b, j_int, coeff, ADD_VALUES); - } - } - } - } - - MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY); - MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY); - VecAssemblyBegin(_b); - VecAssemblyEnd(_b); if(_onlyNeumannBC) //Check that the matrix is symmetric { @@ -396,114 +414,110 @@ double StationaryDiffusionEquation::computeDiffusionMatrixFE(bool & stop){ } double StationaryDiffusionEquation::computeDiffusionMatrixFV(bool & stop){ - long nbFaces = _mesh.getNumberOfFaces(); - Face Fj; - Cell Cell1,Cell2; - string nameOfGroup; - double inv_dxi, inv_dxj; - double barycenterDistance; - Vector normale(_Ndim); - double dn; - PetscInt idm, idn; - std::vector< int > idCells; - MatZeroEntries(_A); - VecZeroEntries(_b); - - for (int j=0; j::iterator it=_dirichletBoundaryValues.find(j); - if( it != _dirichletBoundaryValues.end() ) - { - barycenterDistance=Cell1.getBarryCenter().distance(Fj.getBarryCenter()); - MatSetValue(_A,idm,idm,dn*inv_dxi/barycenterDistance , ADD_VALUES); - VecSetValue(_b,idm, dn*inv_dxi/barycenterDistance*it->second, ADD_VALUES); - } - else - { - nameOfGroup = Fj.getGroupName(); - - if (_limitField[nameOfGroup].bcType==NeumannStationaryDiffusion){ - VecSetValue(_b,idm, -dn*inv_dxi*_limitField[nameOfGroup].normalFlux, ADD_VALUES); - } - else if(_limitField[nameOfGroup].bcType==DirichletStationaryDiffusion){ - barycenterDistance=Cell1.getBarryCenter().distance(Fj.getBarryCenter()); - MatSetValue(_A,idm,idm,dn*inv_dxi/barycenterDistance , ADD_VALUES); - VecSetValue(_b,idm, dn*inv_dxi/barycenterDistance*_limitField[nameOfGroup].T, ADD_VALUES); - } - else { - stop=true ; - cout<<"!!!!!!!!!!!!!!! Error StationaryDiffusionEquation::computeDiffusionMatrixFV !!!!!!!!!!"< idCells; + + for (int j=0; j::iterator it=_dirichletBoundaryValues.find(j); + if( it != _dirichletBoundaryValues.end() ) + { + barycenterDistance=Cell1.getBarryCenter().distance(Fj.getBarryCenter()); + MatSetValue(_A,idm,idm,dn*inv_dxi/barycenterDistance , ADD_VALUES); + VecSetValue(_b,idm, dn*inv_dxi/barycenterDistance*it->second, ADD_VALUES); + } + else + { + nameOfGroup = Fj.getGroupName(); + + if (_limitField[nameOfGroup].bcType==NeumannStationaryDiffusion){ + VecSetValue(_b,idm, -dn*inv_dxi*_limitField[nameOfGroup].normalFlux, ADD_VALUES); + } + else if(_limitField[nameOfGroup].bcType==DirichletStationaryDiffusion){ + barycenterDistance=Cell1.getBarryCenter().distance(Fj.getBarryCenter()); + MatSetValue(_A,idm,idm,dn*inv_dxi/barycenterDistance , ADD_VALUES); + VecSetValue(_b,idm, dn*inv_dxi/barycenterDistance*_limitField[nameOfGroup].T, ADD_VALUES); + } + else { + stop=true ; + cout<<"!!!!!!!!!!!!!!! Error StationaryDiffusionEquation::computeDiffusionMatrixFV !!!!!!!!!!"< 1) + inv_dxj = Fj.getMeasure()/Cell2.getMeasure(); + else + inv_dxj = 1/Cell2.getMeasure(); + + barycenterDistance=Cell1.getBarryCenter().distance(Cell2.getBarryCenter()); + + MatSetValue(_A,idm,idm, dn*inv_dxi/barycenterDistance, ADD_VALUES); + MatSetValue(_A,idm,idn,-dn*inv_dxi/barycenterDistance, ADD_VALUES); + MatSetValue(_A,idn,idn, dn*inv_dxj/barycenterDistance, ADD_VALUES); + MatSetValue(_A,idn,idm,-dn*inv_dxj/barycenterDistance, ADD_VALUES); } - Cell2 = _mesh.getCell(idCells[1]); - idn = idCells[1]; - if (_Ndim > 1) - inv_dxj = Fj.getMeasure()/Cell2.getMeasure(); else - inv_dxj = 1/Cell2.getMeasure(); - - barycenterDistance=Cell1.getBarryCenter().distance(Cell2.getBarryCenter()); - - MatSetValue(_A,idm,idm, dn*inv_dxi/barycenterDistance, ADD_VALUES); - MatSetValue(_A,idm,idn,-dn*inv_dxi/barycenterDistance, ADD_VALUES); - MatSetValue(_A,idn,idn, dn*inv_dxj/barycenterDistance, ADD_VALUES); - MatSetValue(_A,idn,idm,-dn*inv_dxj/barycenterDistance, ADD_VALUES); + { + *_runLogFile<<"StationaryDiffusionEquation::computeDiffusionMatrixFV(): incompatible number of cells around a face"< nodesId; - for (int i=0; i<_Nmailles;i++) - { - Ci=_mesh.getCell(i); - nodesId=Ci.getNodesId(); - for (int j=0; j nodesId; + for (int i=0; i<_Nmailles;i++) + { + Ci=_mesh.getCell(i); + nodesId=Ci.getNodesId(); + for (int j=0; jclose(); - throw CdmathException("StationaryDiffusionEquation::setMesh: mesh has incorrect dimension"); + if(_Ndim != M.getSpaceDimension() or _Ndim!=M.getMeshDimension())//for the moment we must have space dim=mesh dim + { + cout<< "Problem : dimension defined is "<<_Ndim<< " but mesh dimension= "<close(); + throw CdmathException("StationaryDiffusionEquation::setMesh: mesh has incorrect dimension"); + } + + _mesh=M; + _Nmailles = _mesh.getNumberOfCells(); + _Nnodes = _mesh.getNumberOfNodes(); + + cout<<"Mesh has "<< _Nmailles << " cells and " << _Nnodes << " nodes"<close(); + throw CdmathException("StationaryDiffusionEquation::setMesh: mesh has incorrect cell types"); + } + } + else if(_Ndim==3) + { + if( _mesh.isTetrahedral() )//Mesh dim=3 + cout<<"3D Finite element method on tetrahedra"<close(); + throw CdmathException("StationaryDiffusionEquation::setMesh: mesh has incorrect cell types"); + } + } + + _VV=Field ("Temperature", NODES, _mesh, 1); + + _neibMaxNbNodes=_mesh.getMaxNbNeighbours(NODES); + _boundaryNodeIds = _mesh.getBoundaryNodeIds(); + _NboundaryNodes=_boundaryNodeIds.size(); + } } - - _mesh=M; - _Nmailles = _mesh.getNumberOfCells(); - _Nnodes = _mesh.getNumberOfNodes(); - - cout<<"Mesh has "<< _Nmailles << " cells and " << _Nnodes << " nodes"<close(); - throw CdmathException("StationaryDiffusionEquation::setMesh: mesh has incorrect cell types"); - } - } - else if(_Ndim==3) - { - if( _mesh.isTetrahedral() )//Mesh dim=3 - cout<<"3D Finite element method on tetrahedra"<close(); - throw CdmathException("StationaryDiffusionEquation::setMesh: mesh has incorrect cell types"); - } - } - - _VV=Field ("Temperature", NODES, _mesh, 1); - - _neibMaxNbNodes=_mesh.getMaxNbNeighbours(NODES); - _boundaryNodeIds = _mesh.getBoundaryNodeIds(); - _NboundaryNodes=_boundaryNodeIds.size(); - } + + /* MPI distribution parameters */ + int nbVoisinsMax;//Mettre en attribut ? + if(!_FECalculation){ + MPI_Bcast(&_Nmailles , 1, MPI_INT, 0, PETSC_COMM_WORLD); + MPI_Bcast(&_neibMaxNbCells, 1, MPI_INT, 0, PETSC_COMM_WORLD); + nbVoisinsMax = _neibMaxNbCells; + } + else{ + MPI_Bcast(&_Nnodes , 1, MPI_INT, 0, PETSC_COMM_WORLD); + MPI_Bcast(&_neibMaxNbNodes, 1, MPI_INT, 0, PETSC_COMM_WORLD); + nbVoisinsMax = _neibMaxNbNodes; + } + _d_nnz = (nbVoisinsMax+1)*_nVar; + _o_nnz = nbVoisinsMax *_nVar; _meshSet=true; } @@ -742,36 +766,36 @@ bool StationaryDiffusionEquation::solveStationaryProblem() bool stop=false; // Does the Problem want to stop (error) ? bool converged=false; // has the newton scheme converged (end) ? - cout<< "!!! Running test case "<< _fileName << " using "; + PetscPrintf(PETSC_COMM_WORLD,"!!! Running test case %s using ",_fileName); *_runLogFile<< "!!! Running test case "<< _fileName<< " using "; if(!_FECalculation) { - cout<< "Finite volumes method"<close(); throw CdmathException("Failed computing diffusion matrix"); } computeRHS(stop);//For the moment the heat power does not depend on the unknown temperature (linear RHS) if (stop){ - cout << "Error : failed computing right hand side, stopping calculation"<< endl; + PetscPrintf(PETSC_COMM_WORLD,"Error : failed computing right hand side, stopping calculation\n"); *_runLogFile << "Error : failed computing right hand side, stopping calculation"<< endl; throw CdmathException("Failed computing right hand side"); } stop = iterateNewtonStep(converged); if (stop){ - cout << "Error : failed solving linear system, stopping calculation"<< endl; + PetscPrintf(PETSC_COMM_WORLD,"Error : failed solving linear system, stopping calculation\n"); *_runLogFile << "Error : failed linear system, stopping calculation"<< endl; _runLogFile->close(); throw CdmathException("Failed solving linear system"); @@ -813,7 +837,7 @@ bool StationaryDiffusionEquation::solveStationaryProblem() } void StationaryDiffusionEquation::save(){ - cout<< "Saving numerical results"<1){ + VecScatterBegin(_scat,_Tk,_Tk_seq,INSERT_VALUES,SCATTER_FORWARD); + VecScatterEnd( _scat,_Tk,_Tk_seq,INSERT_VALUES,SCATTER_FORWARD); + } - _VV.setInfoOnComponent(0,"Temperature_(K)"); - switch(_saveFormat) - { - case VTK : - _VV.writeVTK(resultFile); - break; - case MED : - _VV.writeMED(resultFile); - break; - case CSV : - _VV.writeCSV(resultFile); - break; - } -} -Field -StationaryDiffusionEquation::getOutputTemperatureField() -{ - if(!_computationCompletedSuccessfully) - throw("Computation not performed yet or failed. No temperature field available"); - else - return _VV; + if(_verbose or _system) + VecView(_Tk,PETSC_VIEWER_STDOUT_WORLD); + + if(_mpi_rank==0){ + double Ti; + if(!_FECalculation) + for(int i=0; i<_Nmailles; i++) + { + if(_mpi_size>1) + VecGetValues(_Tk_seq, 1, &i, &Ti); + else + VecGetValues(_Tk , 1, &i, &Ti); + _VV(i)=Ti; + } + else + { + int globalIndex; + for(int i=0; i<_NunknownNodes; i++) + { + if(_mpi_size>1) + VecGetValues(_Tk_seq, 1, &i, &Ti); + else + VecGetValues(_Tk , 1, &i, &Ti); + globalIndex = DiffusionEquation::globalNodeIndex(i, _dirichletNodeIds); + _VV(globalIndex)=Ti; + } + + Node Ni; + string nameOfGroup; + for(int i=0; i<_NdirichletNodes; i++) + { + Ni=_mesh.getNode(_dirichletNodeIds[i]); + nameOfGroup = Ni.getGroupName(); + _VV(_dirichletNodeIds[i])=_limitField[nameOfGroup].T; + } + } + + _VV.setInfoOnComponent(0,"Temperature_(K)"); + switch(_saveFormat) + { + case VTK : + _VV.writeVTK(resultFile); + break; + case MED : + _VV.writeMED(resultFile); + break; + case CSV : + _VV.writeCSV(resultFile); + break; + } + } } void StationaryDiffusionEquation::terminate() @@ -885,6 +914,15 @@ void StationaryDiffusionEquation::terminate() VecDestroy(&_deltaT); VecDestroy(&_b); MatDestroy(&_A); + if(_mpi_size>1 && _mpi_rank == 0) + VecDestroy(&_Tk_seq); + + PetscBool petscInitialized; + PetscInitialized(&petscInitialized); + if(petscInitialized) + PetscFinalize(); + + delete _runLogFile; } void StationaryDiffusionEquation::setDirichletValues(map< int, double> dirichletBoundaryValues) @@ -953,3 +991,64 @@ StationaryDiffusionEquation::getEigenvectorsField(int nev, EPSWhich which, doubl return my_eigenfield; } + +Field& +StationaryDiffusionEquation::getOutputTemperatureField() +{ + if(!_computationCompletedSuccessfully) + throw("Computation not performed yet or failed. No temperature field available"); + else + return _VV; +} + +Field& +StationaryDiffusionEquation::getRodTemperatureField() +{ + return getOutputTemperatureField(); +} + +vector +StationaryDiffusionEquation::getInputFieldsNames() +{ + vector result(2); + + result[0]="FluidTemperature"; + result[1]="HeatPower"; + + return result; +} +vector +StationaryDiffusionEquation::getOutputFieldsNames() +{ + vector result(1); + + result[0]="RodTemperature"; + + return result; +} + +Field& +StationaryDiffusionEquation::getOutputField(const string& nameField ) +{ + if(nameField=="RodTemperature" || nameField=="RODTEMPERATURE" || nameField=="TEMPERATURECOMBUSTIBLE" || nameField=="TemperatureCombustible" ) + return getRodTemperatureField(); + else + { + cout<<"Error : Field name "<< nameField << " does not exist, call getOutputFieldsNames first" << endl; + throw CdmathException("DiffusionEquation::getOutputField error : Unknown Field name"); + } +} + +void +StationaryDiffusionEquation::setInputField(const string& nameField, Field& inputField ) +{ + if(nameField=="FluidTemperature" || nameField=="FLUIDTEMPERATURE" || nameField=="TemperatureFluide" || nameField=="TEMPERATUREFLUIDE") + return setFluidTemperatureField( inputField) ; + else if(nameField=="HeatPower" || nameField=="HEATPOWER" || nameField=="PuissanceThermique" || nameField=="PUISSANCETHERMIQUE" ) + return setHeatPowerField( inputField ); + else + { + cout<<"Error : Field name "<< nameField << " is not an input field name, call getInputFieldsNames first" << endl; + throw CdmathException("StationaryDiffusionEquation::setInputField error : Unknown Field name"); + } +} diff --git a/CoreFlows/Models/src/TransportEquation.cxx b/CoreFlows/Models/src/TransportEquation.cxx index de715fd..615398a 100755 --- a/CoreFlows/Models/src/TransportEquation.cxx +++ b/CoreFlows/Models/src/TransportEquation.cxx @@ -53,51 +53,74 @@ TransportEquation::TransportEquation(phase fluid, pressureMagnitude pEstimate,ve _dt_transport=0; _dt_src=0; _transportMatrixSet=false; + _FECalculation=false;//Only finite volumes available + _rodTemperatureFieldSet=false; + _rodTemperature=0; } void TransportEquation::initialize() { - if(!_initialDataSet) - throw CdmathException("TransportEquation::initialize() set initial data first"); - else if (_VV.getTypeOfField() != CELLS) - throw CdmathException("TransportEquation::initialize() Initial data should be a field on CELLS, not NODES, neither FACES"); - else - cout<<"Initialising the transport of a fluid enthalpy"<1 && _mpi_rank == 0) + VecCreateSeq(PETSC_COMM_SELF,_globalNbUnknowns,&_Hn_seq);//For saving results on proc 0 + VecScatterCreateToZero(_Hn,&_scat,&_Hn_seq); + //Linear solver KSPCreate(PETSC_COMM_SELF, &_ksp); KSPSetType(_ksp, _ksptype); @@ -113,127 +136,144 @@ void TransportEquation::initialize() double TransportEquation::computeTimeStep(bool & stop){ if(!_transportMatrixSet) _dt_transport=computeTransportMatrix(); + _dt_src=computeRHS(); + if(_verbose or _system) + { + PetscPrintf(PETSC_COMM_WORLD,"Right hand side of the linear system\n"); + VecView(_b,PETSC_VIEWER_STDOUT_WORLD); + } + stop=false; return min(_dt_transport,_dt_src); } double TransportEquation::computeTransportMatrix(){ - long nbFaces = _mesh.getNumberOfFaces(); - Face Fj; - Cell Cell1,Cell2; - string nameOfGroup; - double inv_dxi, inv_dxj; - Vector normale(_Ndim); - double un, hk; - PetscInt idm, idn; - std::vector< int > idCells; MatZeroEntries(_A); VecZeroEntries(_b0); - for (int j=0; j1){ - for(int l=0; l idCells; + for (int j=0; j1){ + for(int l=0; l_maxvp) + _maxvp=abs(un); + + // compute 1/dxi = volume of Ci/area of Fj + if (_Ndim > 1) + inv_dxi = Fj.getMeasure()/Cell1.getMeasure(); + else + inv_dxi = 1/Cell1.getMeasure(); + + // If Fj is on the boundary + if (Fj.getNumberOfCells()==1) { + if(_verbose && (_nbTimeStep-1)%_freqSave ==0) + { + cout << "face numero " << j << " cellule frontiere " << idCells[0] << " ; vecteur normal=("; + for(int p=0; p<_Ndim; p++) + cout << normale[p] << ","; + cout << ") "<0){ + MatSetValue(_A,idm,idm,inv_dxi*un, ADD_VALUES); + } + else{ + hk=_limitField[nameOfGroup].h; + VecSetValue(_b0,idm,-inv_dxi*un*hk, ADD_VALUES); + } + } + else {//bcType=NoneBCTransport + cout<<"!!!!!!!!!!!!!!! Error TransportEquation::computeTransportMatrix() !!!!!!!!!!"< 1) + inv_dxj = Fj.getMeasure()/Cell2.getMeasure(); else - normale[0] = -1; - } - } - //Compute velocity at the face Fj - un=normale*_vitesseTransport; - if(abs(un)>_maxvp) - _maxvp=abs(un); - - // compute 1/dxi = volume of Ci/area of Fj - if (_Ndim > 1) - inv_dxi = Fj.getMeasure()/Cell1.getMeasure(); - else - inv_dxi = 1/Cell1.getMeasure(); - - // If Fj is on the boundary - if (Fj.getNumberOfCells()==1) { - if(_verbose && (_nbTimeStep-1)%_freqSave ==0) - { - cout << "face numero " << j << " cellule frontiere " << idCells[0] << " ; vecteur normal=("; - for(int p=0; p<_Ndim; p++) - cout << normale[p] << ","; - cout << ") "<0){ MatSetValue(_A,idm,idm,inv_dxi*un, ADD_VALUES); + MatSetValue(_A,idn,idm,-inv_dxj*un, ADD_VALUES); } else{ - hk=_limitField[nameOfGroup].h; - VecSetValue(_b0,idm,-inv_dxi*un*hk, ADD_VALUES); + MatSetValue(_A,idm,idn,inv_dxi*un, ADD_VALUES); + MatSetValue(_A,idn,idn,-inv_dxj*un, ADD_VALUES); } } - else {//bcType=NoneBCTransport - cout<<"!!!!!!!!!!!!!!! Error TransportEquation::computeTransportMatrix() !!!!!!!!!!"< 1) - inv_dxj = Fj.getMeasure()/Cell2.getMeasure(); else - inv_dxj = 1/Cell2.getMeasure(); - - if(un>0){ - MatSetValue(_A,idm,idm,inv_dxi*un, ADD_VALUES); - MatSetValue(_A,idn,idm,-inv_dxj*un, ADD_VALUES); - } - else{ - MatSetValue(_A,idm,idn,inv_dxi*un, ADD_VALUES); - MatSetValue(_A,idn,idn,-inv_dxj*un, ADD_VALUES); - } + throw CdmathException("TransportEquation::ComputeTimeStep(): incompatible number of cells around a face"); } - else - throw CdmathException("TransportEquation::ComputeTimeStep(): incompatible number of cells around a face"); } - MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY); - MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY); - VecAssemblyBegin(_b0); - VecAssemblyEnd(_b0); + _transportMatrixSet=true; + + MPI_Bcast(&_maxvp, 1, MPI_DOUBLE, 0, PETSC_COMM_WORLD); + PetscPrintf(PETSC_COMM_WORLD, "Maximum speed is %.2f, CFL = %.2f, Delta x = %.2f\n",_maxvp,_cfl,_minl); + + MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY); + MatAssemblyEnd( _A, MAT_FINAL_ASSEMBLY); + VecAssemblyBegin(_b0); + VecAssemblyEnd( _b0); + if(abs(_maxvp)<_precision) throw CdmathException("TransportEquation::computeTransportMatrix(): maximum eigenvalue for time step is zero"); else @@ -241,57 +281,103 @@ double TransportEquation::computeTransportMatrix(){ } double TransportEquation::computeRHS(){ double rhomin=INFINITY; + VecCopy(_b0,_b); - VecAssemblyBegin(_b); - if(_system) - cout<<"second membre of transport problem"<1){ + VecScatterBegin(_scat,_Hn,_Hn_seq,INSERT_VALUES,SCATTER_FORWARD); + VecScatterEnd( _scat,_Hn,_Hn_seq,INSERT_VALUES,SCATTER_FORWARD); + } + + if(_verbose or _system) + { + PetscPrintf(PETSC_COMM_WORLD,"Unknown of the linear system :\n"); + VecView(_Hn,PETSC_VIEWER_STDOUT_WORLD); + } + + if(_mpi_rank == 0) { - VecGetValues(_Hk, 1, &i, &hi); - _VV(i)=hi; - _TT(i)=temperature(hi); - _Alpha(i)=voidFraction(hi); - _Rho(i)=density(_Alpha(i)); + double hi; + for(int i=0; i<_Nmailles; i++) + { + if(_mpi_size>1) + VecGetValues(_Hn_seq, 1, &i, &hi); + else + VecGetValues(_Hn , 1, &i, &hi); + _VV(i)=hi; + _TT(i)=temperature(hi); + _Alpha(i)=voidFraction(hi); + _Rho(i)=density(_Alpha(i)); + } } } bool TransportEquation::initTimeStep(double dt){ - _dt = dt; if(_verbose && (_nbTimeStep-1)%_freqSave ==0) - MatView(_A,PETSC_VIEWER_STDOUT_SELF); - MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY); - MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY); - - if(_timeScheme == Implicit) - MatShift(_A,1/_dt); + { + PetscPrintf(PETSC_COMM_WORLD,"Matrix of the linear system\n"); + MatView(_A,PETSC_VIEWER_STDOUT_WORLD); + } + if(_dt>0 and dt>0) + { + //Remove the contribution from dt to prepare for new time step. The diffusion matrix is not recomputed + if(_timeScheme == Implicit) + MatShift(_A,-1/_dt+1/dt); + //No need to remove the contribution to the right hand side since it is recomputed from scratch at each time step + } + else if(dt>0)//_dt==0, first time step + { + if(_timeScheme == Implicit) + MatShift(_A,1/dt); + } + else//dt<=0 + { + PetscPrintf(PETSC_COMM_WORLD,"TransportEquation::initTimeStep %.2f = \n",dt); + throw CdmathException("Error TransportEquation::initTimeStep : cannot set time step to zero"); + } + //At this stage _b contains _b0 + power + heat exchange + VecAXPY(_b, 1/dt, _Hn); + + _dt=dt; + return _dt>0; } void TransportEquation::abortTimeStep(){ + //Remove contribution of dt to the RHS VecAXPY(_b, -1/_dt, _Hn); - MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY); - MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY); + //Remove contribution of dt to the matrix if(_timeScheme == Implicit) MatShift(_A,-1/_dt); + _dt = 0; } @@ -310,13 +396,13 @@ bool TransportEquation::iterateTimeStep(bool &converged) VecAXPY(_b, 1/_dt, _Hn); if(_system) { - cout << "Vecteur Hn: " << endl; + PetscPrintf(PETSC_COMM_WORLD,"Vecteur Hn : \n"); VecView(_Hn, PETSC_VIEWER_STDOUT_WORLD); cout << endl; - cout<<"Vecteur _b "<=_maxPetscIts) { - cout<<"Systeme lineaire : pas de convergence de Petsc. Itérations maximales "<<_maxPetscIts<<" atteintes"<=_timeMax || _nbTimeStep>=_maxNbOfTimeStep) { - if(_Alpha(i)>alphamax) - alphamax=_Alpha(i); - if(_Alpha(i)alphamax) + alphamax=_Alpha(i); + if(_Alpha(i) TransportEquation::getOutputFieldsNames() +vector TransportEquation::getInputFieldsNames() { vector result(2); + result[0]="HeatPower"; + result[1]="RodTemperature"; + + return result; +} + +vector TransportEquation::getOutputFieldsNames() +{ + vector result(4); + result[0]="Enthalpy"; result[1]="FluidTemperature"; + result[2]="VoidFraction"; + result[3]="Density"; return result; } Field& TransportEquation::getOutputField(const string& nameField ) { - if(nameField=="FluidTemperature" || nameField=="FLUIDTEMPERATURE" ) + if(nameField=="FluidTemperature" || nameField=="FLUIDTEMPERATURE" || nameField=="TemperatureFluide" || nameField=="TEMPERATUREFLUIDE") return getFluidTemperatureField(); - else if(nameField=="Enthalpy" || nameField=="ENTHALPY" || nameField=="Enthalpie" || nameField=="ENTHALPY" ) + else if(nameField=="Enthalpy" || nameField=="ENTHALPY" || nameField=="Enthalpie" || nameField=="ENTHALPIE" ) return getEnthalpyField(); - else + else if(nameField=="VoidFraction" || nameField=="VOIDFRACTION" || nameField=="TauxDeVide" || nameField=="TAUXDEVIDE") + return getVoidFractionField(); + else if(nameField=="Density" || nameField=="DENSITY" || nameField=="Densité" || nameField=="DENSITE" ) + return getDensityField(); + else { cout<<"Error : Field name "<< nameField << " does not exist, call getOutputFieldsNames first" << endl; throw CdmathException("TransportEquation::getOutputField error : Unknown Field name"); } } +void +TransportEquation::setInputField(const string& nameField, Field& inputField ) +{ + if(nameField=="RodTemperature" || nameField=="RODTEMPERATURE" || nameField=="TemperatureCombustible" || nameField=="TEMPERATURECOMBUSTIBLE") + return setRodTemperatureField( inputField) ; + else if(nameField=="HeatPower" || nameField=="HEATPOWER" || nameField=="PuissanceThermique" || nameField=="PUISSANCETHERMIQUE" ) + return setHeatPowerField( inputField ); + else + { + cout<<"Error : Field name "<< nameField << " is not an input field name, call getInputFieldsNames first" << endl; + throw CdmathException("TransportEquation::setInputField error : Unknown Field name"); + } +} -- 2.39.2