From e41715ceaa8c1336c1092c52eeb1fc58ce3f8503 Mon Sep 17 00:00:00 2001 From: Sambit Das Date: Mon, 7 Oct 2019 22:20:51 -0400 Subject: [PATCH] Modified mesh generator to create consistent refinement across periodic surfaces. --- include/dft.h | 3 +- include/triangulationManager.h | 45 +- src/dft/dft.cc | 27 +- src/dft/moveAtoms.cc | 26 +- ...evOrthogonalizedSubspaceIterationSolver.cc | 16 +- .../triangulationManager/generateMesh.cc | 664 +++++++++++------- .../triangulationManager.cc | 25 +- utils/dftParameters.cc | 4 +- 8 files changed, 504 insertions(+), 306 deletions(-) diff --git a/include/dft.h b/include/dft.h index 8835187ec..809092088 100644 --- a/include/dft.h +++ b/include/dft.h @@ -137,7 +137,8 @@ namespace dftfe { /** * @brief Does KSDFT problem pre-processing steps but without remeshing. */ - void initNoRemesh(bool flag = true); + void initNoRemesh(const bool updateImageKPoints = true, + const bool useSingleAtomSolution = false ); /** * @brief Selects between only electronic field relaxation or combined electronic and geometry relxation diff --git a/include/triangulationManager.h b/include/triangulationManager.h index 7afa07b51..ddc30ede0 100644 --- a/include/triangulationManager.h +++ b/include/triangulationManager.h @@ -197,13 +197,13 @@ namespace dftfe { parallel::distributed::Triangulation<3> & parallelTriangulationB); /** - * @brief generates reset meshes to the last mesh generated by auto mesh approach. This + * @brief generates reset meshes to the last mesh generated by auto mesh approach. This * is used in Gaussian update of meshes during electrostatics */ void generateResetMeshes(const std::vector > & domainBoundingVectors, const bool generateSerialTria, const bool generateElectrostaticsTria); - + /** * @brief serialize the triangulations and the associated solution vectors @@ -277,17 +277,8 @@ namespace dftfe { parallel::distributed::Triangulation<3>& electrostaticsTriangulationRho, parallel::distributed::Triangulation<3>& electrostaticsTriangulationDisp, parallel::distributed::Triangulation<3>& electrostaticsTriangulationForce, - const bool generateElectrostaticsTria); - - /** - * @brief internal function which generates a parallel mesh using a adaptive refinement strategy. - * - */ - void generateMesh(parallel::distributed::Triangulation<3>& parallelTriangulation, - parallel::distributed::Triangulation<3>& electrostaticsTriangulationRho, - parallel::distributed::Triangulation<3>& electrostaticsTriangulationDisp, - parallel::distributed::Triangulation<3>& electrostaticsTriangulationForce, - const bool generateElectrostaticsTria); + const bool generateElectrostaticsTria, + const bool generateSerialTria=false); /** * @brief internal function which generates a coarse mesh which is required for the load function call in @@ -299,8 +290,24 @@ namespace dftfe { /** * @brief internal function which sets refinement flags based on a custom created algorithm * + * @return bool boolean flag is any local cell has refinement flag set + */ + bool refinementAlgorithmA(parallel::distributed::Triangulation<3>& parallelTriangulation, + parallel::distributed::Triangulation<3>& electrostaticsTriangulationRho, + parallel::distributed::Triangulation<3>& electrostaticsTriangulationDisp, + parallel::distributed::Triangulation<3>& electrostaticsTriangulationForce, + const bool generateElectrostaticsTria, + std::vector & locallyOwnedCellsRefineFlags, + std::map & cellIdToCellRefineFlagMapLocal, + const bool smoothenCellsOnPeriodicBoundary=false); + + /** + * @brief internal function which sets refinement flags to have consistent refinement across periodic + * boundary + * + * @return bool boolean flag is any local cell has refinement flag set */ - void refinementAlgorithmA(parallel::distributed::Triangulation<3>& parallelTriangulation, + bool consistentPeriodicBoundaryRefinement(parallel::distributed::Triangulation<3>& parallelTriangulation, parallel::distributed::Triangulation<3>& electrostaticsTriangulationRho, parallel::distributed::Triangulation<3>& electrostaticsTriangulationDisp, parallel::distributed::Triangulation<3>& electrostaticsTriangulationForce, @@ -308,7 +315,13 @@ namespace dftfe { std::vector & locallyOwnedCellsRefineFlags, std::map & cellIdToCellRefineFlagMapLocal); - + /** + * @brief check that triangulation has consistent refinement across periodic boundary including + * for ghost cells + * + */ + bool checkPeriodicSurfaceRefinementConsistency(parallel::distributed::Triangulation<3>& parallelTriangulation); + /** * @brief internal function which refines the serial mesh based on refinement flags from parallel mesh. @@ -351,7 +364,7 @@ namespace dftfe { std::vector > d_atomPositions; std::vector > d_imageAtomPositions; std::vector > d_domainBoundingVectors; - const unsigned int d_max_refinement_steps=20; + const unsigned int d_max_refinement_steps=40; // //parallel objects diff --git a/src/dft/dft.cc b/src/dft/dft.cc index 200e0ec7c..cda4c2692 100644 --- a/src/dft/dft.cc +++ b/src/dft/dft.cc @@ -338,7 +338,7 @@ namespace dftfe { "to half the number of electrons with a 15 percent buffer to avoid convergence issues in" "SCF iterations"<= 1) { @@ -645,10 +645,11 @@ namespace dftfe { } template - void dftClass::initNoRemesh(bool flag) + void dftClass::initNoRemesh(const bool updateImageKPoints, + const bool useSingleAtomSolution) { computingTimerStandard.enter_section("KSDFT problem initialization"); - if(flag) + if(updateImageKPoints) initImageChargesUpdateKPoints(); if (dftParameters::isIonOpt) @@ -658,11 +659,19 @@ namespace dftfe { // initBoundaryConditions(); - // - //rho init (use previous ground state electron density) - // - solveNoSCF(); - noRemeshRhoDataInit(); + if (useSingleAtomSolution) + { + readPSI(); + initRho(); + } + else + { + // + //rho init (use previous ground state electron density) + // + solveNoSCF(); + noRemeshRhoDataInit(); + } // //reinitialize pseudopotential related data structures @@ -1583,7 +1592,7 @@ namespace dftfe { true); d_groundStateEnergy = totalEnergy; - + } MPI_Barrier(interpoolcomm); diff --git a/src/dft/moveAtoms.cc b/src/dft/moveAtoms.cc index b998e7151..58d710364 100644 --- a/src/dft/moveAtoms.cc +++ b/src/dft/moveAtoms.cc @@ -100,7 +100,7 @@ namespace internal{ template void dftClass::updateAtomPositionsAndMoveMesh(const std::vector > & globalAtomsDisplacements, double maximumForceToBeRelaxed) - + { const int numberGlobalAtoms = atomLocations.size(); int numberImageCharges = d_imageIds.size(); @@ -164,7 +164,7 @@ void dftClass::updateAtomPositionsAndMoveMesh(const std::vector::updateAtomPositionsAndMoveMesh(const std::vector::updateAtomPositionsAndMoveMesh(const std::vector::updateAtomPositionsAndMoveMesh(const std::vector::updateAtomPositionsAndMoveMesh(const std::vector::updateAtomPositionsAndMoveMesh(const std::vectorgetGaussianGeneratorParameter()<<" as net max displacement magnitude: "<< maxDispAtom<< " is below " << break1 <<" Bohr"<::updateAtomPositionsAndMoveMesh(const std::vector factor*d_autoMeshMaxJacobianRatio) d_autoMesh=1; MPI_Bcast(&(d_autoMesh), @@ -437,12 +437,12 @@ void dftClass::updateAtomPositionsAndMoveMesh(const std::vector0.06?true:false); if (!dftParameters::reproducible_output) pcout << "...Reinitialization end" << std::endl; } diff --git a/src/solvers/eigenSolvers/chebyshevOrthogonalizedSubspaceIterationSolver.cc b/src/solvers/eigenSolvers/chebyshevOrthogonalizedSubspaceIterationSolver.cc index b1ea1378e..4310343b4 100755 --- a/src/solvers/eigenSolvers/chebyshevOrthogonalizedSubspaceIterationSolver.cc +++ b/src/solvers/eigenSolvers/chebyshevOrthogonalizedSubspaceIterationSolver.cc @@ -33,21 +33,21 @@ namespace dftfe{ else if(upperBoundUnwantedSpectrum > 500 && upperBoundUnwantedSpectrum <= 1000) chebyshevOrder = 30; else if(upperBoundUnwantedSpectrum > 1000 && upperBoundUnwantedSpectrum <= 1500) - chebyshevOrder = 34; + chebyshevOrder = 39; else if(upperBoundUnwantedSpectrum > 1500 && upperBoundUnwantedSpectrum <= 2000) - chebyshevOrder = 38; + chebyshevOrder = 43; else if(upperBoundUnwantedSpectrum > 2000 && upperBoundUnwantedSpectrum <= 3000) - chebyshevOrder = 45; + chebyshevOrder = 50; else if(upperBoundUnwantedSpectrum > 3000 && upperBoundUnwantedSpectrum <= 4000) - chebyshevOrder = 53; + chebyshevOrder = 58; else if(upperBoundUnwantedSpectrum > 4000 && upperBoundUnwantedSpectrum <= 5000) - chebyshevOrder = 60; + chebyshevOrder = 65; else if(upperBoundUnwantedSpectrum > 5000 && upperBoundUnwantedSpectrum <= 9000) - chebyshevOrder = 68; + chebyshevOrder = 73; else if(upperBoundUnwantedSpectrum > 9000 && upperBoundUnwantedSpectrum <= 14000) - chebyshevOrder = 94; + chebyshevOrder = 100; else if(upperBoundUnwantedSpectrum > 14000 && upperBoundUnwantedSpectrum <= 20000) - chebyshevOrder = 109; + chebyshevOrder = 115; else if(upperBoundUnwantedSpectrum > 20000 && upperBoundUnwantedSpectrum <= 30000) chebyshevOrder = 162; else if(upperBoundUnwantedSpectrum > 30000 && upperBoundUnwantedSpectrum <= 50000) diff --git a/src/triangulation/triangulationManager/generateMesh.cc b/src/triangulation/triangulationManager/generateMesh.cc index 9972a4a62..1d3f9851d 100644 --- a/src/triangulation/triangulationManager/generateMesh.cc +++ b/src/triangulation/triangulationManager/generateMesh.cc @@ -235,13 +235,14 @@ namespace dftfe { pcout< & parallelTriangulation, + bool triangulationManager::refinementAlgorithmA(parallel::distributed::Triangulation<3> & parallelTriangulation, parallel::distributed::Triangulation<3> & electrostaticsTriangulationRho, parallel::distributed::Triangulation<3> & electrostaticsTriangulationDisp, parallel::distributed::Triangulation<3> & electrostaticsTriangulationForce, const bool generateElectrostaticsTria, std::vector & locallyOwnedCellsRefineFlags, - std::map & cellIdToCellRefineFlagMapLocal) + std::map & cellIdToCellRefineFlagMapLocal, + const bool smoothenCellsOnPeriodicBoundary) { // //compute magnitudes of domainBounding Vectors @@ -263,6 +264,10 @@ namespace dftfe { cellElectroForce = electrostaticsTriangulationForce.begin_active(); } + std::map cellIdToLocallyOwnedId; + unsigned int locallyOwnedCount=0; + + bool isAnyCellRefined=false; // // // @@ -270,6 +275,9 @@ namespace dftfe { { if(cell->is_locally_owned()) { + cellIdToLocallyOwnedId[cell->id()]=locallyOwnedCount; + locallyOwnedCount++; + const dealii::Point<3> center(cell->center()); double currentMeshSize = cell->minimum_vertex_distance(); @@ -380,6 +388,7 @@ namespace dftfe { locallyOwnedCellsRefineFlags.push_back(1); cellIdToCellRefineFlagMapLocal[cell->id()]=1; cell->set_refine_flag(); + isAnyCellRefined=true; if(generateElectrostaticsTria) { cellElectroRho->set_refine_flag(); @@ -402,8 +411,165 @@ namespace dftfe { } + // + // refine cells on periodic boundary if their length is greater than two times the + // mesh size around atom + // + if (smoothenCellsOnPeriodicBoundary) + { + locallyOwnedCount=0; + cell = parallelTriangulation.begin_active(); + endc = parallelTriangulation.end(); + + if(generateElectrostaticsTria) + { + cellElectroRho = electrostaticsTriangulationRho.begin_active(); + cellElectroDisp = electrostaticsTriangulationDisp.begin_active(); + cellElectroForce = electrostaticsTriangulationForce.begin_active(); + } + + const unsigned int faces_per_cell=dealii::GeometryInfo<3>::faces_per_cell; + + for(;cell != endc; ++cell) + { + if (cell->is_locally_owned()) + { + if(cell->at_boundary() + && cell->minimum_vertex_distance()>2.0*dftParameters::meshSizeOuterBall + && !cell->refine_flag_set() ) + for(unsigned int iFace = 0; iFace < faces_per_cell; ++iFace) + if (cell->has_periodic_neighbor(iFace)) + { + cell->set_refine_flag(); + isAnyCellRefined=true; + locallyOwnedCellsRefineFlags[cellIdToLocallyOwnedId[cell->id()]]=1; + cellIdToCellRefineFlagMapLocal[cell->id()]=1; + if(generateElectrostaticsTria) + { + cellElectroRho->set_refine_flag(); + cellElectroDisp->set_refine_flag(); + cellElectroForce->set_refine_flag(); + } + break; + } + locallyOwnedCount++; + } + if(generateElectrostaticsTria) + { + ++cellElectroRho; + ++cellElectroDisp; + ++cellElectroForce; + } + } + } + + return isAnyCellRefined; + } + + // + // internal function which sets refinement flags to have consistent refinement across periodic boundary + // + bool triangulationManager::consistentPeriodicBoundaryRefinement(parallel::distributed::Triangulation<3> & parallelTriangulation, + parallel::distributed::Triangulation<3> & electrostaticsTriangulationRho, + parallel::distributed::Triangulation<3> & electrostaticsTriangulationDisp, + parallel::distributed::Triangulation<3> & electrostaticsTriangulationForce, + const bool generateElectrostaticsTria, + std::vector & locallyOwnedCellsRefineFlags, + std::map & cellIdToCellRefineFlagMapLocal) + { + locallyOwnedCellsRefineFlags.clear(); + cellIdToCellRefineFlagMapLocal.clear(); + typename parallel::distributed::Triangulation<3>::active_cell_iterator cell, endc, cellElectroRho, cellElectroDisp, cellElectroForce; + cell = parallelTriangulation.begin_active(); + endc = parallelTriangulation.end(); + + std::map cellIdToLocallyOwnedId; + unsigned int locallyOwnedCount=0; + for(;cell != endc; ++cell) + if(cell->is_locally_owned()) + { + cellIdToLocallyOwnedId[cell->id()]=locallyOwnedCount; + locallyOwnedCellsRefineFlags.push_back(0); + cellIdToCellRefineFlagMapLocal[cell->id()]=0; + locallyOwnedCount++; + } + + + cell = parallelTriangulation.begin_active(); + endc = parallelTriangulation.end(); + + if(generateElectrostaticsTria) + { + cellElectroRho = electrostaticsTriangulationRho.begin_active(); + cellElectroDisp = electrostaticsTriangulationDisp.begin_active(); + cellElectroForce = electrostaticsTriangulationForce.begin_active(); + } + + const unsigned int faces_per_cell=dealii::GeometryInfo<3>::faces_per_cell; + + bool isAnyCellRefined=false; + for(;cell != endc; ++cell) + { + if((cell->is_locally_owned() || cell->is_ghost()) && + cell->at_boundary()) + for(unsigned int iFace = 0; iFace < faces_per_cell; ++iFace) + if (cell->has_periodic_neighbor(iFace)) + if (cell->periodic_neighbor_is_coarser(iFace)) + { + typename parallel::distributed::Triangulation<3>::active_cell_iterator periodicCell=cell->periodic_neighbor(iFace); + + if (periodicCell->is_locally_owned()) + { + locallyOwnedCellsRefineFlags[cellIdToLocallyOwnedId[periodicCell->id()]]=1; + cellIdToCellRefineFlagMapLocal[periodicCell->id()]=1; + periodicCell->set_refine_flag(); + + isAnyCellRefined=true; + if(generateElectrostaticsTria) + { + cellElectroRho->periodic_neighbor(iFace)->set_refine_flag(); + cellElectroDisp->periodic_neighbor(iFace)->set_refine_flag(); + cellElectroForce->periodic_neighbor(iFace)->set_refine_flag(); + } + } + } + if(generateElectrostaticsTria) + { + ++cellElectroRho; + ++cellElectroDisp; + ++cellElectroForce; + } + } + return isAnyCellRefined; } + // + // check that triangulation has consistent refinement across periodic boundary including + // + bool triangulationManager::checkPeriodicSurfaceRefinementConsistency(parallel::distributed::Triangulation<3>& parallelTriangulation) + { + typename parallel::distributed::Triangulation<3>::active_cell_iterator cell, endc; + cell = parallelTriangulation.begin_active(); + endc = parallelTriangulation.end(); + + const unsigned int faces_per_cell=dealii::GeometryInfo<3>::faces_per_cell; + + unsigned int notConsistent=0; + for(;cell != endc; ++cell) + if((cell->is_locally_owned() || cell->is_ghost()) && + cell->at_boundary()) + for(unsigned int iFace = 0; iFace < faces_per_cell; ++iFace) + if (cell->has_periodic_neighbor(iFace)) + { + typename parallel::distributed::Triangulation<3>::active_cell_iterator periodicCell=cell->periodic_neighbor(iFace); + if (periodicCell->is_locally_owned() || cell->is_locally_owned()) + if (cell->periodic_neighbor_is_coarser(iFace) || periodicCell->has_children()) + notConsistent=1; + + } + notConsistent=Utilities::MPI::max(notConsistent, mpi_communicator); + return notConsistent==1?false:true; + } // //generate mesh based on a-posteriori estimates @@ -547,209 +713,6 @@ namespace dftfe { } - // - //generate adaptive mesh - // - - void triangulationManager::generateMesh(parallel::distributed::Triangulation<3> & parallelTriangulation, - parallel::distributed::Triangulation<3> & electrostaticsTriangulationRho, - parallel::distributed::Triangulation<3> & electrostaticsTriangulationDisp, - parallel::distributed::Triangulation<3> & electrostaticsTriangulationForce, - const bool generateElectrostaticsTria) - { - if(!dftParameters::meshFileName.empty()) - { - GridIn<3> gridinParallel; - gridinParallel.attach_triangulation(parallelTriangulation); - - // - //Read mesh in UCD format generated from Cubit - // - std::ifstream f1(dftParameters::meshFileName.c_str()); - gridinParallel.read_ucd(f1); - - meshGenUtils::markPeriodicFacesNonOrthogonal(parallelTriangulation,d_domainBoundingVectors); - } - else - { - - generateCoarseMesh(parallelTriangulation); - if(generateElectrostaticsTria) - { - generateCoarseMesh(electrostaticsTriangulationRho); - generateCoarseMesh(electrostaticsTriangulationDisp); - generateCoarseMesh(electrostaticsTriangulationForce); - AssertThrow(parallelTriangulation.n_global_active_cells()==electrostaticsTriangulationRho.n_global_active_cells(),ExcMessage("Number of coarse mesh cells are different in electrostatics triangulations having rho field.")); - AssertThrow(parallelTriangulation.n_global_active_cells()==electrostaticsTriangulationDisp.n_global_active_cells(),ExcMessage("Number of coarse mesh cells are different in electrostatics triangulations having disp field.")); - AssertThrow(parallelTriangulation.n_global_active_cells()==electrostaticsTriangulationForce.n_global_active_cells(),ExcMessage("Number of coarse mesh cells are different in electrostatics triangulation for force computation.")); - } - - d_parallelTriaCurrentRefinement.clear(); - // - //Multilayer refinement - // - unsigned int numLevels=0; - bool refineFlag = true; - - while(refineFlag) - { - refineFlag = false; - - std::vector locallyOwnedCellsRefineFlags; - std::map cellIdToCellRefineFlagMapLocal; - refinementAlgorithmA(parallelTriangulation, - electrostaticsTriangulationRho, - electrostaticsTriangulationDisp, - electrostaticsTriangulationForce, - generateElectrostaticsTria, - locallyOwnedCellsRefineFlags, - cellIdToCellRefineFlagMapLocal); - - //This sets the global refinement sweep flag - refineFlag = std::accumulate(locallyOwnedCellsRefineFlags.begin(), - locallyOwnedCellsRefineFlags.end(), 0)>0?1:0; - refineFlag= Utilities::MPI::max((unsigned int) refineFlag, mpi_communicator); - - if (refineFlag) - { - if(numLevels=4) - pcout<< "refinement in progress, level: "<< numLevels<()); - parallelTriangulation.save_refine_flags(d_parallelTriaCurrentRefinement[numLevels]); - - parallelTriangulation.execute_coarsening_and_refinement(); - if(generateElectrostaticsTria) - { - electrostaticsTriangulationRho.execute_coarsening_and_refinement(); - electrostaticsTriangulationDisp.execute_coarsening_and_refinement(); - electrostaticsTriangulationForce.execute_coarsening_and_refinement(); - } - numLevels++; - } - else - { - refineFlag=false; - } - } - } - // - //compute some adaptive mesh metrics - // - double minElemLength = dftParameters::meshSizeOuterDomain; - double maxElemLength = 0.0; - unsigned int numLocallyOwnedCells=0; - typename parallel::distributed::Triangulation<3>::active_cell_iterator cell, endc, cellDisp, cellForce; - cell = parallelTriangulation.begin_active(); - endc = parallelTriangulation.end(); - for( ; cell != endc; ++cell) - { - if(cell->is_locally_owned()) - { - numLocallyOwnedCells++; - if(cell->minimum_vertex_distance() < minElemLength) - minElemLength = cell->minimum_vertex_distance(); - - if(cell->minimum_vertex_distance() > maxElemLength) - maxElemLength = cell->minimum_vertex_distance(); - } - } - - minElemLength = Utilities::MPI::min(minElemLength, mpi_communicator); - maxElemLength = Utilities::MPI::max(maxElemLength, mpi_communicator); - - // - //print out adaptive mesh metrics - // - if (dftParameters::verbosity>=4) - { - pcout<< "Triangulation generation summary: "<is_locally_owned()) - { - numLocallyOwnedCells++; - if(cell->minimum_vertex_distance() < minElemLengthRho) - minElemLengthRho = cell->minimum_vertex_distance(); - if(cellDisp->minimum_vertex_distance() < minElemLengthDisp) - minElemLengthDisp = cellDisp->minimum_vertex_distance(); - if(cellForce->minimum_vertex_distance() < minElemLengthForce) - minElemLengthForce = cellForce->minimum_vertex_distance(); - } - ++cellDisp; - ++cellForce; - } - - minElemLengthRho = Utilities::MPI::min(minElemLengthRho, mpi_communicator); - minElemLengthDisp = Utilities::MPI::min(minElemLengthDisp, mpi_communicator); - minElemLengthForce = Utilities::MPI::min(minElemLengthForce, mpi_communicator); - // - //print out adaptive electrostatics mesh metrics - // - if (dftParameters::verbosity>=4) - { - pcout<< "Electrostatics Triangulation generation summary: "< & parallelTriangulation, @@ -758,13 +721,15 @@ namespace dftfe { parallel::distributed::Triangulation<3> & electrostaticsTriangulationRho, parallel::distributed::Triangulation<3> & electrostaticsTriangulationDisp, parallel::distributed::Triangulation<3> & electrostaticsTriangulationForce, - const bool generateElectrostaticsTria) + const bool generateElectrostaticsTria, + const bool generateSerialTria) { if(!dftParameters::meshFileName.empty()) { GridIn<3> gridinParallel, gridinSerial; gridinParallel.attach_triangulation(parallelTriangulation); - gridinSerial.attach_triangulation(serialTriangulation); + if (generateSerialTria) + gridinSerial.attach_triangulation(serialTriangulation); // //Read mesh in UCD format generated from Cubit @@ -772,17 +737,22 @@ namespace dftfe { std::ifstream f1(dftParameters::meshFileName.c_str()); std::ifstream f2(dftParameters::meshFileName.c_str()); gridinParallel.read_ucd(f1); - gridinSerial.read_ucd(f2); + if (generateSerialTria) + gridinSerial.read_ucd(f2); meshGenUtils::markPeriodicFacesNonOrthogonal(parallelTriangulation,d_domainBoundingVectors); - meshGenUtils::markPeriodicFacesNonOrthogonal(serialTriangulation,d_domainBoundingVectors); + if (generateSerialTria) + meshGenUtils::markPeriodicFacesNonOrthogonal(serialTriangulation,d_domainBoundingVectors); } else { generateCoarseMesh(parallelTriangulation); - generateCoarseMesh(serialTriangulation); - AssertThrow(parallelTriangulation.n_global_active_cells()==serialTriangulation.n_global_active_cells(),ExcMessage("Number of coarse mesh cells are different in serial and parallel triangulations.")); + if (generateSerialTria) + { + generateCoarseMesh(serialTriangulation); + AssertThrow(parallelTriangulation.n_global_active_cells()==serialTriangulation.n_global_active_cells(),ExcMessage("Number of coarse mesh cells are different in serial and parallel triangulations.")); + } if(generateElectrostaticsTria) { @@ -793,12 +763,16 @@ namespace dftfe { AssertThrow(parallelTriangulation.n_global_active_cells()==electrostaticsTriangulationDisp.n_global_active_cells(),ExcMessage("Number of coarse mesh cells are different in electrostatics triangulations disp field.")); AssertThrow(parallelTriangulation.n_global_active_cells()==electrostaticsTriangulationForce.n_global_active_cells(),ExcMessage("Number of coarse mesh cells are different in electrostatics triangulations for force computation.")); - generateCoarseMesh(serialTriangulationElectrostatics); - AssertThrow(parallelTriangulation.n_global_active_cells()==serialTriangulationElectrostatics.n_global_active_cells(),ExcMessage("Number of coarse mesh cells are different in electrostatics serial triangulation computation.")); + if (generateSerialTria) + { + generateCoarseMesh(serialTriangulationElectrostatics); + AssertThrow(parallelTriangulation.n_global_active_cells()==serialTriangulationElectrostatics.n_global_active_cells(),ExcMessage("Number of coarse mesh cells are different in electrostatics serial triangulation computation.")); + } } d_parallelTriaCurrentRefinement.clear(); - d_serialTriaCurrentRefinement.clear(); + if (generateSerialTria) + d_serialTriaCurrentRefinement.clear(); // //Multilayer refinement @@ -810,19 +784,76 @@ namespace dftfe { refineFlag = false; std::vector locallyOwnedCellsRefineFlags; std::map cellIdToCellRefineFlagMapLocal; - refinementAlgorithmA(parallelTriangulation, - electrostaticsTriangulationRho, - electrostaticsTriangulationDisp, - electrostaticsTriangulationForce, - generateElectrostaticsTria, - locallyOwnedCellsRefineFlags, - cellIdToCellRefineFlagMapLocal); - - - //This sets the global refinement sweep flag - refineFlag = std::accumulate(locallyOwnedCellsRefineFlags.begin(), - locallyOwnedCellsRefineFlags.end(), 0)>0?1:0; - refineFlag= Utilities::MPI::max((unsigned int) refineFlag, mpi_communicator); + if (!dftParameters::reproducible_output) + { + if (numLevels%2==0) + { + refineFlag=refinementAlgorithmA(parallelTriangulation, + electrostaticsTriangulationRho, + electrostaticsTriangulationDisp, + electrostaticsTriangulationForce, + generateElectrostaticsTria, + locallyOwnedCellsRefineFlags, + cellIdToCellRefineFlagMapLocal); + + //This sets the global refinement sweep flag + refineFlag= Utilities::MPI::max((unsigned int) refineFlag, mpi_communicator); + + if (!refineFlag) + { + refineFlag=consistentPeriodicBoundaryRefinement(parallelTriangulation, + electrostaticsTriangulationRho, + electrostaticsTriangulationDisp, + electrostaticsTriangulationForce, + generateElectrostaticsTria, + locallyOwnedCellsRefineFlags, + cellIdToCellRefineFlagMapLocal); + + //This sets the global refinement sweep flag + refineFlag= Utilities::MPI::max((unsigned int) refineFlag, mpi_communicator); + } + } + else + { + refineFlag=consistentPeriodicBoundaryRefinement(parallelTriangulation, + electrostaticsTriangulationRho, + electrostaticsTriangulationDisp, + electrostaticsTriangulationForce, + generateElectrostaticsTria, + locallyOwnedCellsRefineFlags, + cellIdToCellRefineFlagMapLocal); + + //This sets the global refinement sweep flag + refineFlag= Utilities::MPI::max((unsigned int) refineFlag, mpi_communicator); + + if (!refineFlag) + { + refineFlag=refinementAlgorithmA(parallelTriangulation, + electrostaticsTriangulationRho, + electrostaticsTriangulationDisp, + electrostaticsTriangulationForce, + generateElectrostaticsTria, + locallyOwnedCellsRefineFlags, + cellIdToCellRefineFlagMapLocal); + + //This sets the global refinement sweep flag + refineFlag= Utilities::MPI::max((unsigned int) refineFlag, mpi_communicator); + } + } + } + else + { + refineFlag=refinementAlgorithmA(parallelTriangulation, + electrostaticsTriangulationRho, + electrostaticsTriangulationDisp, + electrostaticsTriangulationForce, + generateElectrostaticsTria, + locallyOwnedCellsRefineFlags, + cellIdToCellRefineFlagMapLocal); + + //This sets the global refinement sweep flag + refineFlag= Utilities::MPI::max((unsigned int) refineFlag, mpi_communicator); + } //Refine if (refineFlag) @@ -833,21 +864,25 @@ namespace dftfe { if (dftParameters::verbosity>=4) pcout<< "refinement in progress, level: "<< numLevels<()); - - //First refine serial mesh - refineSerialMesh(cellIdToCellRefineFlagMapLocal, - mpi_communicator, - serialTriangulation, - parallelTriangulation, - d_serialTriaCurrentRefinement[numLevels]) ; + if (generateSerialTria) + { + d_serialTriaCurrentRefinement.push_back(std::vector()); - if(generateElectrostaticsTria) + //First refine serial mesh refineSerialMesh(cellIdToCellRefineFlagMapLocal, mpi_communicator, - serialTriangulationElectrostatics, + serialTriangulation, parallelTriangulation, - d_serialTriaCurrentRefinement[numLevels]); + d_serialTriaCurrentRefinement[numLevels]) ; + + + if(generateElectrostaticsTria) + refineSerialMesh(cellIdToCellRefineFlagMapLocal, + mpi_communicator, + serialTriangulationElectrostatics, + parallelTriangulation, + d_serialTriaCurrentRefinement[numLevels]); + } d_parallelTriaCurrentRefinement.push_back(std::vector()); parallelTriangulation.save_refine_flags(d_parallelTriaCurrentRefinement[numLevels]); @@ -869,6 +904,139 @@ namespace dftfe { } } + + if (!dftParameters::reproducible_output) + { + if (!checkPeriodicSurfaceRefinementConsistency(parallelTriangulation)) + { + refineFlag=true; + while (refineFlag) + { + refineFlag = false; + std::vector locallyOwnedCellsRefineFlags; + std::map cellIdToCellRefineFlagMapLocal; + if (numLevels%2==0) + { + refineFlag=refinementAlgorithmA(parallelTriangulation, + electrostaticsTriangulationRho, + electrostaticsTriangulationDisp, + electrostaticsTriangulationForce, + generateElectrostaticsTria, + locallyOwnedCellsRefineFlags, + cellIdToCellRefineFlagMapLocal, + true); + + //This sets the global refinement sweep flag + refineFlag= Utilities::MPI::max((unsigned int) refineFlag, mpi_communicator); + + if (!refineFlag) + { + refineFlag=consistentPeriodicBoundaryRefinement(parallelTriangulation, + electrostaticsTriangulationRho, + electrostaticsTriangulationDisp, + electrostaticsTriangulationForce, + generateElectrostaticsTria, + locallyOwnedCellsRefineFlags, + cellIdToCellRefineFlagMapLocal); + + //This sets the global refinement sweep flag + refineFlag= Utilities::MPI::max((unsigned int) refineFlag, mpi_communicator); + } + } + else + { + refineFlag=consistentPeriodicBoundaryRefinement(parallelTriangulation, + electrostaticsTriangulationRho, + electrostaticsTriangulationDisp, + electrostaticsTriangulationForce, + generateElectrostaticsTria, + locallyOwnedCellsRefineFlags, + cellIdToCellRefineFlagMapLocal); + + //This sets the global refinement sweep flag + refineFlag= Utilities::MPI::max((unsigned int) refineFlag, mpi_communicator); + + if (!refineFlag) + { + refineFlag=refinementAlgorithmA(parallelTriangulation, + electrostaticsTriangulationRho, + electrostaticsTriangulationDisp, + electrostaticsTriangulationForce, + generateElectrostaticsTria, + locallyOwnedCellsRefineFlags, + cellIdToCellRefineFlagMapLocal, + true); + + //This sets the global refinement sweep flag + refineFlag= Utilities::MPI::max((unsigned int) refineFlag, mpi_communicator); + } + } + + //Refine + if (refineFlag) + { + + if(numLevels=4) + pcout<< "refinement in progress, level: "<< numLevels<()); + + //First refine serial mesh + refineSerialMesh(cellIdToCellRefineFlagMapLocal, + mpi_communicator, + serialTriangulation, + parallelTriangulation, + d_serialTriaCurrentRefinement[numLevels]) ; + + + if(generateElectrostaticsTria) + refineSerialMesh(cellIdToCellRefineFlagMapLocal, + mpi_communicator, + serialTriangulationElectrostatics, + parallelTriangulation, + d_serialTriaCurrentRefinement[numLevels]); + } + + d_parallelTriaCurrentRefinement.push_back(std::vector()); + parallelTriangulation.save_refine_flags(d_parallelTriaCurrentRefinement[numLevels]); + + parallelTriangulation.execute_coarsening_and_refinement(); + if(generateElectrostaticsTria) + { + electrostaticsTriangulationRho.execute_coarsening_and_refinement(); + electrostaticsTriangulationDisp.execute_coarsening_and_refinement(); + electrostaticsTriangulationForce.execute_coarsening_and_refinement(); + } + + numLevels++; + } + else + { + refineFlag=false; + } + } + + } + } + + if (checkPeriodicSurfaceRefinementConsistency(parallelTriangulation)) + { + if (dftParameters::verbosity>=4) + pcout<< "Periodic surface refinement consistency achieved."<=4) + pcout<< "Periodic surface refinement consistency could not be achieved."<=4) - pcout<<" numParallelCells: "<< numberGlobalCellsParallel<<", numSerialCells: "<< numberGlobalCellsSerial<=4) + pcout<<" numParallelCells: "<< numberGlobalCellsParallel<<", numSerialCells: "<< numberGlobalCellsSerial<=4) + pcout<<" numParallelCells: "<< numberGlobalCellsParallel<