Unverified Commit 7e33bc45 authored by Luc Berger-Vergiat's avatar Luc Berger-Vergiat
Browse files

MueLu: adding logic for piece-wise linear operator in regionMG see #5069

The idea is to only look at what happens to point injected from the fine level to the coarse level as these are equivalent to root points. These nodes are easily detected as their corresponding rows in the prolongator contain only one entry.
Also adding a couple of command line options to drive the solver more finely and using command line argument in linear convergence criteria.
parent 28e4bb4a
......@@ -973,9 +973,8 @@ void MakeCoarseLevelMaps2(const int maxRegPerGID,
using MT = typename Teuchos::ScalarTraits<SC>::magnitudeType;
const GO GO_INV = Teuchos::OrdinalTraits<GO>::invalid();
RCP<const Teuchos::Comm<int> > comm = regProlong[1][0]->getRowMap()->getComm();
const GO GO_INV = Teuchos::OrdinalTraits<GO>::invalid();
const int numLevels = regProlong.size();
Teuchos::Array<LO> coarseCompositeToRegionLIDs;
......@@ -1002,7 +1001,7 @@ void MakeCoarseLevelMaps2(const int maxRegPerGID,
}
}
// We want to communicate the coarse GIDs associated with each fine points.
// We gather the coarse GIDs associated with each fine point in the local composite mesh part.
RCP<Xpetra::Vector<MT,LO,GO,NO> > coarseCompositeGIDs
= Xpetra::VectorFactory<MT,LO,GO,NO>::Build(regRowImporters[currentLevel - 1][0]->getSourceMap(), false);
Teuchos::ArrayRCP<MT> coarseCompositeGIDsData = coarseCompositeGIDs->getDataNonConst(0);
......@@ -1013,9 +1012,15 @@ void MakeCoarseLevelMaps2(const int maxRegPerGID,
regProlong[currentLevel][0]->getLocalRowView(compositeToRegionLIDs[compositeNodeIdx],
coarseRegionLID,
dummyData);
coarseCompositeGIDsData[compositeNodeIdx] = regProlong[currentLevel][0]->getColMap()->getGlobalElement(coarseRegionLID[0]);
if(coarseRegionLID.size() == 1) {
coarseCompositeGIDsData[compositeNodeIdx] = regProlong[currentLevel][0]->getColMap()->getGlobalElement(coarseRegionLID[0]);
} else {
coarseCompositeGIDsData[compositeNodeIdx] = -1;
}
}
// We communicate the above GIDs to their duplicate so that we can replace GIDs of the region
// column map and form the quasiregion column map.
std::vector<RCP<Xpetra::Vector<MT, LO, GO, NO> > > coarseQuasiregionGIDs(1), coarseRegionGIDs(1);
compositeToRegional(coarseCompositeGIDs,
coarseQuasiregionGIDs,
......@@ -1045,7 +1050,9 @@ void MakeCoarseLevelMaps2(const int maxRegPerGID,
for(size_t regionIdx = 0; regionIdx < numCoarseRegionNodes; ++regionIdx) {
const GO initialValue = coarseQuasiregRowMapData[regionIdx];
for(size_t duplicateIdx = 0; duplicateIdx < numFineDuplicateNodes; ++duplicateIdx) {
if((initialValue == fineRegionDuplicateCoarseLIDs[duplicateIdx]) && (fineRegionDuplicateCoarseGIDs[duplicateIdx] < coarseQuasiregRowMapData[regionIdx])){
if((initialValue == fineRegionDuplicateCoarseLIDs[duplicateIdx]) &&
(fineRegionDuplicateCoarseGIDs[duplicateIdx] < coarseQuasiregRowMapData[regionIdx]) &&
(-1 < fineRegionDuplicateCoarseGIDs[duplicateIdx])){
coarseQuasiregRowMapData[regionIdx] = fineRegionDuplicateCoarseGIDs[duplicateIdx];
}
}
......
......@@ -12,7 +12,7 @@ IF (${PACKAGE_NAME}_ENABLE_Tpetra AND ${PACKAGE_NAME}_ENABLE_Amesos2)
)
TRIBITS_COPY_FILES_TO_BINARY_DIR(StructuredRegion_cp
SOURCE_FILES structured_1dof.xml
SOURCE_FILES structured_1dof.xml structured_linear_1dof.xml
)
TRIBITS_ADD_TEST(
......@@ -31,5 +31,21 @@ IF (${PACKAGE_NAME}_ENABLE_Tpetra AND ${PACKAGE_NAME}_ENABLE_Amesos2)
NUM_MPI_PROCS 4
)
TRIBITS_ADD_TEST(
StructuredRegionDriver
NAME "Structured_Region_Linear_Star2D_Tpetra"
ARGS "--linAlgebra=Tpetra --xml=structured_linear_1dof.xml --matrixType=Star2D --nx=10 --ny=10"
COMM serial mpi
NUM_MPI_PROCS 1
)
TRIBITS_ADD_TEST(
StructuredRegionDriver
NAME "Structured_Region_Linear_Star2D_Tpetra"
ARGS "--linAlgebra=Tpetra --xml=structured_linear_1dof.xml --matrixType=Star2D --nx=10 --ny=10"
COMM serial mpi
NUM_MPI_PROCS 4
)
ENDIF()
......@@ -159,7 +159,9 @@ int main_(Teuchos::CommandLineProcessor &clp, Xpetra::UnderlyingLib& lib, int ar
std::string xmlFileName = ""; clp.setOption("xml", &xmlFileName, "read parameters from an xml file");
std::string yamlFileName = ""; clp.setOption("yaml", &yamlFileName, "read parameters from a yaml file");
int maxIts = 200; clp.setOption("its", &maxIts, "maximum number of solver iterations");
double tol = 1e-12; clp.setOption("tol", &tol, "solver convergence tolerance");
int smootherIts = 20; clp.setOption("smootherIts", &smootherIts, "number of smoother iterations");
double smootherDamp = 0.67; clp.setOption("smootherDamp", &smootherDamp, "damping parameter for the level smoother");
double tol = 1e-12; clp.setOption("tol", &tol, "solver convergence tolerance");
bool scaleResidualHist = true; clp.setOption("scale", "noscale", &scaleResidualHist, "scaled Krylov residual history");
bool solvePreconditioned = true; clp.setOption("solve-preconditioned","no-solve-preconditioned", &solvePreconditioned, "use MueLu preconditioner in solve");
#ifdef HAVE_MUELU_TPETRA
......@@ -1282,10 +1284,8 @@ int main_(Teuchos::CommandLineProcessor &clp, Xpetra::UnderlyingLib& lib, int ar
/////////////////////////////////////////////////////////////////////////
// define max iteration counts
const int maxVCycle = 200;
const int maxFineIter = 20;
const int maxCoarseIter = 100;
const double omega = 0.67;
// Prepare output of residual norm to file
RCP<std::ofstream> log;
......@@ -1299,7 +1299,7 @@ int main_(Teuchos::CommandLineProcessor &clp, Xpetra::UnderlyingLib& lib, int ar
}
// Richardson iterations
for (int cycle = 0; cycle < maxVCycle; ++cycle) {
for (int cycle = 0; cycle < maxIts; ++cycle) {
// check for convergence
{
////////////////////////////////////////////////////////////////////////
......@@ -1323,7 +1323,7 @@ int main_(Teuchos::CommandLineProcessor &clp, Xpetra::UnderlyingLib& lib, int ar
(*log) << cycle << "\t" << normRes << "\n";
}
if (normRes < 1.0e-12)
if (normRes < tol)
break;
}
......@@ -1332,8 +1332,8 @@ int main_(Teuchos::CommandLineProcessor &clp, Xpetra::UnderlyingLib& lib, int ar
/////////////////////////////////////////////////////////////////////////
// printRegionalObject<Vector>("regB 2", regB, myRank, *fos);
vCycle(0, numLevels, maxFineIter, maxCoarseIter, omega, maxRegPerProc, regX, regB,
regMatrices,
vCycle(0, numLevels, smootherIts, maxCoarseIter, smootherDamp, maxRegPerProc,
regX, regB, regMatrices,
regProlong, compRowMaps, quasiRegRowMaps, regRowMaps, regRowImporters,
regInterfaceScalings, coarseCompOp, coarseCompositeDirectSolver);
}
......
......@@ -94,7 +94,7 @@
<Parameter name="max levels" type="int" value="6"/> <!-- Max number of levels -->
<Parameter name="cycle type" type="string" value="W"/>
<Parameter name="coarse: max size" type="int" value="10"/> <!-- Min number of rows on coarsest level -->
<Parameter name="coarse: max size" type="int" value="20"/> <!-- Min number of rows on coarsest level -->
<Parameter name="verbosity" type="string" value="High"/>
<ParameterList name="All">
......
<ParameterList name="MueLu">
<!-- Configuration of the Xpetra operator (fine level) -->
<ParameterList name="Matrix">
<Parameter name="PDE equations" type="int" value="1"/> <!-- Number of PDE equations at each grid node.-->
</ParameterList>
<!-- Factory collection -->
<ParameterList name="Factories">
<ParameterList name="myCoalesceDropFact">
<Parameter name="factory" type="string" value="CoalesceDropFactory"/>
<Parameter name="lightweight wrap" type="bool" value="true"/>
<Parameter name="aggregation: drop tol" type="double" value="0.00"/>
</ParameterList>
<ParameterList name="myAggregationFact">
<Parameter name="factory" type="string" value="StructuredAggregationFactory"/>
<Parameter name="aggregation: coupling" type="string" value="uncoupled"/>
<Parameter name="aggregation: output type" type="string" value="CrsGraph"/>
<Parameter name="aggregation: coarsening order" type="int" value="1"/>
<Parameter name="aggregation: coarsening rate" type="string" value="{2}"/>
<Parameter name="Graph" type="string" value="myCoalesceDropFact"/>
</ParameterList>
<ParameterList name="myCoarseMapFact">
<Parameter name="factory" type="string" value="CoarseMapFactory"/>
<Parameter name="Aggregates" type="string" value="myAggregationFact"/>
</ParameterList>
<!-- Note that ParameterLists must be defined prior to being used -->
<ParameterList name="myProlongatorFact">
<Parameter name="factory" type="string" value="GeometricInterpolationPFactory"/>
<Parameter name="interp: build coarse coordinates" type="bool" value="true"/>
<Parameter name="interp: interpolation order" type="int" value="1"/>
<Parameter name="prolongatorGraph" type="string" value="myAggregationFact"/>
<Parameter name="coarseCoordinatesFineMap" type="string" value="myAggregationFact"/>
<Parameter name="coarseCoordinatesMap" type="string" value="myAggregationFact"/>
</ParameterList>
<ParameterList name="myCoordTransferFact">
<Parameter name="factory" type="string" value="CoordinatesTransferFactory"/>
<Parameter name="structured aggregation" type="bool" value="true"/>
<Parameter name="numDimensions" type="string" value="myAggregationFact"/>
<Parameter name="lCoarseNodesPerDim" type="string" value="myAggregationFact"/>
</ParameterList>
<ParameterList name="myNullspaceFact">
<Parameter name="factory" type="string" value="NullspaceFactory"/>
<Parameter name="Nullspace" type="string" value="myProlongatorFact"/>
</ParameterList>
<ParameterList name="myRestrictorFact">
<Parameter name="factory" type="string" value="TransPFactory"/>
</ParameterList>
<!-- <ParameterList name="myAggExport"> -->
<!-- <Parameter name="factory" type="string" value="AggregationExportFactory"/> -->
<!-- <Parameter name="Aggregates" type="string" value="myAggregationFact"/> -->
<!-- <Parameter name="aggregation: output filename" type="string" value="structured_aggs"/> -->
<!-- <Parameter name="aggregation: output file: agg style" type="string" value="Jacks"/> -->
<!-- <Parameter name="aggregation: output file: agg style" type="string" value="Convex Hulls"/> -->
<!-- </ParameterList> -->
<ParameterList name="myRAPFact">
<Parameter name="factory" type="string" value="RAPFactory"/>
<Parameter name="P" type="string" value="myProlongatorFact"/>
<Parameter name="R" type="string" value="myRestrictorFact"/>
<ParameterList name="TransferFactories">
<Parameter name="CoordinateTransfer" type="string" value="myCoordTransferFact"/>
<!-- <Parameter name="AggregationExportFactory" type="string" value="myAggExport"/> -->
</ParameterList>
</ParameterList>
<ParameterList name="myILU">
<Parameter name="factory" type="string" value="TrilinosSmoother"/>
<Parameter name="type" type="string" value="RILUK"/>
<ParameterList name="ParameterList">
<Parameter name="schwarz: overlap level" type="int" value="1"/>
<Parameter name="schwarz: combine mode" type="string" value="Zero"/>
<Parameter name="schwarz: use reordering" type="bool" value="false"/>
<Parameter name="fact: iluk level-of-fill" type="int" value="0"/>
<Parameter name="fact: absolute threshold" type="double" value="0."/>
<Parameter name="fact: relative threshold" type="double" value="1."/>
<Parameter name="fact: relax value" type="double" value="0."/>
</ParameterList>
</ParameterList>
</ParameterList>
<!-- Definition of the multigrid preconditioner -->
<ParameterList name="Hierarchy">
<Parameter name="max levels" type="int" value="6"/> <!-- Max number of levels -->
<Parameter name="cycle type" type="string" value="W"/>
<Parameter name="coarse: max size" type="int" value="20"/> <!-- Min number of rows on coarsest level -->
<Parameter name="verbosity" type="string" value="High"/>
<ParameterList name="All">
<Parameter name="PreSmoother" type="string" value="myILU"/>
<Parameter name="PostSmoother" type="string" value="NoSmoother"/>
<Parameter name="Nullspace" type="string" value="myNullspaceFact"/>
<Parameter name="Aggregates" type="string" value="myAggregationFact"/>
<Parameter name="lCoarseNodesPerDim" type="string" value="myAggregationFact"/>
<Parameter name="P" type="string" value="myProlongatorFact"/>
<Parameter name="R" type="string" value="myRestrictorFact"/>
<Parameter name="A" type="string" value="myRAPFact"/>
<!-- <Parameter name="CoarseSolver" type="string" value="DirectSolver"/> -->
<Parameter name="CoarseSolver" type="string" value="myILU"/>
<Parameter name="Coordinates" type="string" value="myProlongatorFact"/>
<Parameter name="lNodesPerDim" type="string" value="myCoordTransferFact"/>
<Parameter name="numDimensions" type="string" value="myCoordTransferFact"/>
</ParameterList>
</ParameterList>
</ParameterList>
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment