Implemented and tested mean value constraints for full periodic elect…

…rostatic total potential solve.

include/poissonSolverProblem.h

@@ -49,7 +49,8 @@ namespace dftfe {
  const unsigned int matrixFreeVectorComponent,
  const std::map<dealii::types::global_dof_index, double> & atoms,
  const std::map<dealii::CellId,std::vector<double> > & rhoValues,
- const bool isComputeDiagonalA=true);
+ const bool isComputeDiagonalA=true,
+ const bool isComputeMeanValueConstraints=false);
 
  /**
  * @brief reinitialize data structures for nuclear electrostatic potential solve
@@ -125,6 +126,36 @@ namespace dftfe {
  */
  void computeDiagonalA();
 
+ /**
+ * @brief Compute mean value constraint which is required in case of fully periodic
+ * boundary conditions.
+ *
+ */
+ void computeMeanValueConstraint();
+
+
+ /**
+ * @brief Compute mean value constraint which is required in case of fully periodic
+ * boundary conditions.
+ *
+ */
+ void meanValueConstraintDistribute(vectorType& vec) const;
+
+ /**
+ * @brief Compute mean value constraint which is required in case of fully periodic
+ * boundary conditions.
+ *
+ */
+ void meanValueConstraintDistributeSlaveToMaster(vectorType& vec) const;
+
+
+ /**
+ * @brief Compute mean value constraint which is required in case of fully periodic
+ * boundary conditions.
+ *
+ */
+ void meanValueConstraintSetZero(vectorType& vec) const;
+
  /**
  * @brief precompute shape function gradient integral.
  *
@@ -159,6 +190,15 @@ namespace dftfe {
 
  bool d_isShapeGradIntegralPrecomputed;
 
+ /// storage for mean value constraint vector
+ vectorType d_meanValueConstraintVec;
+
+ /// boolean flag
+ bool d_isMeanValueConstraintComputed;
+
+ /// mean constrained nodeid
+ dealii::types::global_dof_index d_meanValueConstraintNodeId;
+
  const MPI_Comm mpi_communicator;
  const unsigned int n_mpi_processes;
  const unsigned int this_mpi_process;

src/dft/dft.cc

@@ -1117,7 +1117,9 @@ namespace dftfe {
  *d_constraintsVector[phiTotDofHandlerIndex],
  phiTotDofHandlerIndex,
  d_atomNodeIdToChargeMap,
- *rhoInValues);
+ *rhoInValues,
+ true,
+ dftParameters::periodicX && dftParameters::periodicY && dftParameters::periodicZ && !dftParameters::pinnedNodeForPBC);
 
  dealiiCGSolver.solve(phiTotalSolverProblem,
  dftParameters::absLinearSolverTolerance,
@@ -1129,20 +1131,8 @@ namespace dftfe {
  //
  if(dftParameters::periodicX && dftParameters::periodicY && dftParameters::periodicZ && !dftParameters::pinnedNodeForPBC)
  {
- double integPhi = totalCharge(dofHandler,
- d_phiTotRhoIn);
-
- double volume = computeVolume(dofHandler);
- double shiftingConst = integPhi/volume;
-
- vectorType tempDealiiVec;
- matrix_free_data.initialize_dof_vector(tempDealiiVec);
- tempDealiiVec = shiftingConst;
-
- d_phiTotRhoIn -= tempDealiiVec;
-
  if (dftParameters::verbosity>=2)
- pcout<<"Value of integPhiIn after scaling: "<<totalCharge(dofHandler,d_phiTotRhoIn)<<std::endl;
+ pcout<<"Value of integPhiIn: "<<totalCharge(dofHandler,d_phiTotRhoIn)<<std::endl;
  }
 
  computing_timer.exit_section("phiTot solve");
@@ -1526,20 +1516,8 @@ namespace dftfe {
  //
  if(dftParameters::periodicX && dftParameters::periodicY && dftParameters::periodicZ && !dftParameters::pinnedNodeForPBC)
  {
- double integPhi = totalCharge(dofHandler,
- d_phiTotRhoOut);
-
- double volume = computeVolume(dofHandler);
- double shiftingConst = integPhi/volume;
-
- vectorType tempDealiiVec;
- matrix_free_data.initialize_dof_vector(tempDealiiVec);
- tempDealiiVec = shiftingConst;
-
- d_phiTotRhoOut -= tempDealiiVec;
-
  if(dftParameters::verbosity>=2)
- pcout<<"Value of integPhiOut after scaling: "<<totalCharge(dofHandler,d_phiTotRhoOut);
+ pcout<<"Value of integPhiOut: "<<totalCharge(dofHandler,d_phiTotRhoOut);
  }
 
  computing_timer.exit_section("phiTot solve");

src/dft/electrostaticHRefinedEnergy.cc

@@ -438,9 +438,10 @@ void dftClass<FEOrder>::computeElectrostaticEnergyHRefined(const bool computeFor
  dealii::ConstraintMatrix constraintsForTotalPotential;
  constraintsForTotalPotential.reinit(locallyRelevantDofs);
 
- locatePeriodicPinnedNodes(dofHandlerHRefined,
- constraintsHRefined,
- constraintsForTotalPotential);
+ if (dftParameters::pinnedNodeForPBC)
+ locatePeriodicPinnedNodes(dofHandlerHRefined,
+ constraintsHRefined,
+ constraintsForTotalPotential);
  applyHomogeneousDirichletBC(dofHandlerHRefined,constraintsForTotalPotential);
  constraintsForTotalPotential.close();
 
@@ -554,7 +555,9 @@ void dftClass<FEOrder>::computeElectrostaticEnergyHRefined(const bool computeFor
  *matrixFreeConstraintsInputVector[phiTotDofHandlerIndexHRefined],
  phiTotDofHandlerIndexHRefined,
  atomHRefinedNodeIdToChargeMap,
- rhoOutHRefinedQuadValues);
+ rhoOutHRefinedQuadValues,
+ true,
+ dftParameters::periodicX && dftParameters::periodicY && dftParameters::periodicZ && !dftParameters::pinnedNodeForPBC);
 
  if (dftParameters::verbosity==2)
  pcout<< std::endl<<"Solving for total electrostatic potential (rhoIn+b) on h refined mesh: ";

src/dft/electrostaticPRefinedEnergy.cc

@@ -180,9 +180,11 @@ void dftClass<FEOrder>::computeElectrostaticEnergyPRefined()
  dealii::ConstraintMatrix constraintsForTotalPotential;
  constraintsForTotalPotential.reinit(locallyRelevantDofs);
 
- locatePeriodicPinnedNodes(dofHandlerPRefined,
- constraintsPRefined,
- constraintsForTotalPotential);
+ if (dftParameters::pinnedNodeForPBC)
+ locatePeriodicPinnedNodes(dofHandlerPRefined,
+ constraintsPRefined,
+ constraintsForTotalPotential);
+
  applyHomogeneousDirichletBC(dofHandlerPRefined,constraintsForTotalPotential);
  constraintsForTotalPotential.close();
 
@@ -347,7 +349,9 @@ void dftClass<FEOrder>::computeElectrostaticEnergyPRefined()
  *matrixFreeConstraintsInputVector[phiTotDofHandlerIndexPRefined],
  phiTotDofHandlerIndexPRefined,
  atomPRefinedNodeIdToChargeMap,
- rhoOutPRefinedQuadValues);
+ rhoOutPRefinedQuadValues,
+ true,
+ dftParameters::periodicX && dftParameters::periodicY && dftParameters::periodicZ && !dftParameters::pinnedNodeForPBC);
 
  if (dftParameters::verbosity==2)
  pcout<< std::endl<<"Solving for total electrostatic potential (rhoIn+b) on p refined mesh: ";

src/dft/initBoundaryConditions.cc

@@ -104,7 +104,8 @@ void dftClass<FEOrder>::initBoundaryConditions(){
  d_constraintsForTotalPotential.clear();
  d_constraintsForTotalPotential.reinit(locally_relevant_dofs);
 
- locatePeriodicPinnedNodes(dofHandler,constraintsNone,d_constraintsForTotalPotential);
+ if (dftParameters::pinnedNodeForPBC)
+ locatePeriodicPinnedNodes(dofHandler,constraintsNone,d_constraintsForTotalPotential);
  applyHomogeneousDirichletBC(dofHandler,d_constraintsForTotalPotential);
  d_constraintsForTotalPotential.close ();
 

src/poisson/poissonSolverProblem.cc

@@ -31,6 +31,7 @@ namespace dftfe {
  pcout (std::cout, (dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0))
  {
  d_isShapeGradIntegralPrecomputed=false;
+ d_isMeanValueConstraintComputed=false;
  }
 
  template<unsigned int FEOrder>
@@ -41,7 +42,8 @@ namespace dftfe {
  const unsigned int matrixFreeVectorComponent,
  const std::map<dealii::types::global_dof_index, double> & atoms,
  const std::map<dealii::CellId,std::vector<double> > & rhoValues,
- const bool isComputeDiagonalA)
+ const bool isComputeDiagonalA,
+ const bool isComputeMeanValueConstraint)
  {
  d_matrixFreeDataPtr=&matrixFreeData;
  d_xPtr=&x;
@@ -50,6 +52,12 @@ namespace dftfe {
  d_rhoValuesPtr=&rhoValues;
  d_atomsPtr=&atoms;
 
+ if (isComputeMeanValueConstraint)
+ {
+ computeMeanValueConstraint();
+ d_isMeanValueConstraintComputed=true;
+ }
+
  if (isComputeDiagonalA)
  computeDiagonalA();
  }
@@ -85,6 +93,9 @@ namespace dftfe {
  void poissonSolverProblem<FEOrder>::distributeX()
  {
  d_constraintMatrixPtr->distribute(*d_xPtr);
+
+ if (d_isMeanValueConstraintComputed)
+ meanValueConstraintDistribute(*d_xPtr);
  }
 
  template<unsigned int FEOrder>
@@ -139,6 +150,7 @@ namespace dftfe {
  {
 
  rhs.reinit(*d_xPtr);
+ rhs=0;
 
  const dealii::DoFHandler<3> & dofHandler=
  d_matrixFreeDataPtr->get_dof_handler(d_matrixFreeVectorComponent);
@@ -228,6 +240,9 @@ namespace dftfe {
  //MPI operation to sync data
  rhs.compress(dealii::VectorOperation::add);
 
+ if (d_isMeanValueConstraintComputed)
+ meanValueConstraintDistributeSlaveToMaster(rhs);
+
  //FIXME: check if this is really required
  d_constraintMatrixPtr->set_zero(rhs);
 
@@ -243,10 +258,122 @@ namespace dftfe {
  dst.scale(d_diagonalA);
  }
 
+ template<unsigned int FEOrder>
+ void poissonSolverProblem<FEOrder>::meanValueConstraintDistribute(vectorType& vec) const
+ {
+ const double constrainedNodeValue=d_meanValueConstraintVec*vec;
+ if (dealii::Utilities::MPI::this_mpi_process(mpi_communicator) ==0)
+ vec[d_meanValueConstraintNodeId]=constrainedNodeValue;
+ }
+
+ template<unsigned int FEOrder>
+ void poissonSolverProblem<FEOrder>::meanValueConstraintDistributeSlaveToMaster(vectorType& vec) const
+ {
+ double constrainedNodeValue=0;
+ if (dealii::Utilities::MPI::this_mpi_process(mpi_communicator) ==0)
+ constrainedNodeValue=vec[d_meanValueConstraintNodeId];
+
+ MPI_Bcast(&constrainedNodeValue,
+ 1,
+ MPI_DOUBLE,
+ 0,
+ mpi_communicator);
+
+ vec.add(constrainedNodeValue,d_meanValueConstraintVec);
+
+ meanValueConstraintSetZero(vec);
+ }
+
+ template<unsigned int FEOrder>
+ void poissonSolverProblem<FEOrder>::meanValueConstraintSetZero(vectorType& vec) const
+ {
+ if (d_isMeanValueConstraintComputed)
+ if (dealii::Utilities::MPI::this_mpi_process(mpi_communicator) ==0)
+ vec[d_meanValueConstraintNodeId]=0;
+ }
+
+ template<unsigned int FEOrder>
+ void poissonSolverProblem<FEOrder>::computeMeanValueConstraint()
+ {
+ d_meanValueConstraintVec.reinit(*d_xPtr);
+ d_meanValueConstraintVec=0;
+
+ const dealii::DoFHandler<3> & dofHandler=
+ d_matrixFreeDataPtr->get_dof_handler(d_matrixFreeVectorComponent);
+
+ dealii::QGauss<3> quadrature(C_num1DQuad<FEOrder>());
+ dealii::FEValues<3> fe_values (dofHandler.get_fe(), quadrature, dealii::update_values| dealii::update_JxW_values);
+ const unsigned int dofs_per_cell = dofHandler.get_fe().dofs_per_cell;
+ const unsigned int num_quad_points = quadrature.size();
+ dealii::Vector<double> elementalValues(dofs_per_cell);
+ std::vector<dealii::types::global_dof_index> local_dof_indices (dofs_per_cell);
+
+ //parallel loop over all elements
+ typename dealii::DoFHandler<3>::active_cell_iterator cell = dofHandler.begin_active(), endc = dofHandler.end();
+ for(; cell!=endc; ++cell)
+ if (cell->is_locally_owned())
+ {
+ fe_values.reinit (cell);
+
+ cell->get_dof_indices (local_dof_indices);
+
+ elementalValues=0.0;
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ for (unsigned int q_point = 0; q_point < num_quad_points; ++q_point)
+ elementalValues(i) += fe_values.shape_value (i, q_point)*fe_values.JxW(q_point);
+
+ d_constraintMatrixPtr->distribute_local_to_global(elementalValues,
+ local_dof_indices,
+ d_meanValueConstraintVec);
+ }
+
+ //MPI operation to sync data
+ d_meanValueConstraintVec.compress(dealii::VectorOperation::add);
+
+ dealii::IndexSet locallyOwnedElements=d_meanValueConstraintVec.locally_owned_elements();
+
+ dealii::IndexSet locallyRelevantElements;
+ dealii::DoFTools::extract_locally_relevant_dofs(dofHandler, locallyRelevantElements);
+
+ dealii::IndexSet allIndicesTouchedByConstraints(d_meanValueConstraintVec.size());
+ for (dealii::IndexSet::ElementIterator it=locallyRelevantElements.begin();
+ it<locallyRelevantElements.end();it++)
+ if (d_constraintMatrixPtr->is_constrained(*it))
+ {
+ const dealii::types::global_dof_index lineDof = *it;
+ const std::vector<std::pair<dealii::types::global_dof_index, double > > * rowData
+ =d_constraintMatrixPtr->get_constraint_entries(lineDof);
+ allIndicesTouchedByConstraints.add_index(lineDof);
+ for(unsigned int j = 0; j < rowData->size();++j)
+ allIndicesTouchedByConstraints.add_index((*rowData)[j].first);
+ }
+ for (std::map<dealii::types::global_dof_index, double>::const_iterator it=(*d_atomsPtr).begin(); it!=(*d_atomsPtr).end(); ++it)
+ allIndicesTouchedByConstraints.add_index(it->first);
+
+ locallyOwnedElements.subtract_set(allIndicesTouchedByConstraints);
+ d_meanValueConstraintNodeId=*locallyOwnedElements.begin();
+ AssertThrow(!d_constraintMatrixPtr->is_constrained(d_meanValueConstraintNodeId),dealii::ExcMessage("DFT-FE Error: BUG."));
+
+ //pcout<<"Mean value constraint nodeid: "<<d_meanValueConstraintNodeId<<std::endl;
+
+ double valueAtConstraintNode=d_meanValueConstraintVec[d_meanValueConstraintNodeId];
+ MPI_Bcast(&valueAtConstraintNode,
+ 1,
+ MPI_DOUBLE,
+ 0,
+ mpi_communicator);
+
+ d_meanValueConstraintVec/=-valueAtConstraintNode;
+ if (dealii::Utilities::MPI::this_mpi_process(mpi_communicator) ==0)
+ d_meanValueConstraintVec[d_meanValueConstraintNodeId]=0;
+ }
+
+
  template<unsigned int FEOrder>
  void poissonSolverProblem<FEOrder>::computeDiagonalA()
  {
  d_diagonalA.reinit(*d_xPtr);
+ d_diagonalA=0;
 
  const dealii::DoFHandler<3> & dofHandler=
  d_matrixFreeDataPtr->get_dof_handler(d_matrixFreeVectorComponent);
@@ -320,7 +447,18 @@ namespace dftfe {
  void poissonSolverProblem<FEOrder>::vmult(vectorType &Ax,const vectorType &x) const
  {
  Ax=0.0;
- d_matrixFreeDataPtr->cell_loop (&poissonSolverProblem<FEOrder>::AX, this, Ax, x);
+
+ if (d_isMeanValueConstraintComputed)
+ {
+ vectorType tempVec=x;
+ meanValueConstraintDistribute(tempVec);
+
+ d_matrixFreeDataPtr->cell_loop (&poissonSolverProblem<FEOrder>::AX, this, Ax, tempVec);
+
+ meanValueConstraintDistributeSlaveToMaster(Ax);
+ }
+ else
+ d_matrixFreeDataPtr->cell_loop (&poissonSolverProblem<FEOrder>::AX, this, Ax, x);
  }