diff --git a/include/dft.h b/include/dft.h index 98b9d1903..0b5d9dd1f 100644 --- a/include/dft.h +++ b/include/dft.h @@ -215,6 +215,25 @@ namespace dftfe distributedCPUVec & fvSpin0, distributedCPUVec & fvSpin1); + void copyDensityToVector( std::shared_ptr>> &rhoValues, + std::vector &rhoValuesVector); + + void copyDensityFromVector( std::shared_ptr>> &rhoValues, + std::vector &rhoValuesVector); + + void copyGradDensityToVector( std::shared_ptr>> &gradRhoValues, + std::vector &gradRhoValuesVector); + + void copyGradDensityFromVector( std::shared_ptr>> &gradRhoValues, + std::vector &gradRhoValuesVector); + + void computeTotalDensityFromSpinPolarised(std::shared_ptr>> &rhoValues, + std::shared_ptr>> &rhoSpinValues); + + void computeTotalGradDensityFromSpinPolarised(std::shared_ptr>> &gradRhoValues, + std::shared_ptr>> &gradRhoSpinValues); + + void computeJxWForRho(std::vector &vecJxW); void initializeKohnShamDFTOperator(const bool initializeCublas = true); diff --git a/include/mixingClass.h b/include/mixingClass.h index f31075a2c..48a8ee540 100644 --- a/include/mixingClass.h +++ b/include/mixingClass.h @@ -24,102 +24,77 @@ namespace dftfe { + + enum mixingVariable + { + rho, + gradRho + }; class MixingScheme { public: - MixingScheme(const dftParameters &dftParam, - const MPI_Comm &mpi_comm_domain); - - void init(const dealii::MatrixFree<3, double> & matrixFreeData, - const unsigned int matrixFreeVectorComponent, - const unsigned int matrixFreeQuadratureComponent, - excManager * excManagerPtr); - - void dotProduct(const std::deque> &inHist, - const std::deque> &outHist, - std::vector &A, - std::vector &c, - std::string fieldName); - - void copyDensityToInHist(std::shared_ptr>> &rhoInValues); - - void copyDensityToOutHist(std::shared_ptr>> &rhoOutValues); + MixingScheme(const MPI_Comm &mpi_comm_domain); unsigned int lengthOfHistory(); void computeAndersonMixingCoeff(); + void popOldHistory(); - void popRhoInHist(); - - void popRhoOutHist(); - - void copyDensityFromInHist(std::shared_ptr>> &rhoInValues); - void copyDensityFromOutHist(std::shared_ptr>> &rhoOutValues); - - void copyGradDensityFromInHist(std::shared_ptr>> &gradInput); - void copyGradDensityFromOutHist(std::shared_ptr>> &gradOutput); - - void copySpinGradDensityFromInHist(std::shared_ptr>> &gradInputSpin); - void copySpinGradDensityFromOutHist(std::shared_ptr>> &gradOutputSpin); - - void copyGradDensityToInHist(std::shared_ptr>> &gradInput); - void copyGradDensityToOutHist(std::shared_ptr>> &gradOutput); - - void copySpinGradDensityToInHist(std::shared_ptr>> &gradInputSpin); - void copySpinGradDensityToOutHist(std::shared_ptr>> &gradOutputSpin); + void clearHistory(); - void computeMixingMatricesDensity(const std::deque> &inHist, - const std::deque> &outHist, - std::vector &A, - std::vector &c); + void addMixingVariable(mixingVariable &mixingVariableList, + std::vector &weightDotProducts, + bool performMPIReduce, + double mixingValue); - double mixDensity(std::shared_ptr>> &rhoInValues, - std::shared_ptr>> &rhoOutValues, - std::shared_ptr>> &rhoInValuesSpinPolarized, - std::shared_ptr>> &rhoOutValuesSpinPolarized, - std::shared_ptr>> &gradRhoInValues, - std::shared_ptr>> &gradRhoOutValues, - std::shared_ptr>> &gradRhoInValuesSpinPolarized, - std::shared_ptr>> &gradRhoOutValuesSpinPolarized); + void addVariableToInHist(mixingVariable &mixingVariableName, + std::vector &inputVariableToInHist); - void popOldHistory(); + void addVariableToOutHist(mixingVariable &mixingVariableName, + std::vector &inputVariableToOutHist); - void clearHistory(); + double mixVariable(mixingVariable &mixingVariableName, + std::vector &outputVariable); - void popGradRhoInHist(); - void popGradRhoOutHist(); - void popGradRhoSpinInHist(); - void popGradRhoSpinOutHist(); private: + void computeMixingMatrices(const std::deque> &inHist, + const std::deque> &outHist, + const std::vector &weightDotProducts, + const bool isMPIAllReduce, + std::vector &A, + std::vector &c); + std::vector d_A, d_c; double d_cFinal; - std::deque> d_rhoInVals, - d_rhoOutVals, d_rhoInValsSpinPolarized, d_rhoOutValsSpinPolarized; - - std::deque>>> d_gradRhoInVals, - d_gradRhoInValsSpinPolarized, d_gradRhoOutVals, d_gradRhoOutValsSpinPolarized; - const dealii::MatrixFree<3,double> *d_matrixFreeDataPtr; - unsigned int d_matrixFreeVectorComponent,d_matrixFreeQuadratureComponent; - - unsigned int d_numberQuadraturePointsPerCell; - unsigned int d_numCells; - const dealii::DoFHandler<3> * d_dofHandler; - - std::vector d_vecJxW; - - const excManager * d_excManagerPtr; + std::map>> d_variableHistoryIn, d_variableHistoryOut ; + std:map> d_vectorDotProductWeights; + std:map d_performMPIReduce; const MPI_Comm d_mpi_comm_domain; - const dftParameters *d_dftParamsPtr; + std:map d_mixingParameter; }; } // end of namespace dftfe #endif // DFTFE_MIXINGCLASS_H + +// +//double mixDensity(std::shared_ptr>> &rhoInValues, +// std::shared_ptr>> &rhoOutValues, +// std::shared_ptr>> &rhoInValuesSpinPolarized, +// std::shared_ptr>> &rhoOutValuesSpinPolarized, +// std::shared_ptr>> &gradRhoInValues, +// std::shared_ptr>> &gradRhoOutValues, +// std::shared_ptr>> &gradRhoInValuesSpinPolarized, +// std::shared_ptr>> &gradRhoOutValuesSpinPolarized); + +// void copyDensityToInHist(std::shared_ptr>> &rhoInValues); +// +// void copyDensityToOutHist(std::shared_ptr>> &rhoOutValues); \ No newline at end of file diff --git a/src/dft/dft.cc b/src/dft/dft.cc index 0c11efdad..59fa9d8a6 100644 --- a/src/dft/dft.cc +++ b/src/dft/dft.cc @@ -146,8 +146,7 @@ namespace dftfe 0.0, 0.0, dftParams) - , d_mixingScheme(dftParams, - mpi_comm_domain) + , d_mixingScheme(mpi_comm_domain) #ifdef DFTFE_WITH_DEVICE , d_subspaceIterationSolverDevice(mpi_comm_parent, mpi_comm_domain, @@ -2165,11 +2164,24 @@ namespace dftfe 1e-3 : d_dftParamsPtr->chebyshevTolerance; - // initialise the mixing scheme - d_mixingScheme.init(matrix_free_data, - d_densityDofHandlerIndex, - d_densityQuadratureId, - d_excManagerPtr); + // initialise the variables in the mixing scheme + std::vector rhoJxW; + computeJxWForRho(rhoJxW); + d_mixingScheme.addMixingVariable(mixingVariable::rho, + rhoJxW, + true, // call MPI REDUCE while computing dot products + d_dftParamsPtr->mixingParameter); + + if(d_excManagerPtr->getDensityBasedFamilyType() == + densityFamilyType::GGA) + { + std::vector gradRhoJxW; + gradRhoJxW.resize(rhoJxW.size()*3); + d_mixingScheme.addMixingVariable(mixingVariable::gradRho, + gradRhoJxW, // this is just a dummy variable to amke it compatible with rho + false, // call MPI REDUCE while computing dot products + d_dftParamsPtr->mixingParameter); + } // // Begin SCF iteration @@ -2212,13 +2224,21 @@ namespace dftfe { if (d_dftParamsPtr->mixingMethod == "ANDERSON") { - d_mixingScheme.copyDensityToInHist(rhoInValuesSpinPolarized); - d_mixingScheme.copyDensityToOutHist(rhoOutValuesSpinPolarized); + std::vector rhoInOld, rhoOutOld; + + copyDensityToVector(rhoInValuesSpinPolarized,rhoInOld); + copyDensityToVector(rhoOutValuesSpinPolarized,rhoOutOld); + d_mixingScheme.addVariableToInHist(mixingVariable::rho,rhoInValues); + d_mixingScheme.addVariableToOutHist(mixingVariable::rho,rhoOutValues); if(d_excManagerPtr->getDensityBasedFamilyType() == densityFamilyType::GGA) { - d_mixingScheme.copySpinGradDensityToInHist(gradRhoInValuesSpinPolarized); - d_mixingScheme.copySpinGradDensityToOutHist(gradRhoOutValuesSpinPolarized); + std::vector gradRhoInOld, gradRhoOutOld; + + copyGradDensityToVector(gradRhoInValuesSpinPolarized,gradRhoInOld); + copyGradDensityToVector(gradRhoOutValuesSpinPolarized,gradRhoOutOld); + d_mixingScheme.addVariableToInHist(mixingVariable::gradRho,gradRhoInOld); + d_mixingScheme.addVariableToOutHist(mixingVariable::gradRho,gradRhoOutOld); } } norm = mixing_simple_spinPolarized(); @@ -2244,13 +2264,21 @@ namespace dftfe { if (d_dftParamsPtr->mixingMethod == "ANDERSON") { - d_mixingScheme.copyDensityToInHist(rhoInValues); - d_mixingScheme.copyDensityToOutHist(rhoOutValues); + std::vector rhoInOld, rhoOutOld; + + copyDensityToVector(rhoInValues,rhoInOld); + copyDensityToVector(rhoOutValues,rhoOutOld); + d_mixingScheme.addVariableToInHist(mixingVariable::rho,rhoInValues); + d_mixingScheme.addVariableToOutHist(mixingVariable::rho,rhoOutValues); if(d_excManagerPtr->getDensityBasedFamilyType() == densityFamilyType::GGA) { - d_mixingScheme.copyGradDensityToInHist(gradRhoInValues); - d_mixingScheme.copyGradDensityToOutHist(gradRhoOutValues); + std::vector gradRhoInOld, gradRhoOutOld; + + copyGradDensityToVector(gradRhoInValues,gradRhoInOld); + copyGradDensityToVector(gradRhoOutValues,gradRhoOutOld); + d_mixingScheme.addVariableToInHist(mixingVariable::gradRho,gradRhoInOld); + d_mixingScheme.addVariableToOutHist(mixingVariable::gradRho,gradRhoOutOld); } } norm = mixing_simple(); @@ -2270,26 +2298,45 @@ namespace dftfe { if (d_dftParamsPtr->mixingMethod == "ANDERSON") { - d_mixingScheme.copyDensityToInHist(rhoInValuesSpinPolarized); - d_mixingScheme.copyDensityToOutHist(rhoOutValuesSpinPolarized); - if(d_excManagerPtr->getDensityBasedFamilyType() == - densityFamilyType::GGA) - { - d_mixingScheme.copySpinGradDensityToInHist(gradRhoInValuesSpinPolarized); - d_mixingScheme.copySpinGradDensityToOutHist(gradRhoOutValuesSpinPolarized); - } + + std::vector rhoInOld, rhoOutOld; + + copyDensityToVector(rhoInValuesSpinPolarized,rhoInOld); + copyDensityToVector(rhoOutValuesSpinPolarized,rhoOutOld); + d_mixingScheme.addVariableToInHist(mixingVariable::rho,rhoInValues); + d_mixingScheme.addVariableToOutHist(mixingVariable::rho,rhoOutValues); + if(d_excManagerPtr->getDensityBasedFamilyType() == + densityFamilyType::GGA) + { + std::vector gradRhoInOld, gradRhoOutOld; + + copyGradDensityToVector(gradRhoInValuesSpinPolarized,gradRhoInOld); + copyGradDensityToVector(gradRhoOutValuesSpinPolarized,gradRhoOutOld); + d_mixingScheme.addVariableToInHist(mixingVariable::gradRho,gradRhoInOld); + d_mixingScheme.addVariableToOutHist(mixingVariable::gradRho,gradRhoOutOld); + } d_mixingScheme.popOldHistory(); d_mixingScheme.computeAndersonMixingCoeff(); - norm = d_mixingScheme.mixDensity(rhoInValues, - rhoOutValues, - rhoInValuesSpinPolarized, - rhoOutValuesSpinPolarized, - gradRhoInValues, - gradRhoOutValues, - gradRhoInValuesSpinPolarized, - gradRhoOutValuesSpinPolarized); + std::vector rhoInNew; + norm = d_mixingScheme.mixVariable(mixingVariable::rho, + rhoInNew); + + copyDensityFromVector(rhoInValuesSpinPolarized,rhoInNew); + computeTotalDensityFromSpinPolarised(rhoInValues,rhoInValuesSpinPolarized); + + if(d_excManagerPtr->getDensityBasedFamilyType() == + densityFamilyType::GGA) + { + std::vector gradRhoInNew; + d_mixingScheme.mixVariable(mixingVariable::gradRho, + gradRhoInNew); + copyGradDensityFromVector(gradRhoInValuesSpinPolarized,gradRhoInNew); + computeTotalGradDensityFromSpinPolarised(gradRhoInValues,gradRhoInValuesSpinPolarized); + } + + // double normRhoInAfterMixing = 0.0; // double normGradRhoInAfterMixing = 0.0; // double normSpinRhoInAfterMixing = 0.0; @@ -2376,26 +2423,39 @@ namespace dftfe { if (d_dftParamsPtr->mixingMethod == "ANDERSON") { - d_mixingScheme.copyDensityToInHist(rhoInValues); - d_mixingScheme.copyDensityToOutHist(rhoOutValues); - if(d_excManagerPtr->getDensityBasedFamilyType() == - densityFamilyType::GGA) - { - d_mixingScheme.copyGradDensityToInHist(gradRhoInValues); - d_mixingScheme.copyGradDensityToOutHist(gradRhoOutValues); - } + std::vector rhoInOld, rhoOutOld; + + copyDensityToVector(rhoInValues,rhoInOld); + copyDensityToVector(rhoOutValues,rhoOutOld); + d_mixingScheme.addVariableToInHist(mixingVariable::rho,rhoInValues); + d_mixingScheme.addVariableToOutHist(mixingVariable::rho,rhoOutValues); + if(d_excManagerPtr->getDensityBasedFamilyType() == + densityFamilyType::GGA) + { + std::vector gradRhoInOld, gradRhoOutOld; + + copyGradDensityToVector(gradRhoInValues,gradRhoInOld); + copyGradDensityToVector(gradRhoOutValues,gradRhoOutOld); + d_mixingScheme.addVariableToInHist(mixingVariable::gradRho,gradRhoInOld); + d_mixingScheme.addVariableToOutHist(mixingVariable::gradRho,gradRhoOutOld); + } d_mixingScheme.popOldHistory(); d_mixingScheme.computeAndersonMixingCoeff(); - norm = d_mixingScheme.mixDensity(rhoInValues, - rhoOutValues, - rhoInValuesSpinPolarized, - rhoOutValuesSpinPolarized, - gradRhoInValues, - gradRhoOutValues, - gradRhoInValuesSpinPolarized, - gradRhoOutValuesSpinPolarized); + std::vector rhoInNew; + norm = d_mixingScheme.mixVariable(mixingVariable::rho, + rhoInNew); + + copyDensityFromVector(rhoInValues,rhoInNew); + if(d_excManagerPtr->getDensityBasedFamilyType() == + densityFamilyType::GGA) + { + std::vector gradRhoInNew; + d_mixingScheme.mixVariable(mixingVariable::gradRho, + gradRhoInNew); + copyGradDensityFromVector(gradRhoInValues,gradRhoInNew); + } // double normRhoInAfterMixing = 0.0; // double normGradRhoInAfterMixing = 0.0; // const dealii::Quadrature<3> &quadrature = @@ -4768,6 +4828,317 @@ namespace dftfe << std::endl; } + template + void dftClass::copyDensityToVector( std::shared_ptr>> &rhoValues, + std::vector &rhoValuesVector) + { + unsigned int numCells = matrix_free_data.n_physical_cells(); + + const dealii::Quadrature<3> &quadratureDensity = + matrix_free_data.get_quadrature( + d_densityQuadratureId); + + unsigned int numQuadPoints = quadratureDensity.size(); + + if (d_dftParamsPtr->spinPolarized == 1) + { + numQuadPoints = 2*numQuadPoints; + } + + rhoValuesVector.resize(numCells*numQuadPoints); + std::fill(rhoValuesVector.begin(),rhoValuesVector.end(),0.0); + + unsigned int iElem = 0; + + const dealii::DoFHandler<3> * dofHandler = + &matrix_free_data.get_dof_handler(d_densityDofHandlerIndex); + + typename dealii::DoFHandler<3>::active_cell_iterator + cell = dofHandler->begin_active(), + endc = dofHandler->end(); + for (; cell != endc; ++cell) + { + if (cell->is_locally_owned()) + { + for(unsigned int iQuad = 0 ;iQuad < numQuadPoints; iQuad++) + { + rhoValuesVector[iElem*numQuadPoints + + iQuad] = + (*rhoValues)[cell->id()][iQuad]; + } + iElem++; + } + } + } + + template + void dftClass::copyDensityFromVector( std::shared_ptr>> &rhoValues, + std::vector &rhoValuesVector) + { + unsigned int numCells = matrix_free_data.n_physical_cells(); + + const dealii::Quadrature<3> &quadratureDensity = + matrix_free_data.get_quadrature( + d_densityQuadratureId); + + unsigned int numQuadPoints = quadratureDensity.size(); + + if (d_dftParamsPtr->spinPolarized == 1) + { + numQuadPoints = 2*numQuadPoints; + } + + unsigned int iElem = 0; + + const dealii::DoFHandler<3> * dofHandler = + &matrix_free_data.get_dof_handler(d_densityDofHandlerIndex); + + typename dealii::DoFHandler<3>::active_cell_iterator + cell = dofHandler->begin_active(), + endc = dofHandler->end(); + for (; cell != endc; ++cell) + { + if (cell->is_locally_owned()) + { + for(unsigned int iQuad = 0 ;iQuad < numQuadPoints; iQuad++) + { + (*rhoValues)[cell->id()][iQuad] = + rhoValuesVector[iElem*numQuadPoints + + iQuad]; + } + iElem++; + } + } + } + + template + void dftClass::copyGradDensityToVector( std::shared_ptr>> &gradRhoValues, + std::vector &gradRhoValuesVector) + { + unsigned int numCells = matrix_free_data.n_physical_cells(); + + const dealii::Quadrature<3> &quadratureDensity = + matrix_free_data.get_quadrature( + d_densityQuadratureId); + + unsigned int numQuadPoints = 3*quadratureDensity.size(); + + if (d_dftParamsPtr->spinPolarized == 1) + { + numQuadPoints = 2*numQuadPoints; + } + + gradRhoValuesVector.resize(numCells*numQuadPoints); + std::fill(gradRhoValuesVector.begin(),gradRhoValuesVector.end(),0.0); + + unsigned int iElem = 0; + + const dealii::DoFHandler<3> * dofHandler = + &matrix_free_data.get_dof_handler(d_densityDofHandlerIndex); + + typename dealii::DoFHandler<3>::active_cell_iterator + cell = dofHandler->begin_active(), + endc = dofHandler->end(); + for (; cell != endc; ++cell) + { + if (cell->is_locally_owned()) + { + for(unsigned int iQuad = 0 ;iQuad < numQuadPoints; iQuad++) + { + gradRhoValuesVector[iElem*numQuadPoints + + iQuad] = + (*gradRhoValues)[cell->id()][iQuad] ; + } + iElem++; + } + } + } + + template + void dftClass::copyGradDensityFromVector( std::shared_ptr>> &gradRhoValues, + std::vector &gradRhoValuesVector) + { + unsigned int numCells = matrix_free_data.n_physical_cells(); + + const dealii::Quadrature<3> &quadratureDensity = + matrix_free_data.get_quadrature( + d_densityQuadratureId); + + unsigned int numQuadPoints = 3*quadratureDensity.size(); + + if (d_dftParamsPtr->spinPolarized == 1) + { + numQuadPoints = 2*numQuadPoints; + } + + unsigned int iElem = 0; + + const dealii::DoFHandler<3> * dofHandler = + &matrix_free_data.get_dof_handler(d_densityDofHandlerIndex); + + typename dealii::DoFHandler<3>::active_cell_iterator + cell = dofHandler->begin_active(), + endc = dofHandler->end(); + for (; cell != endc; ++cell) + { + if (cell->is_locally_owned()) + { + for(unsigned int iQuad = 0 ;iQuad < numQuadPoints; iQuad++) + { + (*gradRhoValues)[cell->id()][iQuad] = + gradRhoValuesVector[iElem*numQuadPoints + + iQuad]; + } + iElem++; + } + } + } + + template + void dftClass::computeTotalDensityFromSpinPolarised(std::shared_ptr>> &rhoValues, + std::shared_ptr>> &rhoSpinValues) + { + unsigned int numCells = matrix_free_data.n_physical_cells(); + + const dealii::Quadrature<3> &quadratureDensity = + matrix_free_data.get_quadrature( + d_densityQuadratureId); + + unsigned int numQuadPoints = quadratureDensity.size(); + + unsigned int iElem = 0; + + const dealii::DoFHandler<3> * dofHandler = + &matrix_free_data.get_dof_handler(d_densityDofHandlerIndex); + + typename dealii::DoFHandler<3>::active_cell_iterator + cell = dofHandler->begin_active(), + endc = dofHandler->end(); + for (; cell != endc; ++cell) + { + if (cell->is_locally_owned()) + { + for(unsigned int iQuad = 0 ;iQuad < numQuadPoints; iQuad++) + { + (*rhoValues)[cell->id()][iQuad] = + (*rhoSpinValues)[cell->id()][2*iQuad+0] + + (*rhoSpinValues)[cell->id()][2*iQuad+0]; + } + iElem++; + } + } + } + + template + void dftClass::computeTotalGradDensityFromSpinPolarised(std::shared_ptr>> &gradRhoValues, + std::shared_ptr>> &gradRhoSpinValues) + { + unsigned int numCells = matrix_free_data.n_physical_cells(); + + const dealii::Quadrature<3> &quadratureDensity = + matrix_free_data.get_quadrature( + d_densityQuadratureId); + + unsigned int numQuadPoints = quadratureDensity.size(); + + unsigned int iElem = 0; + + const dealii::DoFHandler<3> * dofHandler = + &matrix_free_data.get_dof_handler(d_densityDofHandlerIndex); + + typename dealii::DoFHandler<3>::active_cell_iterator + cell = dofHandler->begin_active(), + endc = dofHandler->end(); + for (; cell != endc; ++cell) + { + if (cell->is_locally_owned()) + { + for(unsigned int iQuad = 0 ;iQuad < numQuadPoints; iQuad++) + { + ((*gradRhoValues)[cell->id()][3 * q_point + 0]) = + ((*gradRhoSpinValues)[cell->id()] + [6 * q_point + 0]) + + ((*gradRhoSpinValues)[cell->id()] + [6 * q_point + 3]); + ((*gradRhoValues)[cell->id()][3 * q_point + 1]) = + ((*gradRhoSpinValues)[cell->id()] + [6 * q_point + 1]) + + ((*gradRhoInValuesSpinPolarized)[cell->id()] + [6 * q_point + 4]); + ((*gradRhoValues)[cell->id()][3 * q_point + 2]) = + ((*gradRhoSpinValues)[cell->id()] + [6 * q_point + 2]) + + ((*gradRhoSpinValues)[cell->id()] + [6 * q_point + 5]); + } + iElem++; + } + } + } + + template + void dftClass::computeJxWForRho(std::vector &vecJxW) + { + unsigned int numCells = matrix_free_data.n_physical_cells(); + + const dealii::Quadrature<3> &quadratureDensity = + matrix_free_data.get_quadrature( + d_densityQuadratureId); + + unsigned int numQuadPoints = quadratureDensity.size(); + + + + const dealii::DoFHandler<3> * dofHandler = + &matrix_free_data.get_dof_handler(d_densityDofHandlerIndex); + + dealii::FEValues<3> fe_values(dofHandler->get_fe(), + quadratureDensity, + dealii::update_JxW_values); + + typename dealii::DoFHandler<3>::active_cell_iterator + cell = dofHandler->begin_active(), + endc = dofHandler->end(); + unsigned int iElem = 0; + if (d_dftParamsPtr->spinPolarized == 1) + { + vecJxW.reinit(numCells*numQuadPoints*2); + std::fill(vecJxW.begin(),vecJxW.end(),0.0); + for (; cell != endc; ++cell) + { + if (cell->is_locally_owned()) + { + fe_values.reinit(cell); + for(unsigned int iQuad = 0 ;iQuad < numQuadPoints; iQuad++) + { + vecJxW[iElem*numQuadPoints*2 + 2*iQuad + 0 ] = fe_values.JxW(iQuad); + vecJxW[iElem*numQuadPoints*2 + 2*iQuad + 1 ] = + vecJxW[iElem*numQuadPoints*2 + 2*iQuad + 0 ]; + + } + iElem++; + } + } + } + else + { + vecJxW.reinit(numCells*numQuadPoints); + std::fill(vecJxW.begin(),vecJxW.end(),0.0); + for (; cell != endc; ++cell) + { + if (cell->is_locally_owned()) + { + fe_values.reinit(cell); + for(unsigned int iQuad = 0 ;iQuad < numQuadPoints; iQuad++) + { + vecJxW[iElem*numQuadPoints + iQuad ] = fe_values.JxW(iQuad); + } + iElem++; + } + } + } + } + #include "dft.inst.cc" } // namespace dftfe diff --git a/src/dft/mixingClass.cc b/src/dft/mixingClass.cc index 7a03ad653..df42882a0 100644 --- a/src/dft/mixingClass.cc +++ b/src/dft/mixingClass.cc @@ -23,200 +23,32 @@ namespace dftfe { - MixingScheme::MixingScheme(const dftParameters &dftParam, - const MPI_Comm &mpi_comm_domain): - d_dftParamsPtr(&dftParam), + MixingScheme::MixingScheme(const MPI_Comm &mpi_comm_domain): d_mpi_comm_domain(mpi_comm_domain) { } - void MixingScheme::copyDensityFromInHist(std::shared_ptr>> &rhoInValues) + void MixingScheme::addMixingVariable(mixingVariable &mixingVariableList, + std::vector &weightDotProducts, + bool performMPIReduce, + double mixingValue) { - unsigned int numQuadPointsPerCell = d_numberQuadraturePointsPerCell; - if (d_dftParamsPtr->spinPolarized == 1) - { - numQuadPointsPerCell = 2*numQuadPointsPerCell; - } - - unsigned int hist =d_rhoInVals.size(); - - unsigned int iElem = 0; - - typename dealii::DoFHandler<3>::active_cell_iterator - cell = d_dofHandler->begin_active(), - endc = d_dofHandler->end(); - for (; cell != endc; ++cell) - { - if (cell->is_locally_owned()) - { - for(unsigned int iQuad = 0 ;iQuad < numQuadPointsPerCell; iQuad++) - { - (*rhoInValues)[cell->id()][iQuad] = - d_rhoInVals[hist-1][iElem*numQuadPointsPerCell + - iQuad]; - } - iElem++; - } - } - - } - void MixingScheme::copyDensityFromOutHist(std::shared_ptr>> &rhoOutValues) - { - unsigned int numQuadPointsPerCell = d_numberQuadraturePointsPerCell; - if (d_dftParamsPtr->spinPolarized == 1) - { - numQuadPointsPerCell = 2*numQuadPointsPerCell; - } - - unsigned int hist =d_rhoOutVals.size(); - - unsigned int iElem = 0; - - d_dofHandler = - &d_matrixFreeDataPtr->get_dof_handler(d_matrixFreeVectorComponent); + d_variableHistoryIn.insert({mixingVariableList, std::deque>()}); - typename dealii::DoFHandler<3>::active_cell_iterator - cell = d_dofHandler->begin_active(), - endc = d_dofHandler->end(); - for (; cell != endc; ++cell) - { - if (cell->is_locally_owned()) - { - for(unsigned int iQuad = 0 ;iQuad < numQuadPointsPerCell; iQuad++) - { - (*rhoOutValues)[cell->id()][iQuad] = - d_rhoOutVals[hist-1][iElem*numQuadPointsPerCell + - iQuad]; - } - iElem++; - } - } + d_variableHistoryOut.insert({mixingVariableList, std::deque>()}); - } + d_vectorDotProductWeights.insert({mixingVariableList,weightDotProducts}); - void MixingScheme::copyGradDensityFromInHist(std::shared_ptr>> &gradInput) - { - gradInput = d_gradRhoInVals.back(); - } + d_performMPIReduce.insert({mixingVariableList,performMPIReduce}); - void MixingScheme::copyGradDensityFromOutHist(std::shared_ptr>> &gradOutput) - { - gradOutput = d_gradRhoOutVals.back(); + d_mixingParameter.insert({mixingVariableList,mixingValue}); } - void MixingScheme::copySpinGradDensityFromInHist(std::shared_ptr>> &gradInputSpin) - { - gradInputSpin = d_gradRhoInValsSpinPolarized.back(); - } - - void MixingScheme::copySpinGradDensityFromOutHist(std::shared_ptr>> &gradOutputSpin) - { - gradOutputSpin = d_gradRhoOutValsSpinPolarized.back(); - } - - void MixingScheme::copyGradDensityToInHist(std::shared_ptr>> &gradInput) - { - d_gradRhoInVals.push_back(gradInput); - } - void MixingScheme::copyGradDensityToOutHist(std::shared_ptr>> &gradOutput) - { - d_gradRhoOutVals.push_back(gradOutput); - } - - void MixingScheme::copySpinGradDensityToInHist(std::shared_ptr>> &gradInputSpin) - { - d_gradRhoInValsSpinPolarized.push_back(gradInputSpin); - } - void MixingScheme::copySpinGradDensityToOutHist(std::shared_ptr>> &gradOutputSpin) - { - d_gradRhoOutValsSpinPolarized.push_back(gradOutputSpin); - } - - - void MixingScheme::init(const dealii::MatrixFree<3, double> & matrixFreeData, - const unsigned int matrixFreeVectorComponent, - const unsigned int matrixFreeQuadratureComponent, - excManager * excManagerPtr) - { - - d_matrixFreeDataPtr = &matrixFreeData ; - d_matrixFreeVectorComponent = matrixFreeVectorComponent; - d_matrixFreeQuadratureComponent = matrixFreeQuadratureComponent; - d_excManagerPtr = excManagerPtr; - const dealii::Quadrature<3> &quadratureRhs = - d_matrixFreeDataPtr->get_quadrature( - d_matrixFreeQuadratureComponent); - - d_numberQuadraturePointsPerCell = quadratureRhs.size(); - - d_dofHandler = - &d_matrixFreeDataPtr->get_dof_handler(d_matrixFreeVectorComponent); - - dealii::FEValues<3> fe_values(d_dofHandler->get_fe(), - quadratureRhs, - dealii::update_JxW_values); - - unsigned int iElem = 0; - - d_numCells = d_matrixFreeDataPtr->n_physical_cells(); - - if(d_dftParamsPtr->spinPolarized == 0) - { - d_vecJxW.resize(d_numberQuadraturePointsPerCell* - d_numCells); - typename dealii::DoFHandler<3>::active_cell_iterator - cell = d_dofHandler->begin_active(), - endc = d_dofHandler->end(); - for (; cell != endc; ++cell) - { - if (cell->is_locally_owned()) - { - fe_values.reinit(cell); - for(unsigned int iQuad = 0 ;iQuad < d_numberQuadraturePointsPerCell; iQuad++) - { - d_vecJxW[iElem*d_numberQuadraturePointsPerCell + - iQuad] = - fe_values.JxW(iQuad); - } - iElem++; - } - } - } - else - { - d_vecJxW.resize(2* - d_numberQuadraturePointsPerCell* - d_numCells,0.0); - typename dealii::DoFHandler<3>::active_cell_iterator - cell = d_dofHandler->begin_active(), - endc = d_dofHandler->end(); - for (; cell != endc; ++cell) - { - if (cell->is_locally_owned()) - { - fe_values.reinit(cell); - for(unsigned int iQuad = 0 ;iQuad < d_numberQuadraturePointsPerCell; iQuad++) - { - d_vecJxW[2*iElem*d_numberQuadraturePointsPerCell + - 2*iQuad - +0] = - fe_values.JxW(iQuad); - - d_vecJxW[2*iElem*d_numberQuadraturePointsPerCell + - 2*iQuad - +1] = - fe_values.JxW(iQuad); - } - iElem++; - } - } - } -// std::cout<<" size of JxW = "<> &inHist, + void MixingScheme::computeMixingMatrices(const std::deque> &inHist, const std::deque> &outHist, + const std::vector &weightDotProducts, + const bool isMPIAllReduce, std::vector &A, std::vector &c) { @@ -228,59 +60,58 @@ namespace dftfe cDensity.resize(c.size()); std::fill(cDensity.begin(),cDensity.end(),0.0); - int N = d_rhoOutVals.size() - 1; + int N = inHist.size() - 1; // std::cout<<" size of hist = "< 0) numQuadPoints = inHist[0].size(); - if( numQuadPoints != d_vecJxW.size()) + if( numQuadPoints != weightDotProducts.size()) { std::cout<<" ERROR in vec size in mixing anderson \n"; } - double vecJxWSum = 0.0; for( unsigned int iQuad = 0; iQuad < numQuadPoints; iQuad++) { - vecJxWSum += d_vecJxW[iQuad]; - double Fn = d_rhoOutVals[N][iQuad] - d_rhoInVals[N][iQuad]; + double Fn = outHist[N][iQuad] - inHist[N][iQuad]; for (int m = 0; m < N; m++) { - double Fnm = d_rhoOutVals[N - 1 - m][iQuad] - - d_rhoInVals[N - 1 - m][iQuad]; + double Fnm = outHist[N - 1 - m][iQuad] - + inHist[N - 1 - m][iQuad]; for (int k = 0; k < N; k++) { - double Fnk = d_rhoOutVals[N - 1 - k][iQuad] - - d_rhoInVals[N - 1 - k][iQuad]; + double Fnk = outHist[N - 1 - k][iQuad] - + inHist[N - 1 - k][iQuad]; Adensity[k * N + m] += (Fn - Fnm) * (Fn - Fnk) * - d_vecJxW[iQuad]; // (m,k)^th entry + weightDotProducts[iQuad]; // (m,k)^th entry } cDensity[m] += - (Fn - Fnm) * (Fn)*d_vecJxW[iQuad]; // (m)^th entry + (Fn - Fnm) * (Fn)*weightDotProducts[iQuad]; // (m)^th entry } } - MPI_Allreduce(MPI_IN_PLACE, - &vecJxWSum, - 1, - MPI_DOUBLE, - MPI_SUM, - d_mpi_comm_domain); - -// std::cout<<" Sum of Jxw = "< ATotal(aSize), cTotal(cSize); std::fill(ATotal.begin(),ATotal.end(),0.0); std::fill(cTotal.begin(),cTotal.end(),0.0); - MPI_Allreduce( - &Adensity[0], &ATotal[0], aSize, MPI_DOUBLE, MPI_SUM, d_mpi_comm_domain); - MPI_Allreduce( - &cDensity[0], &cTotal[0], cSize, MPI_DOUBLE, MPI_SUM, d_mpi_comm_domain); + if (isMPIAllReduce) + { + MPI_Allreduce( + &Adensity[0], &ATotal[0], aSize, MPI_DOUBLE, MPI_SUM, d_mpi_comm_domain); + MPI_Allreduce( + &cDensity[0], &cTotal[0], cSize, MPI_DOUBLE, MPI_SUM, d_mpi_comm_domain); + } + else + { + ATotal = Adensity; + cTotal = cDensity; + } + + for (unsigned int i = 0 ; i < aSize; i++) { @@ -295,86 +126,17 @@ namespace dftfe } } - void MixingScheme::copyDensityToInHist(std::shared_ptr>> &rhoInValues) - { - std::vector latestHistRhoIn; - unsigned int numQuadPointsPerCell = d_numberQuadraturePointsPerCell; - if (d_dftParamsPtr->spinPolarized == 1) - { - numQuadPointsPerCell = 2*numQuadPointsPerCell; - } - - latestHistRhoIn.resize(d_numCells*numQuadPointsPerCell); - std::fill(latestHistRhoIn.begin(),latestHistRhoIn.end(),0.0); - - unsigned int iElem = 0; - typename dealii::DoFHandler<3>::active_cell_iterator - cell = d_dofHandler->begin_active(), - endc = d_dofHandler->end(); - for (; cell != endc; ++cell) - { - if (cell->is_locally_owned()) - { - for(unsigned int iQuad = 0 ;iQuad < numQuadPointsPerCell; iQuad++) - { - latestHistRhoIn[iElem*numQuadPointsPerCell + - iQuad] = - (*rhoInValues)[cell->id()][iQuad]; - } - iElem++; - } - } - - d_rhoInVals.push_back(latestHistRhoIn); - } - - void MixingScheme::copyDensityToOutHist( std::shared_ptr>> &rhoOutValues) - { - std::vector latestHistRhoOut; - unsigned int numQuadPointsPerCell = d_numberQuadraturePointsPerCell; - if (d_dftParamsPtr->spinPolarized == 1) - { - numQuadPointsPerCell = 2*numQuadPointsPerCell; - } - - latestHistRhoOut.resize(d_numCells*numQuadPointsPerCell); - std::fill(latestHistRhoOut.begin(),latestHistRhoOut.end(),0.0); - - unsigned int iElem = 0; - - d_dofHandler = - &d_matrixFreeDataPtr->get_dof_handler(d_matrixFreeVectorComponent); - - typename dealii::DoFHandler<3>::active_cell_iterator - cell = d_dofHandler->begin_active(), - endc = d_dofHandler->end(); - for (; cell != endc; ++cell) - { - if (cell->is_locally_owned()) - { - for(unsigned int iQuad = 0 ;iQuad < numQuadPointsPerCell; iQuad++) - { - latestHistRhoOut[iElem*numQuadPointsPerCell + - iQuad] = - (*rhoOutValues)[cell->id()][iQuad]; - } - iElem++; - } - } - - d_rhoOutVals.push_back(latestHistRhoOut); - } unsigned int MixingScheme::lengthOfHistory() { - return d_rhoInVals.size(); + return variableHistoryIn[mixingVariable::rho].size(); } void MixingScheme::computeAndersonMixingCoeff() { // initialize data structures - int N = d_rhoOutVals.size() - 1; + int N = variableHistoryIn[mixingVariable::rho].size() - 1; // pcout << "\nN:" << N << "\n"; int NRHS = 1, lda = N, ldb = N, info; std::vector ipiv(N); @@ -385,10 +147,15 @@ namespace dftfe for (int i = 0; i < ldb * NRHS; i++) d_c[i] = 0.0; - computeMixingMatricesDensity(d_rhoInVals, - d_rhoOutVals, - d_A, - d_c); + for (const auto& [key, value] : variableHistoryIn) + { + computeMixingMatrices(variableHistoryIn[key], + variableHistoryOut[key], + vectorDotProductWeights[key], + performMPIReduce[key], + d_A, + d_c); + } dgesv_(&N, &NRHS, &d_A[0], &lda, &ipiv[0], &d_c[0], &ldb, &info); @@ -404,627 +171,872 @@ namespace dftfe // std::cout<<"cFinal = "<>> &rhoInValues, - std::shared_ptr>> &rhoOutValues, - std::shared_ptr>> &rhoInValuesSpinPolarized, - std::shared_ptr>> &rhoOutValuesSpinPolarized, - std::shared_ptr>> &gradRhoInValues, - std::shared_ptr>> &gradRhoOutValues, - std::shared_ptr>> &gradRhoInValuesSpinPolarized, - std::shared_ptr>> &gradRhoOutValuesSpinPolarized) + void MixingScheme::addVariableToInHist(mixingVariable &mixingVariableName, + std::vector &inputVariableToInHist) { - double normValue = 0.0; - -// std::cout<<" size of rhoInHist = "< rhoInHistNorm(N+1,0.0), rhoOutHistNorm(N+1,0.0), gradRhoInHistNorm(N+1,0.0), gradRhoOutHistNorm(N+1,0.0), gradSpinRhoInHistNorm(N+1,0.0), gradSpinRhoOutHistNorm(N+1,0.0); - double rhoOutputNorm = 0.0, gradRhoOutputNorm = 0.0; - - std::fill(rhoInHistNorm.begin(),rhoInHistNorm.end(),0.0); - std::fill(rhoOutHistNorm.begin(),rhoOutHistNorm.end(),0.0); - std::fill(gradRhoInHistNorm.begin(),gradRhoInHistNorm.end(),0.0); - std::fill(gradRhoOutHistNorm.begin(),gradRhoOutHistNorm.end(),0.0); - std::fill(gradSpinRhoInHistNorm.begin(),gradSpinRhoInHistNorm.end(),0.0); - std::fill(gradSpinRhoOutHistNorm.begin(),gradSpinRhoOutHistNorm.end(),0.0); - - std::shared_ptr>> rhoInputValues, rhoOutputValues; - // create new rhoValue tables - std::map> rhoInValuesOld = *rhoInValues; - + variableHistoryIn[mixingVariableName].push_back(inputVariableToInHist); + } - unsigned int numQuadPointsPerCell = d_numberQuadraturePointsPerCell; - if(d_dftParamsPtr->spinPolarized == 0) - { - rhoInValuesOld = *rhoInValues; - rhoInputValues = rhoInValues; - rhoOutputValues = rhoOutValues; - } - if (d_dftParamsPtr->spinPolarized == 1) - { - numQuadPointsPerCell = 2*numQuadPointsPerCell; + void MixingScheme::addVariableToOutHist(mixingVariable &mixingVariableName, + std::vector &inputVariableToOutHist) + { + variableHistoryOut[mixingVariableName].push_back(inputVariableToOutHist); + } - rhoInValuesOld = *rhoInValuesSpinPolarized; - rhoInputValues = rhoInValuesSpinPolarized; - rhoOutputValues = rhoOutValuesSpinPolarized; - } + double MixingScheme::mixVariable(mixingVariable &mixingVariableName, + std::vector &outputVariable) + { + double normValue = 0.0; + int N = variableHistoryIn[mixingVariableName].size() - 1; + unsigned int lenVar = variableHistoryIn[mixingVariableName][0].size(); + outputVariable.resize(lenVar); + std::fill(outputVariable.begin(),outputVariable.end(),0.0); - // implement anderson mixing - typename dealii::DoFHandler<3>::active_cell_iterator - cell = d_dofHandler->begin_active(), - endc = d_dofHandler->end(); - unsigned int iElem = 0; - for (; cell != endc; ++cell) + for( unsigned int iQuad = 0; iQuad < lenVar; iQuad++) { - if (cell->is_locally_owned()) - { - (*rhoInputValues)[cell->id()].resize(numQuadPointsPerCell); - - for (unsigned int q_point = 0; q_point < numQuadPointsPerCell; ++q_point) - { - - for(unsigned int iHist = 0; iHist < N+1; iHist++) - { - rhoInHistNorm[iHist] += std::abs(d_rhoInVals[iHist][iElem*numQuadPointsPerCell + q_point]); - rhoOutHistNorm[iHist] += std::abs(d_rhoOutVals[iHist][iElem*numQuadPointsPerCell + q_point]); - } - // Compute (rhoIn-rhoOut)^2 - normValue += std::pow((rhoInValuesOld)[cell->id()][q_point] - - (*rhoOutputValues)[cell->id()][q_point], - 2.0) * - d_vecJxW[iElem*numQuadPointsPerCell + q_point]; - // Anderson mixing scheme - // double rhoOutBar=cn*(rhoOutVals[N])[cell->id()][q_point]; - // double rhoInBar=cn*(rhoInVals[N])[cell->id()][q_point]; - double rhoOutBar = d_cFinal * d_rhoOutVals[N][iElem*numQuadPointsPerCell + q_point]; - double rhoInBar = d_cFinal * d_rhoInVals[N][iElem*numQuadPointsPerCell + q_point]; - - for (int i = 0; i < N; i++) - { - rhoOutBar += d_c[i] * d_rhoOutVals[N - 1 - i][iElem*numQuadPointsPerCell + q_point]; - rhoInBar += d_c[i] * d_rhoInVals[N - 1 - i][iElem*numQuadPointsPerCell + q_point]; - } - (*rhoInputValues)[cell->id()][q_point] = - ((1 - d_dftParamsPtr->mixingParameter) * rhoInBar + - d_dftParamsPtr->mixingParameter * rhoOutBar); - rhoOutputNorm += (*rhoInputValues)[cell->id()][q_point]; - } - - if(d_dftParamsPtr->spinPolarized == 1) - { - for (unsigned int q_point = 0; q_point < d_numberQuadraturePointsPerCell; ++q_point) - { - (*rhoInValues)[cell->id()][q_point] = (*rhoInputValues)[cell->id()][2*q_point + 0] + - (*rhoInputValues)[cell->id()][2*q_point + 1]; - } - } - iElem++; - } - } - - MPI_Allreduce(MPI_IN_PLACE, - &rhoOutputNorm, - 1, - MPI_DOUBLE, - MPI_SUM, - d_mpi_comm_domain); - - MPI_Allreduce(MPI_IN_PLACE, - &rhoInHistNorm[0], - N+1, - MPI_DOUBLE, - MPI_SUM, - d_mpi_comm_domain); - MPI_Allreduce(MPI_IN_PLACE, - &rhoOutHistNorm[0], - N+1, - MPI_DOUBLE, - MPI_SUM, - d_mpi_comm_domain); - -// std::cout<<" norm of output rho = "<getDensityBasedFamilyType() == densityFamilyType::GGA) - { - if(d_dftParamsPtr->spinPolarized == 0) - { - std::vector> gradRhoOutTemp( - N + 1, std::vector(3 * d_numberQuadraturePointsPerCell, 0.0)); - - std::vector> gradRhoInTemp( - N + 1, std::vector(3 * d_numberQuadraturePointsPerCell, 0.0)); - - std::map> gradRhoInValuesOld = - *gradRhoInValues; - gradRhoInValues = std::make_shared>>(); - cell = d_dofHandler->begin_active(); - for (; cell != endc; ++cell) - { - if (cell->is_locally_owned()) - { - (*gradRhoInValues)[cell->id()] = - std::vector(3 * d_numberQuadraturePointsPerCell); - - - for (int hist = 0; hist < N + 1; hist++) - { - gradRhoOutTemp[hist] = (*(d_gradRhoOutVals[hist]))[cell->id()]; - gradRhoInTemp[hist] = (*(d_gradRhoInVals[hist]))[cell->id()]; - } - - - for (unsigned int q_point = 0; q_point < d_numberQuadraturePointsPerCell; - ++q_point) - { - for( unsigned int iHist = 0 ;iHist < N+1;iHist++) - { - gradRhoInHistNorm[iHist] += std::abs(gradRhoInTemp[iHist][3*q_point+0]); - gradRhoInHistNorm[iHist] += std::abs(gradRhoInTemp[iHist][3*q_point+1]); - gradRhoInHistNorm[iHist] += std::abs(gradRhoInTemp[iHist][3*q_point+2]); - gradRhoOutHistNorm[iHist] += std::abs(gradRhoOutTemp[iHist][3*q_point+0]); - gradRhoOutHistNorm[iHist] += std::abs(gradRhoOutTemp[iHist][3*q_point+1]); - gradRhoOutHistNorm[iHist] += std::abs(gradRhoOutTemp[iHist][3*q_point+2]); - - } - // - // Anderson mixing scheme - // - double gradRhoXOutBar = - d_cFinal * - gradRhoOutTemp - [N][3 * q_point + - 0]; // cn*(gradRhoOutVals[N])[cell->id()][3*q_point - // + 0]; - double gradRhoYOutBar = - d_cFinal * - gradRhoOutTemp - [N][3 * q_point + - 1]; // cn*(gradRhoOutVals[N])[cell->id()][3*q_point - // + 1]; - double gradRhoZOutBar = - d_cFinal * - gradRhoOutTemp - [N][3 * q_point + - 2]; // cn*(gradRhoOutVals[N])[cell->id()][3*q_point - // + 2]; - - double gradRhoXInBar = - d_cFinal * - gradRhoInTemp - [N][3 * q_point + - 0]; // cn*(gradRhoInVals[N])[cell->id()][3*q_point - // + 0]; - double gradRhoYInBar = - d_cFinal * - gradRhoInTemp - [N][3 * q_point + - 1]; // cn*(gradRhoInVals[N])[cell->id()][3*q_point - // + 1]; - double gradRhoZInBar = - d_cFinal * - gradRhoInTemp - [N][3 * q_point + - 2]; // cn*(gradRhoInVals[N])[cell->id()][3*q_point - // + 2]; + normValue += std::pow(variableHistoryOut[mixingVariableName][N][iQuad] - + variableHistoryIn[mixingVariableName][N][iQuad], + 2.0) * + vectorDotProductWeights[mixingVariableName][iQuad]; - for (int i = 0; i < N; i++) - { - gradRhoXOutBar += - d_c[i] * - gradRhoOutTemp - [N - 1 - i] - [3 * q_point + - 0]; // cTotal[i]*(gradRhoOutVals[N-1-i])[cell->id()][3*q_point - // + 0]; - gradRhoYOutBar += - d_c[i] * - gradRhoOutTemp - [N - 1 - i] - [3 * q_point + - 1]; // cTotal[i]*(gradRhoOutVals[N-1-i])[cell->id()][3*q_point - // + 1]; - gradRhoZOutBar += - d_c[i] * - gradRhoOutTemp - [N - 1 - i] - [3 * q_point + - 2]; // cTotal[i]*(gradRhoOutVals[N-1-i])[cell->id()][3*q_point - // + 2]; + double varOutBar = d_cFinal * variableHistoryOut[mixingVariableName][N][iQuad]; + double varInBar = d_cFinal * variableHistoryIn[mixingVariableName][N][iQuad]; - gradRhoXInBar += - d_c[i] * - gradRhoInTemp - [N - 1 - i] - [3 * q_point + - 0]; // cTotal[i]*(gradRhoInVals[N-1-i])[cell->id()][3*q_point - // + 0]; - gradRhoYInBar += - d_c[i] * - gradRhoInTemp - [N - 1 - i] - [3 * q_point + - 1]; // cTotal[i]*(gradRhoInVals[N-1-i])[cell->id()][3*q_point - // + 1]; - gradRhoZInBar += - d_c[i] * - gradRhoInTemp - [N - 1 - i] - [3 * q_point + - 2]; // cTotal[i]*(gradRhoInVals[N-1-i])[cell->id()][3*q_point - // + 2]; - } - - (*gradRhoInValues)[cell->id()][3 * q_point + 0] = - ((1 - d_dftParamsPtr->mixingParameter) * gradRhoXInBar + - d_dftParamsPtr->mixingParameter * gradRhoXOutBar); - (*gradRhoInValues)[cell->id()][3 * q_point + 1] = - ((1 - d_dftParamsPtr->mixingParameter) * gradRhoYInBar + - d_dftParamsPtr->mixingParameter * gradRhoYOutBar); - (*gradRhoInValues)[cell->id()][3 * q_point + 2] = - ((1 - d_dftParamsPtr->mixingParameter) * gradRhoZInBar + - d_dftParamsPtr->mixingParameter * gradRhoZOutBar); - - gradRhoOutputNorm += std::abs((*gradRhoInValues)[cell->id()][3 * q_point + 0]); - gradRhoOutputNorm += std::abs((*gradRhoInValues)[cell->id()][3 * q_point + 1]); - gradRhoOutputNorm += std::abs((*gradRhoInValues)[cell->id()][3 * q_point + 2]); - } - } - } - } - else + for (int i = 0; i < N; i++) { - std::vector> gradRhoOutTemp( - N + 1, std::vector(3 * d_numberQuadraturePointsPerCell, 0.0)); - - std::vector> gradRhoInTemp( - N + 1, std::vector(3 * d_numberQuadraturePointsPerCell, 0.0)); - - std::vector> gradRhoOutSpinPolarizedTemp( - N + 1, std::vector(6 * d_numberQuadraturePointsPerCell, 0.0)); - - std::vector> gradRhoInSpinPolarizedTemp( - N + 1, std::vector(6 * d_numberQuadraturePointsPerCell, 0.0)); - - std::map> gradRhoInValuesOld = - *gradRhoInValues; - - gradRhoInValues = std::make_shared>>(); - - gradRhoInValuesSpinPolarized = std::make_shared>>(); - cell = d_dofHandler->begin_active(); - for (; cell != endc; ++cell) - { - if (cell->is_locally_owned()) - { - (*gradRhoInValues)[cell->id()] = - std::vector(3 * d_numberQuadraturePointsPerCell); - (*gradRhoInValuesSpinPolarized)[cell->id()] = - std::vector(6 * d_numberQuadraturePointsPerCell); - // - - for (int hist = 0; hist < N + 1; hist++) - { - gradRhoOutSpinPolarizedTemp[hist] = - (*(d_gradRhoOutValsSpinPolarized[hist]))[cell->id()]; - gradRhoInSpinPolarizedTemp[hist] = - (*(d_gradRhoInValsSpinPolarized[hist]))[cell->id()]; - } - - for (unsigned int q_point = 0; q_point < d_numberQuadraturePointsPerCell; - ++q_point) - { - for( unsigned int iHist = 0 ;iHist < N+1;iHist++) - { -// gradRhoInHistNorm[iHist] += std::abs(gradRhoInTemp[iHist][3*q_point+0]); -// gradRhoInHistNorm[iHist] += std::abs(gradRhoInTemp[iHist][3*q_point+1]); -// gradRhoInHistNorm[iHist] += std::abs(gradRhoInTemp[iHist][3*q_point+2]); -// gradRhoOutHistNorm[iHist] += std::abs(gradRhoOutTemp[iHist][3*q_point+0]); -// gradRhoOutHistNorm[iHist] += std::abs(gradRhoOutTemp[iHist][3*q_point+1]); -// gradRhoOutHistNorm[iHist] += std::abs(gradRhoOutTemp[iHist][3*q_point+2]); - - gradSpinRhoInHistNorm[iHist] += std::abs(gradRhoInSpinPolarizedTemp[iHist][6*q_point+0]); - gradSpinRhoInHistNorm[iHist] += std::abs(gradRhoInSpinPolarizedTemp[iHist][6*q_point+1]); - gradSpinRhoInHistNorm[iHist] += std::abs(gradRhoInSpinPolarizedTemp[iHist][6*q_point+2]); - gradSpinRhoInHistNorm[iHist] += std::abs(gradRhoInSpinPolarizedTemp[iHist][6*q_point+3]); - gradSpinRhoInHistNorm[iHist] += std::abs(gradRhoInSpinPolarizedTemp[iHist][6*q_point+4]); - gradSpinRhoInHistNorm[iHist] += std::abs(gradRhoInSpinPolarizedTemp[iHist][6*q_point+5]); - - gradSpinRhoOutHistNorm[iHist] += std::abs(gradRhoOutSpinPolarizedTemp[iHist][6*q_point+0]); - gradSpinRhoOutHistNorm[iHist] += std::abs(gradRhoOutSpinPolarizedTemp[iHist][6*q_point+1]); - gradSpinRhoOutHistNorm[iHist] += std::abs(gradRhoOutSpinPolarizedTemp[iHist][6*q_point+2]); - gradSpinRhoOutHistNorm[iHist] += std::abs(gradRhoOutSpinPolarizedTemp[iHist][6*q_point+3]); - gradSpinRhoOutHistNorm[iHist] += std::abs(gradRhoOutSpinPolarizedTemp[iHist][6*q_point+4]); - gradSpinRhoOutHistNorm[iHist] += std::abs(gradRhoOutSpinPolarizedTemp[iHist][6*q_point+5]); - - - } - // - // Anderson mixing scheme spin up - // - double gradRhoXOutBar1 = - d_cFinal * - gradRhoOutSpinPolarizedTemp - [N] - [6 * q_point + - 0]; // cn*(gradRhoOutValsSpinPolarized[N])[cell->id()][6*q_point - // + 0]; - double gradRhoYOutBar1 = - d_cFinal * - gradRhoOutSpinPolarizedTemp - [N] - [6 * q_point + - 1]; // cn*(gradRhoOutValsSpinPolarized[N])[cell->id()][6*q_point - // + 1]; - double gradRhoZOutBar1 = - d_cFinal * - gradRhoOutSpinPolarizedTemp - [N] - [6 * q_point + - 2]; // cn*(gradRhoOutValsSpinPolarized[N])[cell->id()][6*q_point - // + 2]; - - double gradRhoXInBar1 = - d_cFinal * gradRhoInSpinPolarizedTemp[N][6 * q_point + 0]; - double gradRhoYInBar1 = - d_cFinal * gradRhoInSpinPolarizedTemp[N][6 * q_point + 1]; - double gradRhoZInBar1 = - d_cFinal * gradRhoInSpinPolarizedTemp[N][6 * q_point + 2]; - - for (int i = 0; i < N; i++) - { - gradRhoXOutBar1 += - d_c[i] * - gradRhoOutSpinPolarizedTemp[N - 1 - i] - [6 * q_point + 0]; - gradRhoYOutBar1 += - d_c[i] * - gradRhoOutSpinPolarizedTemp[N - 1 - i] - [6 * q_point + 1]; - gradRhoZOutBar1 += - d_c[i] * - gradRhoOutSpinPolarizedTemp[N - 1 - i] - [6 * q_point + 2]; - - gradRhoXInBar1 += - d_c[i] * - gradRhoInSpinPolarizedTemp[N - 1 - i] - [6 * q_point + 0]; - gradRhoYInBar1 += - d_c[i] * - gradRhoInSpinPolarizedTemp[N - 1 - i] - [6 * q_point + 1]; - gradRhoZInBar1 += - d_c[i] * - gradRhoInSpinPolarizedTemp[N - 1 - i] - [6 * q_point + 2]; - } - // - // Anderson mixing scheme spin down - // - double gradRhoXOutBar2 = - d_cFinal * gradRhoOutSpinPolarizedTemp[N][6 * q_point + 3]; - double gradRhoYOutBar2 = - d_cFinal * gradRhoOutSpinPolarizedTemp[N][6 * q_point + 4]; - double gradRhoZOutBar2 = - d_cFinal * gradRhoOutSpinPolarizedTemp[N][6 * q_point + 5]; - - double gradRhoXInBar2 = - d_cFinal * gradRhoInSpinPolarizedTemp[N][6 * q_point + 3]; - double gradRhoYInBar2 = - d_cFinal * gradRhoInSpinPolarizedTemp[N][6 * q_point + 4]; - double gradRhoZInBar2 = - d_cFinal * gradRhoInSpinPolarizedTemp[N][6 * q_point + 5]; - - for (int i = 0; i < N; i++) - { - gradRhoXOutBar2 += - d_c[i] * - gradRhoOutSpinPolarizedTemp[N - 1 - i] - [6 * q_point + 3]; - gradRhoYOutBar2 += - d_c[i] * - gradRhoOutSpinPolarizedTemp[N - 1 - i] - [6 * q_point + 4]; - gradRhoZOutBar2 += - d_c[i] * - gradRhoOutSpinPolarizedTemp[N - 1 - i] - [6 * q_point + 5]; - - gradRhoXInBar2 += - d_c[i] * - gradRhoInSpinPolarizedTemp[N - 1 - i] - [6 * q_point + 3]; - gradRhoYInBar2 += - d_c[i] * - gradRhoInSpinPolarizedTemp[N - 1 - i] - [6 * q_point + 4]; - gradRhoZInBar2 += - d_c[i] * - gradRhoInSpinPolarizedTemp[N - 1 - i] - [6 * q_point + 5]; - } - // - (*gradRhoInValuesSpinPolarized)[cell->id()][6 * q_point + - 0] = - ((1 - d_dftParamsPtr->mixingParameter) * gradRhoXInBar1 + - d_dftParamsPtr->mixingParameter * gradRhoXOutBar1); - (*gradRhoInValuesSpinPolarized)[cell->id()][6 * q_point + - 1] = - ((1 - d_dftParamsPtr->mixingParameter) * gradRhoYInBar1 + - d_dftParamsPtr->mixingParameter * gradRhoYOutBar1); - (*gradRhoInValuesSpinPolarized)[cell->id()][6 * q_point + - 2] = - ((1 - d_dftParamsPtr->mixingParameter) * gradRhoZInBar1 + - d_dftParamsPtr->mixingParameter * gradRhoZOutBar1); - (*gradRhoInValuesSpinPolarized)[cell->id()][6 * q_point + - 3] = - ((1 - d_dftParamsPtr->mixingParameter) * gradRhoXInBar2 + - d_dftParamsPtr->mixingParameter * gradRhoXOutBar2); - (*gradRhoInValuesSpinPolarized)[cell->id()][6 * q_point + - 4] = - ((1 - d_dftParamsPtr->mixingParameter) * gradRhoYInBar2 + - d_dftParamsPtr->mixingParameter * gradRhoYOutBar2); - (*gradRhoInValuesSpinPolarized)[cell->id()][6 * q_point + - 5] = - ((1 - d_dftParamsPtr->mixingParameter) * gradRhoZInBar2 + - d_dftParamsPtr->mixingParameter * gradRhoZOutBar2); - - ((*gradRhoInValues)[cell->id()][3 * q_point + 0]) = - ((*gradRhoInValuesSpinPolarized)[cell->id()] - [6 * q_point + 0]) + - ((*gradRhoInValuesSpinPolarized)[cell->id()] - [6 * q_point + 3]); - ((*gradRhoInValues)[cell->id()][3 * q_point + 1]) = - ((*gradRhoInValuesSpinPolarized)[cell->id()] - [6 * q_point + 1]) + - ((*gradRhoInValuesSpinPolarized)[cell->id()] - [6 * q_point + 4]); - ((*gradRhoInValues)[cell->id()][3 * q_point + 2]) = - ((*gradRhoInValuesSpinPolarized)[cell->id()] - [6 * q_point + 2]) + - ((*gradRhoInValuesSpinPolarized)[cell->id()] - [6 * q_point + 5]); - - gradRhoOutputNorm += std::abs((*gradRhoInValuesSpinPolarized)[cell->id()][6 * q_point + 0]); - gradRhoOutputNorm += std::abs((*gradRhoInValuesSpinPolarized)[cell->id()][6 * q_point + 1]); - gradRhoOutputNorm += std::abs((*gradRhoInValuesSpinPolarized)[cell->id()][6 * q_point + 2]); - gradRhoOutputNorm += std::abs((*gradRhoInValuesSpinPolarized)[cell->id()][6 * q_point + 3]); - gradRhoOutputNorm += std::abs((*gradRhoInValuesSpinPolarized)[cell->id()][6 * q_point + 4]); - gradRhoOutputNorm += std::abs((*gradRhoInValuesSpinPolarized)[cell->id()][6 * q_point + 5]); - } - } - } + varOutBar += d_c[i] * variableHistoryOut[mixingVariableName][N - 1 - i][iQuad]; + varInBar += d_c[i] * variableHistoryIn[mixingVariableName][N - 1 - i][iQuad]; } - } - - MPI_Allreduce(MPI_IN_PLACE, - &gradRhoOutputNorm, - 1, - MPI_DOUBLE, - MPI_SUM, - d_mpi_comm_domain); - - MPI_Allreduce(MPI_IN_PLACE, - &gradRhoInHistNorm[0], - N+1, - MPI_DOUBLE, - MPI_SUM, - d_mpi_comm_domain); - MPI_Allreduce(MPI_IN_PLACE, - &gradRhoOutHistNorm[0], - N+1, - MPI_DOUBLE, - MPI_SUM, - d_mpi_comm_domain); - - MPI_Allreduce(MPI_IN_PLACE, - &gradSpinRhoInHistNorm[0], - N+1, - MPI_DOUBLE, - MPI_SUM, - d_mpi_comm_domain); - MPI_Allreduce(MPI_IN_PLACE, - &gradSpinRhoOutHistNorm[0], - N+1, - MPI_DOUBLE, - MPI_SUM, - d_mpi_comm_domain); - -// std::cout<<" norm of output grad rho = "<= d_dftParamsPtr->mixingHistory) + if (variableHistoryIn[mixingVariable::rho].size() >= d_dftParamsPtr->mixingHistory) { - d_rhoInVals.pop_front(); - d_rhoOutVals.pop_front(); - - if (d_excManagerPtr->getDensityBasedFamilyType() == - densityFamilyType::GGA) // GGA - { - d_gradRhoInVals.pop_front(); - d_gradRhoOutVals.pop_front(); - } - if (d_dftParamsPtr->spinPolarized == 1 && - d_excManagerPtr->getDensityBasedFamilyType() == - densityFamilyType::GGA) + for (const auto& [key, value] : variableHistoryIn) { - d_gradRhoInValsSpinPolarized.pop_front(); - d_gradRhoOutValsSpinPolarized.pop_front(); + variableHistoryIn[key].pop_front(); + variableHistoryOut[key].pop_front(); } - } } - void MixingScheme::popRhoInHist() - { - d_rhoInVals.pop_front(); - } - - void MixingScheme::popRhoOutHist() - { - d_rhoOutVals.pop_front(); - } - - void MixingScheme::popGradRhoInHist() - { - d_gradRhoInVals.pop_front(); - } - - void MixingScheme::popGradRhoOutHist() - { - d_gradRhoOutVals.pop_front(); - } - - void MixingScheme::popGradRhoSpinInHist() - { - d_gradRhoInValsSpinPolarized.pop_front(); - } +} // end of namespace - void MixingScheme::popGradRhoSpinOutHist() - { - d_gradRhoOutValsSpinPolarized.pop_front(); - } -} // end of namespace +// void MixingScheme::popRhoInHist() +// { +// d_rhoInVals.pop_front(); +// } +// +// void MixingScheme::popRhoOutHist() +// { +// d_rhoOutVals.pop_front(); +// } +// +// void MixingScheme::popGradRhoInHist() +// { +// d_gradRhoInVals.pop_front(); +// } +// +// void MixingScheme::popGradRhoOutHist() +// { +// d_gradRhoOutVals.pop_front(); +// } +// +// void MixingScheme::popGradRhoSpinInHist() +// { +// d_gradRhoInValsSpinPolarized.pop_front(); +// } +// +// void MixingScheme::popGradRhoSpinOutHist() +// { +// d_gradRhoOutValsSpinPolarized.pop_front(); +// } + +//double MixingScheme::mixDensity(std::shared_ptr>> &rhoInValues, +// std::shared_ptr>> &rhoOutValues, +// std::shared_ptr>> &rhoInValuesSpinPolarized, +// std::shared_ptr>> &rhoOutValuesSpinPolarized, +// std::shared_ptr>> &gradRhoInValues, +// std::shared_ptr>> &gradRhoOutValues, +// std::shared_ptr>> &gradRhoInValuesSpinPolarized, +// std::shared_ptr>> &gradRhoOutValuesSpinPolarized) +//{ +// double normValue = 0.0; +// +// // std::cout<<" size of rhoInHist = "< rhoInHistNorm(N+1,0.0), rhoOutHistNorm(N+1,0.0), gradRhoInHistNorm(N+1,0.0), gradRhoOutHistNorm(N+1,0.0), gradSpinRhoInHistNorm(N+1,0.0), gradSpinRhoOutHistNorm(N+1,0.0); +// double rhoOutputNorm = 0.0, gradRhoOutputNorm = 0.0; +// +// std::fill(rhoInHistNorm.begin(),rhoInHistNorm.end(),0.0); +// std::fill(rhoOutHistNorm.begin(),rhoOutHistNorm.end(),0.0); +// std::fill(gradRhoInHistNorm.begin(),gradRhoInHistNorm.end(),0.0); +// std::fill(gradRhoOutHistNorm.begin(),gradRhoOutHistNorm.end(),0.0); +// std::fill(gradSpinRhoInHistNorm.begin(),gradSpinRhoInHistNorm.end(),0.0); +// std::fill(gradSpinRhoOutHistNorm.begin(),gradSpinRhoOutHistNorm.end(),0.0); +// +// std::shared_ptr>> rhoInputValues, rhoOutputValues; +// // create new rhoValue tables +// std::map> rhoInValuesOld = *rhoInValues; +// +// +// unsigned int numQuadPointsPerCell = d_numberQuadraturePointsPerCell; +// if(d_dftParamsPtr->spinPolarized == 0) +// { +// rhoInValuesOld = *rhoInValues; +// rhoInputValues = rhoInValues; +// rhoOutputValues = rhoOutValues; +// } +// if (d_dftParamsPtr->spinPolarized == 1) +// { +// numQuadPointsPerCell = 2*numQuadPointsPerCell; +// +// rhoInValuesOld = *rhoInValuesSpinPolarized; +// rhoInputValues = rhoInValuesSpinPolarized; +// rhoOutputValues = rhoOutValuesSpinPolarized; +// } +// +// // implement anderson mixing +// typename dealii::DoFHandler<3>::active_cell_iterator +// cell = d_dofHandler->begin_active(), +// endc = d_dofHandler->end(); +// unsigned int iElem = 0; +// for (; cell != endc; ++cell) +// { +// if (cell->is_locally_owned()) +// { +// (*rhoInputValues)[cell->id()].resize(numQuadPointsPerCell); +// +// for (unsigned int q_point = 0; q_point < numQuadPointsPerCell; ++q_point) +// { +// +// for(unsigned int iHist = 0; iHist < N+1; iHist++) +// { +// rhoInHistNorm[iHist] += std::abs(d_rhoInVals[iHist][iElem*numQuadPointsPerCell + q_point]); +// rhoOutHistNorm[iHist] += std::abs(d_rhoOutVals[iHist][iElem*numQuadPointsPerCell + q_point]); +// } +// // Compute (rhoIn-rhoOut)^2 +// normValue += std::pow((rhoInValuesOld)[cell->id()][q_point] - +// (*rhoOutputValues)[cell->id()][q_point], +// 2.0) * +// d_vecJxW[iElem*numQuadPointsPerCell + q_point]; +// // Anderson mixing scheme +// // double rhoOutBar=cn*(rhoOutVals[N])[cell->id()][q_point]; +// // double rhoInBar=cn*(rhoInVals[N])[cell->id()][q_point]; +// double rhoOutBar = d_cFinal * d_rhoOutVals[N][iElem*numQuadPointsPerCell + q_point]; +// double rhoInBar = d_cFinal * d_rhoInVals[N][iElem*numQuadPointsPerCell + q_point]; +// +// for (int i = 0; i < N; i++) +// { +// rhoOutBar += d_c[i] * d_rhoOutVals[N - 1 - i][iElem*numQuadPointsPerCell + q_point]; +// rhoInBar += d_c[i] * d_rhoInVals[N - 1 - i][iElem*numQuadPointsPerCell + q_point]; +// } +// (*rhoInputValues)[cell->id()][q_point] = +// ((1 - d_dftParamsPtr->mixingParameter) * rhoInBar + +// d_dftParamsPtr->mixingParameter * rhoOutBar); +// rhoOutputNorm += (*rhoInputValues)[cell->id()][q_point]; +// } +// +// if(d_dftParamsPtr->spinPolarized == 1) +// { +// for (unsigned int q_point = 0; q_point < d_numberQuadraturePointsPerCell; ++q_point) +// { +// (*rhoInValues)[cell->id()][q_point] = (*rhoInputValues)[cell->id()][2*q_point + 0] + +// (*rhoInputValues)[cell->id()][2*q_point + 1]; +// } +// } +// iElem++; +// } +// } +// +// MPI_Allreduce(MPI_IN_PLACE, +// &rhoOutputNorm, +// 1, +// MPI_DOUBLE, +// MPI_SUM, +// d_mpi_comm_domain); +// +// MPI_Allreduce(MPI_IN_PLACE, +// &rhoInHistNorm[0], +// N+1, +// MPI_DOUBLE, +// MPI_SUM, +// d_mpi_comm_domain); +// MPI_Allreduce(MPI_IN_PLACE, +// &rhoOutHistNorm[0], +// N+1, +// MPI_DOUBLE, +// MPI_SUM, +// d_mpi_comm_domain); +// +// // std::cout<<" norm of output rho = "<getDensityBasedFamilyType() == densityFamilyType::GGA) +// { +// if(d_dftParamsPtr->spinPolarized == 0) +// { +// std::vector> gradRhoOutTemp( +// N + 1, std::vector(3 * d_numberQuadraturePointsPerCell, 0.0)); +// +// std::vector> gradRhoInTemp( +// N + 1, std::vector(3 * d_numberQuadraturePointsPerCell, 0.0)); +// +// std::map> gradRhoInValuesOld = +// *gradRhoInValues; +// gradRhoInValues = std::make_shared>>(); +// cell = d_dofHandler->begin_active(); +// for (; cell != endc; ++cell) +// { +// if (cell->is_locally_owned()) +// { +// (*gradRhoInValues)[cell->id()] = +// std::vector(3 * d_numberQuadraturePointsPerCell); +// +// +// for (int hist = 0; hist < N + 1; hist++) +// { +// gradRhoOutTemp[hist] = (*(d_gradRhoOutVals[hist]))[cell->id()]; +// gradRhoInTemp[hist] = (*(d_gradRhoInVals[hist]))[cell->id()]; +// } +// +// +// for (unsigned int q_point = 0; q_point < d_numberQuadraturePointsPerCell; +// ++q_point) +// { +// for( unsigned int iHist = 0 ;iHist < N+1;iHist++) +// { +// gradRhoInHistNorm[iHist] += std::abs(gradRhoInTemp[iHist][3*q_point+0]); +// gradRhoInHistNorm[iHist] += std::abs(gradRhoInTemp[iHist][3*q_point+1]); +// gradRhoInHistNorm[iHist] += std::abs(gradRhoInTemp[iHist][3*q_point+2]); +// gradRhoOutHistNorm[iHist] += std::abs(gradRhoOutTemp[iHist][3*q_point+0]); +// gradRhoOutHistNorm[iHist] += std::abs(gradRhoOutTemp[iHist][3*q_point+1]); +// gradRhoOutHistNorm[iHist] += std::abs(gradRhoOutTemp[iHist][3*q_point+2]); +// +// } +// // +// // Anderson mixing scheme +// // +// double gradRhoXOutBar = +// d_cFinal * +// gradRhoOutTemp +// [N][3 * q_point + +// 0]; // cn*(gradRhoOutVals[N])[cell->id()][3*q_point +// // + 0]; +// double gradRhoYOutBar = +// d_cFinal * +// gradRhoOutTemp +// [N][3 * q_point + +// 1]; // cn*(gradRhoOutVals[N])[cell->id()][3*q_point +// // + 1]; +// double gradRhoZOutBar = +// d_cFinal * +// gradRhoOutTemp +// [N][3 * q_point + +// 2]; // cn*(gradRhoOutVals[N])[cell->id()][3*q_point +// // + 2]; +// +// double gradRhoXInBar = +// d_cFinal * +// gradRhoInTemp +// [N][3 * q_point + +// 0]; // cn*(gradRhoInVals[N])[cell->id()][3*q_point +// // + 0]; +// double gradRhoYInBar = +// d_cFinal * +// gradRhoInTemp +// [N][3 * q_point + +// 1]; // cn*(gradRhoInVals[N])[cell->id()][3*q_point +// // + 1]; +// double gradRhoZInBar = +// d_cFinal * +// gradRhoInTemp +// [N][3 * q_point + +// 2]; // cn*(gradRhoInVals[N])[cell->id()][3*q_point +// // + 2]; +// +// for (int i = 0; i < N; i++) +// { +// gradRhoXOutBar += +// d_c[i] * +// gradRhoOutTemp +// [N - 1 - i] +// [3 * q_point + +// 0]; // cTotal[i]*(gradRhoOutVals[N-1-i])[cell->id()][3*q_point +// // + 0]; +// gradRhoYOutBar += +// d_c[i] * +// gradRhoOutTemp +// [N - 1 - i] +// [3 * q_point + +// 1]; // cTotal[i]*(gradRhoOutVals[N-1-i])[cell->id()][3*q_point +// // + 1]; +// gradRhoZOutBar += +// d_c[i] * +// gradRhoOutTemp +// [N - 1 - i] +// [3 * q_point + +// 2]; // cTotal[i]*(gradRhoOutVals[N-1-i])[cell->id()][3*q_point +// // + 2]; +// +// gradRhoXInBar += +// d_c[i] * +// gradRhoInTemp +// [N - 1 - i] +// [3 * q_point + +// 0]; // cTotal[i]*(gradRhoInVals[N-1-i])[cell->id()][3*q_point +// // + 0]; +// gradRhoYInBar += +// d_c[i] * +// gradRhoInTemp +// [N - 1 - i] +// [3 * q_point + +// 1]; // cTotal[i]*(gradRhoInVals[N-1-i])[cell->id()][3*q_point +// // + 1]; +// gradRhoZInBar += +// d_c[i] * +// gradRhoInTemp +// [N - 1 - i] +// [3 * q_point + +// 2]; // cTotal[i]*(gradRhoInVals[N-1-i])[cell->id()][3*q_point +// // + 2]; +// } +// +// (*gradRhoInValues)[cell->id()][3 * q_point + 0] = +// ((1 - d_dftParamsPtr->mixingParameter) * gradRhoXInBar + +// d_dftParamsPtr->mixingParameter * gradRhoXOutBar); +// (*gradRhoInValues)[cell->id()][3 * q_point + 1] = +// ((1 - d_dftParamsPtr->mixingParameter) * gradRhoYInBar + +// d_dftParamsPtr->mixingParameter * gradRhoYOutBar); +// (*gradRhoInValues)[cell->id()][3 * q_point + 2] = +// ((1 - d_dftParamsPtr->mixingParameter) * gradRhoZInBar + +// d_dftParamsPtr->mixingParameter * gradRhoZOutBar); +// +// gradRhoOutputNorm += std::abs((*gradRhoInValues)[cell->id()][3 * q_point + 0]); +// gradRhoOutputNorm += std::abs((*gradRhoInValues)[cell->id()][3 * q_point + 1]); +// gradRhoOutputNorm += std::abs((*gradRhoInValues)[cell->id()][3 * q_point + 2]); +// } +// } +// } +// } +// else +// { +// std::vector> gradRhoOutTemp( +// N + 1, std::vector(3 * d_numberQuadraturePointsPerCell, 0.0)); +// +// std::vector> gradRhoInTemp( +// N + 1, std::vector(3 * d_numberQuadraturePointsPerCell, 0.0)); +// +// std::vector> gradRhoOutSpinPolarizedTemp( +// N + 1, std::vector(6 * d_numberQuadraturePointsPerCell, 0.0)); +// +// std::vector> gradRhoInSpinPolarizedTemp( +// N + 1, std::vector(6 * d_numberQuadraturePointsPerCell, 0.0)); +// +// std::map> gradRhoInValuesOld = +// *gradRhoInValues; +// +// gradRhoInValues = std::make_shared>>(); +// +// gradRhoInValuesSpinPolarized = std::make_shared>>(); +// cell = d_dofHandler->begin_active(); +// for (; cell != endc; ++cell) +// { +// if (cell->is_locally_owned()) +// { +// (*gradRhoInValues)[cell->id()] = +// std::vector(3 * d_numberQuadraturePointsPerCell); +// (*gradRhoInValuesSpinPolarized)[cell->id()] = +// std::vector(6 * d_numberQuadraturePointsPerCell); +// // +// +// for (int hist = 0; hist < N + 1; hist++) +// { +// gradRhoOutSpinPolarizedTemp[hist] = +// (*(d_gradRhoOutValsSpinPolarized[hist]))[cell->id()]; +// gradRhoInSpinPolarizedTemp[hist] = +// (*(d_gradRhoInValsSpinPolarized[hist]))[cell->id()]; +// } +// +// for (unsigned int q_point = 0; q_point < d_numberQuadraturePointsPerCell; +// ++q_point) +// { +// for( unsigned int iHist = 0 ;iHist < N+1;iHist++) +// { +// // gradRhoInHistNorm[iHist] += std::abs(gradRhoInTemp[iHist][3*q_point+0]); +// // gradRhoInHistNorm[iHist] += std::abs(gradRhoInTemp[iHist][3*q_point+1]); +// // gradRhoInHistNorm[iHist] += std::abs(gradRhoInTemp[iHist][3*q_point+2]); +// // gradRhoOutHistNorm[iHist] += std::abs(gradRhoOutTemp[iHist][3*q_point+0]); +// // gradRhoOutHistNorm[iHist] += std::abs(gradRhoOutTemp[iHist][3*q_point+1]); +// // gradRhoOutHistNorm[iHist] += std::abs(gradRhoOutTemp[iHist][3*q_point+2]); +// +// gradSpinRhoInHistNorm[iHist] += std::abs(gradRhoInSpinPolarizedTemp[iHist][6*q_point+0]); +// gradSpinRhoInHistNorm[iHist] += std::abs(gradRhoInSpinPolarizedTemp[iHist][6*q_point+1]); +// gradSpinRhoInHistNorm[iHist] += std::abs(gradRhoInSpinPolarizedTemp[iHist][6*q_point+2]); +// gradSpinRhoInHistNorm[iHist] += std::abs(gradRhoInSpinPolarizedTemp[iHist][6*q_point+3]); +// gradSpinRhoInHistNorm[iHist] += std::abs(gradRhoInSpinPolarizedTemp[iHist][6*q_point+4]); +// gradSpinRhoInHistNorm[iHist] += std::abs(gradRhoInSpinPolarizedTemp[iHist][6*q_point+5]); +// +// gradSpinRhoOutHistNorm[iHist] += std::abs(gradRhoOutSpinPolarizedTemp[iHist][6*q_point+0]); +// gradSpinRhoOutHistNorm[iHist] += std::abs(gradRhoOutSpinPolarizedTemp[iHist][6*q_point+1]); +// gradSpinRhoOutHistNorm[iHist] += std::abs(gradRhoOutSpinPolarizedTemp[iHist][6*q_point+2]); +// gradSpinRhoOutHistNorm[iHist] += std::abs(gradRhoOutSpinPolarizedTemp[iHist][6*q_point+3]); +// gradSpinRhoOutHistNorm[iHist] += std::abs(gradRhoOutSpinPolarizedTemp[iHist][6*q_point+4]); +// gradSpinRhoOutHistNorm[iHist] += std::abs(gradRhoOutSpinPolarizedTemp[iHist][6*q_point+5]); +// +// +// } +// // +// // Anderson mixing scheme spin up +// // +// double gradRhoXOutBar1 = +// d_cFinal * +// gradRhoOutSpinPolarizedTemp +// [N] +// [6 * q_point + +// 0]; // cn*(gradRhoOutValsSpinPolarized[N])[cell->id()][6*q_point +// // + 0]; +// double gradRhoYOutBar1 = +// d_cFinal * +// gradRhoOutSpinPolarizedTemp +// [N] +// [6 * q_point + +// 1]; // cn*(gradRhoOutValsSpinPolarized[N])[cell->id()][6*q_point +// // + 1]; +// double gradRhoZOutBar1 = +// d_cFinal * +// gradRhoOutSpinPolarizedTemp +// [N] +// [6 * q_point + +// 2]; // cn*(gradRhoOutValsSpinPolarized[N])[cell->id()][6*q_point +// // + 2]; +// +// double gradRhoXInBar1 = +// d_cFinal * gradRhoInSpinPolarizedTemp[N][6 * q_point + 0]; +// double gradRhoYInBar1 = +// d_cFinal * gradRhoInSpinPolarizedTemp[N][6 * q_point + 1]; +// double gradRhoZInBar1 = +// d_cFinal * gradRhoInSpinPolarizedTemp[N][6 * q_point + 2]; +// +// for (int i = 0; i < N; i++) +// { +// gradRhoXOutBar1 += +// d_c[i] * +// gradRhoOutSpinPolarizedTemp[N - 1 - i] +// [6 * q_point + 0]; +// gradRhoYOutBar1 += +// d_c[i] * +// gradRhoOutSpinPolarizedTemp[N - 1 - i] +// [6 * q_point + 1]; +// gradRhoZOutBar1 += +// d_c[i] * +// gradRhoOutSpinPolarizedTemp[N - 1 - i] +// [6 * q_point + 2]; +// +// gradRhoXInBar1 += +// d_c[i] * +// gradRhoInSpinPolarizedTemp[N - 1 - i] +// [6 * q_point + 0]; +// gradRhoYInBar1 += +// d_c[i] * +// gradRhoInSpinPolarizedTemp[N - 1 - i] +// [6 * q_point + 1]; +// gradRhoZInBar1 += +// d_c[i] * +// gradRhoInSpinPolarizedTemp[N - 1 - i] +// [6 * q_point + 2]; +// } +// // +// // Anderson mixing scheme spin down +// // +// double gradRhoXOutBar2 = +// d_cFinal * gradRhoOutSpinPolarizedTemp[N][6 * q_point + 3]; +// double gradRhoYOutBar2 = +// d_cFinal * gradRhoOutSpinPolarizedTemp[N][6 * q_point + 4]; +// double gradRhoZOutBar2 = +// d_cFinal * gradRhoOutSpinPolarizedTemp[N][6 * q_point + 5]; +// +// double gradRhoXInBar2 = +// d_cFinal * gradRhoInSpinPolarizedTemp[N][6 * q_point + 3]; +// double gradRhoYInBar2 = +// d_cFinal * gradRhoInSpinPolarizedTemp[N][6 * q_point + 4]; +// double gradRhoZInBar2 = +// d_cFinal * gradRhoInSpinPolarizedTemp[N][6 * q_point + 5]; +// +// for (int i = 0; i < N; i++) +// { +// gradRhoXOutBar2 += +// d_c[i] * +// gradRhoOutSpinPolarizedTemp[N - 1 - i] +// [6 * q_point + 3]; +// gradRhoYOutBar2 += +// d_c[i] * +// gradRhoOutSpinPolarizedTemp[N - 1 - i] +// [6 * q_point + 4]; +// gradRhoZOutBar2 += +// d_c[i] * +// gradRhoOutSpinPolarizedTemp[N - 1 - i] +// [6 * q_point + 5]; +// +// gradRhoXInBar2 += +// d_c[i] * +// gradRhoInSpinPolarizedTemp[N - 1 - i] +// [6 * q_point + 3]; +// gradRhoYInBar2 += +// d_c[i] * +// gradRhoInSpinPolarizedTemp[N - 1 - i] +// [6 * q_point + 4]; +// gradRhoZInBar2 += +// d_c[i] * +// gradRhoInSpinPolarizedTemp[N - 1 - i] +// [6 * q_point + 5]; +// } +// // +// (*gradRhoInValuesSpinPolarized)[cell->id()][6 * q_point + +// 0] = +// ((1 - d_dftParamsPtr->mixingParameter) * gradRhoXInBar1 + +// d_dftParamsPtr->mixingParameter * gradRhoXOutBar1); +// (*gradRhoInValuesSpinPolarized)[cell->id()][6 * q_point + +// 1] = +// ((1 - d_dftParamsPtr->mixingParameter) * gradRhoYInBar1 + +// d_dftParamsPtr->mixingParameter * gradRhoYOutBar1); +// (*gradRhoInValuesSpinPolarized)[cell->id()][6 * q_point + +// 2] = +// ((1 - d_dftParamsPtr->mixingParameter) * gradRhoZInBar1 + +// d_dftParamsPtr->mixingParameter * gradRhoZOutBar1); +// (*gradRhoInValuesSpinPolarized)[cell->id()][6 * q_point + +// 3] = +// ((1 - d_dftParamsPtr->mixingParameter) * gradRhoXInBar2 + +// d_dftParamsPtr->mixingParameter * gradRhoXOutBar2); +// (*gradRhoInValuesSpinPolarized)[cell->id()][6 * q_point + +// 4] = +// ((1 - d_dftParamsPtr->mixingParameter) * gradRhoYInBar2 + +// d_dftParamsPtr->mixingParameter * gradRhoYOutBar2); +// (*gradRhoInValuesSpinPolarized)[cell->id()][6 * q_point + +// 5] = +// ((1 - d_dftParamsPtr->mixingParameter) * gradRhoZInBar2 + +// d_dftParamsPtr->mixingParameter * gradRhoZOutBar2); +// +// ((*gradRhoInValues)[cell->id()][3 * q_point + 0]) = +// ((*gradRhoInValuesSpinPolarized)[cell->id()] +// [6 * q_point + 0]) + +// ((*gradRhoInValuesSpinPolarized)[cell->id()] +// [6 * q_point + 3]); +// ((*gradRhoInValues)[cell->id()][3 * q_point + 1]) = +// ((*gradRhoInValuesSpinPolarized)[cell->id()] +// [6 * q_point + 1]) + +// ((*gradRhoInValuesSpinPolarized)[cell->id()] +// [6 * q_point + 4]); +// ((*gradRhoInValues)[cell->id()][3 * q_point + 2]) = +// ((*gradRhoInValuesSpinPolarized)[cell->id()] +// [6 * q_point + 2]) + +// ((*gradRhoInValuesSpinPolarized)[cell->id()] +// [6 * q_point + 5]); +// +// gradRhoOutputNorm += std::abs((*gradRhoInValuesSpinPolarized)[cell->id()][6 * q_point + 0]); +// gradRhoOutputNorm += std::abs((*gradRhoInValuesSpinPolarized)[cell->id()][6 * q_point + 1]); +// gradRhoOutputNorm += std::abs((*gradRhoInValuesSpinPolarized)[cell->id()][6 * q_point + 2]); +// gradRhoOutputNorm += std::abs((*gradRhoInValuesSpinPolarized)[cell->id()][6 * q_point + 3]); +// gradRhoOutputNorm += std::abs((*gradRhoInValuesSpinPolarized)[cell->id()][6 * q_point + 4]); +// gradRhoOutputNorm += std::abs((*gradRhoInValuesSpinPolarized)[cell->id()][6 * q_point + 5]); +// } +// } +// } +// } +// } +// +// MPI_Allreduce(MPI_IN_PLACE, +// &gradRhoOutputNorm, +// 1, +// MPI_DOUBLE, +// MPI_SUM, +// d_mpi_comm_domain); +// +// MPI_Allreduce(MPI_IN_PLACE, +// &gradRhoInHistNorm[0], +// N+1, +// MPI_DOUBLE, +// MPI_SUM, +// d_mpi_comm_domain); +// MPI_Allreduce(MPI_IN_PLACE, +// &gradRhoOutHistNorm[0], +// N+1, +// MPI_DOUBLE, +// MPI_SUM, +// d_mpi_comm_domain); +// +// MPI_Allreduce(MPI_IN_PLACE, +// &gradSpinRhoInHistNorm[0], +// N+1, +// MPI_DOUBLE, +// MPI_SUM, +// d_mpi_comm_domain); +// MPI_Allreduce(MPI_IN_PLACE, +// &gradSpinRhoOutHistNorm[0], +// N+1, +// MPI_DOUBLE, +// MPI_SUM, +// d_mpi_comm_domain); +// +// // std::cout<<" norm of output grad rho = "< & matrixFreeData, +// const unsigned int matrixFreeVectorComponent, +// const unsigned int matrixFreeQuadratureComponent, +// excManager * excManagerPtr) +//{ +// +// d_matrixFreeDataPtr = &matrixFreeData ; +// d_matrixFreeVectorComponent = matrixFreeVectorComponent; +// d_matrixFreeQuadratureComponent = matrixFreeQuadratureComponent; +// d_excManagerPtr = excManagerPtr; +// const dealii::Quadrature<3> &quadratureRhs = +// d_matrixFreeDataPtr->get_quadrature( +// d_matrixFreeQuadratureComponent); +// +// d_numberQuadraturePointsPerCell = quadratureRhs.size(); +// +// d_dofHandler = +// &d_matrixFreeDataPtr->get_dof_handler(d_matrixFreeVectorComponent); +// +// dealii::FEValues<3> fe_values(d_dofHandler->get_fe(), +// quadratureRhs, +// dealii::update_JxW_values); +// +// unsigned int iElem = 0; +// +// d_numCells = d_matrixFreeDataPtr->n_physical_cells(); +// +// if(d_dftParamsPtr->spinPolarized == 0) +// { +// d_vecJxW.resize(d_numberQuadraturePointsPerCell* +// d_numCells); +// typename dealii::DoFHandler<3>::active_cell_iterator +// cell = d_dofHandler->begin_active(), +// endc = d_dofHandler->end(); +// for (; cell != endc; ++cell) +// { +// if (cell->is_locally_owned()) +// { +// fe_values.reinit(cell); +// for(unsigned int iQuad = 0 ;iQuad < d_numberQuadraturePointsPerCell; iQuad++) +// { +// d_vecJxW[iElem*d_numberQuadraturePointsPerCell + +// iQuad] = +// fe_values.JxW(iQuad); +// } +// iElem++; +// } +// } +// } +// else +// { +// d_vecJxW.resize(2* +// d_numberQuadraturePointsPerCell* +// d_numCells,0.0); +// typename dealii::DoFHandler<3>::active_cell_iterator +// cell = d_dofHandler->begin_active(), +// endc = d_dofHandler->end(); +// for (; cell != endc; ++cell) +// { +// if (cell->is_locally_owned()) +// { +// fe_values.reinit(cell); +// for(unsigned int iQuad = 0 ;iQuad < d_numberQuadraturePointsPerCell; iQuad++) +// { +// d_vecJxW[2*iElem*d_numberQuadraturePointsPerCell + +// 2*iQuad +// +0] = +// fe_values.JxW(iQuad); +// +// d_vecJxW[2*iElem*d_numberQuadraturePointsPerCell + +// 2*iQuad +// +1] = +// fe_values.JxW(iQuad); +// } +// iElem++; +// } +// } +// } +// // std::cout<<" size of JxW = "<>> &gradInput) +//{ +// gradInput = d_gradRhoInVals.back(); +//} +// +//void MixingScheme::copyGradDensityFromOutHist(std::shared_ptr>> &gradOutput) +//{ +// gradOutput = d_gradRhoOutVals.back(); +//} +// +//void MixingScheme::copySpinGradDensityFromInHist(std::shared_ptr>> &gradInputSpin) +//{ +// gradInputSpin = d_gradRhoInValsSpinPolarized.back(); +//} +// +//void MixingScheme::copySpinGradDensityFromOutHist(std::shared_ptr>> &gradOutputSpin) +//{ +// gradOutputSpin = d_gradRhoOutValsSpinPolarized.back(); +//} +// +//void MixingScheme::copyGradDensityToInHist(std::shared_ptr>> &gradInput) +//{ +// d_gradRhoInVals.push_back(gradInput); +//} +//void MixingScheme::copyGradDensityToOutHist(std::shared_ptr>> &gradOutput) +//{ +// d_gradRhoOutVals.push_back(gradOutput); +//} +// +//void MixingScheme::copySpinGradDensityToInHist(std::shared_ptr>> &gradInputSpin) +//{ +// d_gradRhoInValsSpinPolarized.push_back(gradInputSpin); +//} +//void MixingScheme::copySpinGradDensityToOutHist(std::shared_ptr>> &gradOutputSpin) +//{ +// d_gradRhoOutValsSpinPolarized.push_back(gradOutputSpin); +//} + +//void MixingScheme::copyDensityFromInHist(std::shared_ptr>> &rhoInValues) +//{ +// unsigned int numQuadPointsPerCell = d_numberQuadraturePointsPerCell; +// if (d_dftParamsPtr->spinPolarized == 1) +// { +// numQuadPointsPerCell = 2*numQuadPointsPerCell; +// } +// +// unsigned int hist =d_rhoInVals.size(); +// +// unsigned int iElem = 0; +// +// typename dealii::DoFHandler<3>::active_cell_iterator +// cell = d_dofHandler->begin_active(), +// endc = d_dofHandler->end(); +// for (; cell != endc; ++cell) +// { +// if (cell->is_locally_owned()) +// { +// for(unsigned int iQuad = 0 ;iQuad < numQuadPointsPerCell; iQuad++) +// { +// (*rhoInValues)[cell->id()][iQuad] = +// d_rhoInVals[hist-1][iElem*numQuadPointsPerCell + +// iQuad]; +// } +// iElem++; +// } +// } +// +//} +//void MixingScheme::copyDensityFromOutHist(std::shared_ptr>> &rhoOutValues) +//{ +// unsigned int numQuadPointsPerCell = d_numberQuadraturePointsPerCell; +// if (d_dftParamsPtr->spinPolarized == 1) +// { +// numQuadPointsPerCell = 2*numQuadPointsPerCell; +// } +// +// unsigned int hist =d_rhoOutVals.size(); +// +// unsigned int iElem = 0; +// +// d_dofHandler = +// &d_matrixFreeDataPtr->get_dof_handler(d_matrixFreeVectorComponent); +// +// typename dealii::DoFHandler<3>::active_cell_iterator +// cell = d_dofHandler->begin_active(), +// endc = d_dofHandler->end(); +// for (; cell != endc; ++cell) +// { +// if (cell->is_locally_owned()) +// { +// for(unsigned int iQuad = 0 ;iQuad < numQuadPointsPerCell; iQuad++) +// { +// (*rhoOutValues)[cell->id()][iQuad] = +// d_rhoOutVals[hist-1][iElem*numQuadPointsPerCell + +// iQuad]; +// } +// iElem++; +// } +// } +// +//} +// +//void MixingScheme::copyDensityToInHist(std::shared_ptr>> &rhoInValues) +//{ +// std::vector latestHistRhoIn; +// unsigned int numQuadPointsPerCell = d_numberQuadraturePointsPerCell; +// if (d_dftParamsPtr->spinPolarized == 1) +// { +// numQuadPointsPerCell = 2*numQuadPointsPerCell; +// } +// +// latestHistRhoIn.resize(d_numCells*numQuadPointsPerCell); +// std::fill(latestHistRhoIn.begin(),latestHistRhoIn.end(),0.0); +// +// unsigned int iElem = 0; +// +// typename dealii::DoFHandler<3>::active_cell_iterator +// cell = d_dofHandler->begin_active(), +// endc = d_dofHandler->end(); +// for (; cell != endc; ++cell) +// { +// if (cell->is_locally_owned()) +// { +// for(unsigned int iQuad = 0 ;iQuad < numQuadPointsPerCell; iQuad++) +// { +// latestHistRhoIn[iElem*numQuadPointsPerCell + +// iQuad] = +// (*rhoInValues)[cell->id()][iQuad]; +// } +// iElem++; +// } +// } +// +// d_rhoInVals.push_back(latestHistRhoIn); +//}