diff --git a/src/cantorfaultnetworks/cantorfaultnetwork.cc b/src/cantorfaultnetworks/cantorfaultnetwork.cc index 82e7e419515a28360dc7e2a28d547369553f046a..870b2a2537f23f3ad5e37462a8afa3433dff4314 100644 --- a/src/cantorfaultnetworks/cantorfaultnetwork.cc +++ b/src/cantorfaultnetworks/cantorfaultnetwork.cc @@ -43,6 +43,7 @@ #include <dune/faultnetworks/utils/levelinterfacenetworkwriter.hh> #include <dune/faultnetworks/preconditioners/multilevelpatchpreconditioner.hh> +#include <dune/faultnetworks/preconditioners/genericmultilevelpreconditioner.hh> #include <dune/faultnetworks/solvers/osccgsolver.hh> #include <dune/faultnetworks/solvers/richardsonstep.hh> #include <dune/faultnetworks/solvers/cgstep.hh> @@ -177,6 +178,7 @@ int main(int argc, char** argv) { try // build multilevel cantor fault network CantorFaultFactory<GridType> faultFactory(levelToCantorLevel, coarseResolution, exactLevelIdx, maxCantorLevel); const InterfaceNetwork<GridType>& interfaceNetwork = faultFactory.interfaceNetwork(); + auto Ck = [](int k) { return std::pow(2.0, k); }; /* if (problemCount==0) { std::vector<int> writeLevels {1, 2}; @@ -220,7 +222,7 @@ int main(int argc, char** argv) { try std::vector<std::shared_ptr<LocalInterfaceAssembler>> coarseLocalInterfaceAssemblers(coarseLevelInterfaceNetwork.size()); for (size_t i=0; i<coarseLocalInterfaceAssemblers.size(); i++) { const int k = minCantorResolution + coarseLevelInterfaceNetwork.getIntersectionsLevels().at(i); - const double pen = (1+1.0/c) * std::pow(1.0+c, k) * std::pow(2.0, k-1); + const double pen = (1+1.0/c) * std::pow(1.0+c, k) * Ck(k); coarseLocalInterfaceAssemblers[i] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen); } @@ -266,7 +268,7 @@ int main(int argc, char** argv) { try std::vector<std::shared_ptr<LocalInterfaceAssembler>> exactLocalInterfaceAssemblers(exactLevelInterfaceNetwork.size()); for (size_t i=0; i<exactLocalInterfaceAssemblers.size(); i++) { const int k = minCantorResolution + exactLevelInterfaceNetwork.getIntersectionsLevels().at(i); - const double pen = (1+1.0/c) * std::pow(1.0+c, k) * std::pow(2.0, k-1); + const double pen = (1+1.0/c) * std::pow(1.0+c, k) * Ck(k); exactLocalInterfaceAssemblers[i] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen); } @@ -295,7 +297,7 @@ int main(int argc, char** argv) { try std::vector<std::shared_ptr<LocalInterfaceAssembler>> fineLocalInterfaceAssemblers(fineLevelInterfaceNetwork.size()); for (size_t i=0; i<fineLocalInterfaceAssemblers.size(); i++) { const int k = minCantorResolution + fineLevelInterfaceNetwork.getIntersectionsLevels().at(i); - const double pen = (1+1.0/c) * std::pow(1.0+c, k) * std::pow(2.0, k-1); + const double pen = (1+1.0/c) * std::pow(1.0+c, k) * Ck(k); fineLocalInterfaceAssemblers[i] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen); } @@ -326,7 +328,7 @@ int main(int argc, char** argv) { try levelLocalIntersectionAssemblers.resize(levelInterfaceNetwork.size()); for (size_t j=0; j<levelLocalIntersectionAssemblers.size(); j++) { const int k = minCantorResolution + levelInterfaceNetwork.getIntersectionsLevels().at(j); - const double pen = (1+1.0/c) * std::pow(1.0+c, k) * std::pow(2.0, k-1); + const double pen = (1+1.0/c) * std::pow(1.0+c, k) * Ck(k); levelLocalIntersectionAssemblers[j] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen); } } @@ -374,7 +376,7 @@ int main(int argc, char** argv) { try typedef DGMGTransfer<DGBasis, VectorType, BitVector, MatrixType> TransferImplementation; std::vector<std::shared_ptr<TransferImplementation>> transfer(0); - for (size_t i=0; i<activeLevels.size(); i++) { + /* for (size_t i=0; i<activeLevels.size(); i++) { if (activeLevels[i][0] && i+1<activeLevels.size()) { // init global problem on coarsest level const LevelInterfaceNetwork<GV>& coarseLevelNetwork = interfaceNetwork.levelInterfaceNetwork(i); @@ -385,7 +387,7 @@ int main(int argc, char** argv) { try auto t = std::make_shared<DGMGTransfer<DGBasis, VectorType, BitVector, MatrixType>>(coarseBasis, fineBasis); transfer.push_back(t); } - } + }*/ // set up smoothers and basesolver auto presmoother = @@ -419,14 +421,14 @@ int main(int argc, char** argv) { try // setup loop solver Dune::Solvers::LoopSolver<VectorType> - solver (mgStep, maxIterations, solverTolerance, + solver (cgStep, maxIterations, solverTolerance, fineEnergyNorm, Solver::FULL, true, &fineSol); if (useExact) { DGMGTransfer<DGBasis, VectorType, BitVector, MatrixType> discMGTransfer(*fineBasis, *exactBasis); - auto discAccuracyCriterion = discretizationAccuracy(mgStep, fineSol, exactSol, fineEnergyNorm, exactEnergyNorm, discMGTransfer); + auto discAccuracyCriterion = discretizationAccuracy(cgStep, fineSol, exactSol, fineEnergyNorm, exactEnergyNorm, discMGTransfer); solver.addCriterion(discAccuracyCriterion); } //solver.historyBuffer_ = outPath + "solutions/level_" + std::to_string(fineResolution) + ".vec"; diff --git a/src/cantorfaultnetworks/cantorfaultnetwork.parset b/src/cantorfaultnetworks/cantorfaultnetwork.parset index 0073f395969c81ada600d149d8614262675565ee..970d16d1ab1a2d5bbaeff2fbfaff1eb0962ea616 100644 --- a/src/cantorfaultnetworks/cantorfaultnetwork.parset +++ b/src/cantorfaultnetworks/cantorfaultnetwork.parset @@ -1,9 +1,9 @@ path = /home/joscha/software/dune/dune-faultnetworks/src/data/ -resultPath = /home/joscha/software/dune/dune-faultnetworks/src/cantorfaultnetworks/results/cantor/mg/ +resultPath = /home/joscha/software/dune/dune-faultnetworks/src/cantorfaultnetworks/results/cantor/cell-GS [preconditioner] -patch = SUPPORT # CELL , SUPPORT -mode = MULTIPLICATIVE # ADDITIVE , MULTIPLICATIVE +patch = CELL # CELL , SUPPORT +mode = MULTIPLICATIVE # ADDITIVE , MULTIPLICATIVE multDirection = SYMMETRIC # SYMMETRIC , FORWARD , BACKWARD patchDepth = 1 @@ -14,8 +14,8 @@ oscDataFile = oscDataLaplace4.mat # level resolution in 2^(-...) coarseResolution = 0 -fineResolution = 3 -exactResolution = 4 +fineResolution = 9 +exactResolution = 10 minCantorResolution = 0 penaltyFactor = 1 @@ -23,6 +23,61 @@ penaltyFactor = 1 # cg solver nestedIteration = 0 tol = 1e-12 -maxIterations = 8 +maxIterations = 10 ########################################### + +[problem1] +oscDataFile = oscDataLaplace4.mat + +# level resolution in 2^(-...) +coarseResolution = 0 +fineResolution = 9 +exactResolution = 10 +minCantorResolution = 0 + +penaltyFactor = 1 + +# cg solver +nestedIteration = 0 +tol = 1e-12 +maxIterations = 10 + +########################################### + +[problem2] +oscDataFile = oscDataLaplace4.mat + +# level resolution in 2^(-...) +coarseResolution = 0 +fineResolution = 9 +exactResolution = 10 +minCantorResolution = 0 + +penaltyFactor = 1 + +# cg solver +nestedIteration = 0 +tol = 1e-12 +maxIterations = 10 + +########################################### + +[problem3] +oscDataFile = oscDataLaplace4.mat + +# level resolution in 2^(-...) +coarseResolution = 0 +fineResolution = 9 +exactResolution = 10 +minCantorResolution = 0 + +penaltyFactor = 1 + +# cg solver +nestedIteration = 0 +tol = 1e-12 +maxIterations = 10 + +########################################### + diff --git a/src/cantorfaultnetworks/sparsecantorfaultnetwork.cc b/src/cantorfaultnetworks/sparsecantorfaultnetwork.cc index 932680b5f4937bb684d3af2debc6aedfa717ed5f..07c1c85434c98f4f8dbb4598314228d202e9a97c 100644 --- a/src/cantorfaultnetworks/sparsecantorfaultnetwork.cc +++ b/src/cantorfaultnetworks/sparsecantorfaultnetwork.cc @@ -178,8 +178,10 @@ int main(int argc, char** argv) { try // build multilevel cantor fault network SparseCantorFaultFactory<GridType> faultFactory(minCantorLevel, std::max(fineCantorLevel, maxCantorLevel)); const auto& interfaceNetwork = faultFactory.interfaceNetwork(); + auto Ck = [](int k) { return std::pow(2.0, k+1); }; - /*if (problemCount==0) { + + if (problemCount==0) { std::vector<int> writeLevels(0); for (size_t i=0; i<=exactLevelIdx; i++, i++) { writeLevels.push_back(i); @@ -193,7 +195,7 @@ int main(int argc, char** argv) { try LevelInterfaceNetworkWriter networkWriter(outPath + "levelinterfacenetwork_" + std::to_string(writeLevel) + ".tikz"); networkWriter.write(interfaceNetwork.levelInterfaceNetwork(writeLevel), writeGrid); } - }*/ + } const GridType& grid = faultFactory.grid(); @@ -226,7 +228,7 @@ int main(int argc, char** argv) { try std::vector<std::shared_ptr<LocalInterfaceAssembler>> coarseLocalInterfaceAssemblers(coarseLevelInterfaceNetwork.size()); for (size_t i=0; i<coarseLocalInterfaceAssemblers.size(); i++) { const int k = minCantorLevel + coarseLevelInterfaceNetwork.getIntersectionsLevels().at(i)/2; - const double pen = (1+1.0/c) * std::pow(1.0+c, k) * std::pow(2.0, k); + const double pen = (1+1.0/c) * std::pow(1.0+c, k) * Ck(k); coarseLocalInterfaceAssemblers[i] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen); } @@ -273,7 +275,7 @@ int main(int argc, char** argv) { try std::vector<std::shared_ptr<LocalInterfaceAssembler>> exactLocalInterfaceAssemblers(exactLevelInterfaceNetwork.size()); for (size_t i=0; i<exactLocalInterfaceAssemblers.size(); i++) { const int k = minCantorLevel + exactLevelInterfaceNetwork.getIntersectionsLevels().at(i)/2; - const double pen = (1+1.0/c) * std::pow(1.0+c, k) * std::pow(2.0, k); + const double pen = (1+1.0/c) * std::pow(1.0+c, k) * Ck(k); exactLocalInterfaceAssemblers[i] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen); } @@ -302,7 +304,7 @@ int main(int argc, char** argv) { try std::vector<std::shared_ptr<LocalInterfaceAssembler>> fineLocalInterfaceAssemblers(fineLevelInterfaceNetwork.size()); for (size_t i=0; i<fineLocalInterfaceAssemblers.size(); i++) { const int k = minCantorLevel + fineLevelInterfaceNetwork.getIntersectionsLevels().at(i)/2; - const double pen = (1+1.0/c) * std::pow(1.0+c, k) * std::pow(2.0, k); + const double pen = (1+1.0/c) * std::pow(1.0+c, k) * Ck(k); fineLocalInterfaceAssemblers[i] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen); } @@ -338,7 +340,7 @@ int main(int argc, char** argv) { try levelLocalIntersectionAssemblers.resize(levelInterfaceNetwork.size()); for (size_t j=0; j<levelLocalIntersectionAssemblers.size(); j++) { const int k = minCantorLevel + levelInterfaceNetwork.getIntersectionsLevels().at(j)/2; - const double pen = (1+1.0/c) * std::pow(1.0+c, k) * std::pow(2.0, k); + const double pen = (1+1.0/c) * std::pow(1.0+c, k) * Ck(k); levelLocalIntersectionAssemblers[j] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen); } } diff --git a/src/cantorfaultnetworks/sparsecantorfaultnetwork.parset b/src/cantorfaultnetworks/sparsecantorfaultnetwork.parset index 10255cbb4e1c69daf5efa40dd35a194350152fa9..ba932e801c605890f46cb189b4370d4133661805 100644 --- a/src/cantorfaultnetworks/sparsecantorfaultnetwork.parset +++ b/src/cantorfaultnetworks/sparsecantorfaultnetwork.parset @@ -1,9 +1,9 @@ path = /home/joscha/software/dune/dune-faultnetworks/src/data/ -resultPath = /home/joscha/software/dune/dune-faultnetworks/src/cantorfaultnetworks/results/sparse/cell/add/ +resultPath = /home/joscha/software/dune/dune-faultnetworks/src/cantorfaultnetworks/results/sparse/patch-GS/ [preconditioner] -patch = CELL # CELL , SUPPORT -mode = ADDITIVE # ADDITIVE , MULTIPLICATIVE +patch = SUPPORT # CELL , SUPPORT +mode = MULTIPLICATIVE # ADDITIVE , MULTIPLICATIVE multDirection = SYMMETRIC # SYMMETRIC , FORWARD , BACKWARD patchDepth = 1 @@ -15,15 +15,16 @@ oscDataFile = oscDataLaplace4.mat # level resolution in 2^(-...) coarseCantorLevel = 1 fineCantorLevel = 2 -#maxCantorLevel = 3 +maxCantorLevel = 3 penaltyFactor = 1 # cg solver nestedIteration = 0 -tol = 1e-12 +tol = 1e-15 maxIterations = 10 +########################################### [problem1] oscDataFile = oscDataLaplace4.mat @@ -31,15 +32,16 @@ oscDataFile = oscDataLaplace4.mat # level resolution in 2^(-...) coarseCantorLevel = 1 fineCantorLevel = 3 -#maxCantorLevel = 3 +maxCantorLevel = 4 penaltyFactor = 1 # cg solver nestedIteration = 0 -tol = 1e-12 +tol = 1e-15 maxIterations = 10 +########################################### [problem2] oscDataFile = oscDataLaplace4.mat @@ -47,15 +49,16 @@ oscDataFile = oscDataLaplace4.mat # level resolution in 2^(-...) coarseCantorLevel = 1 fineCantorLevel = 4 -#maxCantorLevel = 3 +maxCantorLevel = 5 penaltyFactor = 1 # cg solver nestedIteration = 0 -tol = 1e-12 +tol = 1e-15 maxIterations = 10 +########################################### [problem3] oscDataFile = oscDataLaplace4.mat @@ -63,11 +66,11 @@ oscDataFile = oscDataLaplace4.mat # level resolution in 2^(-...) coarseCantorLevel = 1 fineCantorLevel = 5 -#maxCantorLevel = 3 +maxCantorLevel = 5 penaltyFactor = 1 # cg solver nestedIteration = 0 -tol = 1e-12 +tol = 1e-15 maxIterations = 10 diff --git a/src/geofaultnetworks/rockfaultnetwork.cc b/src/geofaultnetworks/rockfaultnetwork.cc index 4ab84e118763a48fb2525d549094e3a717cce8c6..03e7f294eed79c7017fb45a76a21578c21e862ed 100644 --- a/src/geofaultnetworks/rockfaultnetwork.cc +++ b/src/geofaultnetworks/rockfaultnetwork.cc @@ -35,6 +35,10 @@ // dune-solver includes #include <dune/solvers/norms/energynorm.hh> #include <dune/solvers/norms/twonorm.hh> +#include <dune/solvers/solvers/loopsolver.hh> +#include <dune/solvers/transferoperators/compressedmultigridtransfer.hh> +#include <dune/solvers/iterationsteps/blockgssteps.hh> +#include <dune/solvers/iterationsteps/multigridstep.hh> // dune-faultnetworks includes #include <dune/faultnetworks/utils/debugutils.hh> @@ -43,7 +47,11 @@ #include <dune/faultnetworks/utils/levelinterfacenetworkwriter.hh> #include <dune/faultnetworks/preconditioners/multilevelpatchpreconditioner.hh> +#include <dune/faultnetworks/preconditioners/genericmultilevelpreconditioner.hh> #include <dune/faultnetworks/solvers/osccgsolver.hh> +#include <dune/faultnetworks/solvers/richardsonstep.hh> +#include <dune/faultnetworks/solvers/cgstep.hh> +#include <dune/faultnetworks/solvers/criterion.hh> #include <dune/faultnetworks/assemblers/intersectionjumpmassassembler.hh> #include <dune/faultnetworks/assemblers/osclocalassembler.hh> @@ -60,9 +68,14 @@ #include <dune/faultnetworks/faultfactories/rockfaultfactory.hh> -#include <dune/faultnetworks/solvers/richardsonstep.hh> -#include <dune/faultnetworks/solvers/criterion.hh> -#include <dune/solvers/solvers/loopsolver.hh> +const std::string sourcePath = "/home/joscha/software/dune/dune-faultnetworks/src/geofaultnetworks/"; + +Dune::ParameterTree getParameters(int argc, char *argv[]) { + Dune::ParameterTree parset; + Dune::ParameterTreeParser::readINITree(sourcePath + "rockfaultnetwork.parset", parset); + Dune::ParameterTreeParser::readOptions(argc, argv, parset); + return parset; +} int main(int argc, char** argv) { try { @@ -91,43 +104,47 @@ int main(int argc, char** argv) { try typedef OscLocalAssembler<GridType, DGFE, DGFE> LocalOscAssembler; typedef IntersectionJumpMassAssembler<GV, DGFE, DGFE> LocalInterfaceAssembler; - // parse parameter file - ParameterTree parameterSet; - if (argc==2) - ParameterTreeParser::readINITree(argv[1], parameterSet); - else - ParameterTreeParser::readINITree("rockfaultnetwork.parset", parameterSet); + // parse parameter file + auto const parset = getParameters(argc, argv); - const std::string path = parameterSet.get<std::string>("path"); - const std::string resultPath = parameterSet.get<std::string>("resultPath"); + // set io + const auto dataPath = parset.get<std::string>("path"); + const auto outPath = parset.get<std::string>("resultPath"); int problemCount = 0; - while (parameterSet.hasSub("problem" + std::to_string(problemCount))) { + while (parset.hasSub("problem" + std::to_string(problemCount))) { - ParameterTree& problemParameters = parameterSet.sub("problem" + std::to_string(problemCount)); + auto problemParset = parset.sub("problem" + std::to_string(problemCount)); + + const std::string oscDataFile = problemParset.get<std::string>("oscDataFile"); - const std::string oscDataFile = problemParameters.get<std::string>("oscDataFile"); + const int coarseResolution = problemParset.get<int>("coarseResolution"); + const int fineResolution = problemParset.get<int>("fineResolution"); - const int coarseResolution = problemParameters.get<int>("coarseResolution"); - const int fineResolution = problemParameters.get<int>("fineResolution"); - const int exactResolution = problemParameters.get<int>("exactResolution"); + bool useExact = problemParset.hasKey("exactResolution"); + const int exactResolution = useExact ? problemParset.get<int>("exactResolution") : fineResolution; - const double c = problemParameters.get<double>("penaltyFactor"); + const double c = problemParset.get<double>("penaltyFactor"); - const int maxIterations = problemParameters.get<int>("maxIterations"); - const double solverTolerance = problemParameters.get<double>("tol"); - const bool nestedIteration = problemParameters.get<bool>("nestedIteration"); + const int maxIterations = problemParset.get<int>("maxIterations"); + const double solverTolerance = problemParset.get<double>("tol"); + const bool nestedIteration = problemParset.get<bool>("nestedIteration"); const int fineGridN = std::pow(2, fineResolution); - std::ofstream out(resultPath + oscDataFile + "_" + std::to_string(coarseResolution) + "_" + std::to_string(fineResolution) + "_" + std::to_string(exactResolution) + ".log"); + std::ofstream out(outPath + oscDataFile + "_" + std::to_string(coarseResolution) + "_" + std::to_string(fineResolution) + "_" + std::to_string(exactResolution) + ".log"); std::streambuf *coutbuf = std::cout.rdbuf(); //save old buf std::cout.rdbuf(out.rdbuf()); //redirect std::cout to out.txt - std::cout << std::endl; std::cout << "Parameter file read successfully!" << std::endl; + std::cout << "preconditioner:" << std::endl; + std::cout << "patch: " << parset.get<std::string>("preconditioner.patch") << std::endl; + std::cout << "mode: " << parset.get<std::string>("preconditioner.mode") << std::endl; + std::cout << "multDirection: " << parset.get<std::string>("preconditioner.multDirection") << std::endl; + std::cout << "-----------------------------" << std::endl; + typedef Dune::Matrix<Dune::FieldMatrix<double,1,1>, std::allocator<Dune::FieldMatrix<double,1,1>>> OscMatrixType; OscMatrixType matrix; @@ -138,7 +155,7 @@ int main(int argc, char** argv) { try B[0][1] = -1; B[1][0] = -1; - MatrixReader<OscMatrixType> matrixReader(path + oscDataFile); + MatrixReader<OscMatrixType> matrixReader(dataPath + oscDataFile); matrixReader.read(matrix); const int oscGridN = matrix.N(); @@ -153,7 +170,7 @@ int main(int argc, char** argv) { try const int exactLevelIdx = exactResolution - coarseResolution; // build multilevel fault network - const auto& networkParset = parameterSet.sub("network"); + const auto& networkParset = parset.sub("network"); const unsigned int randomSeed = networkParset.get<unsigned int>("randomSeed"); const double maxAngle = networkParset.get<double>("maxAngle"); @@ -162,10 +179,11 @@ int main(int argc, char** argv) { try std::srand(randomSeed); RockFaultFactory<GridType> faultFactory(coarseResolution, exactLevelIdx, maxAngle); const auto& interfaceNetwork = faultFactory.interfaceNetwork(); + auto Ck = [](int k) { return std::pow(2.0, k); }; std::cout << "done!" << std::endl; - if (problemCount==0) { + /*if (problemCount==0) { std::vector<int> writeLevels(exactLevelIdx+1); // {0, 1, 2}; for (size_t i=0; i<writeLevels.size(); i++) { writeLevels[i] = i; @@ -175,10 +193,10 @@ int main(int argc, char** argv) { try int writeLevel = writeLevels[i]; bool writeGrid = false; //coarseResolution+writeLevel<7; - LevelInterfaceNetworkWriter networkWriter(resultPath + "levelinterfacenetwork_" + std::to_string(writeLevel) + ".tikz"); + LevelInterfaceNetworkWriter networkWriter(outPath + "levelinterfacenetwork_" + std::to_string(writeLevel) + ".tikz"); networkWriter.write(interfaceNetwork.levelInterfaceNetwork(writeLevel), writeGrid); } - } + }*/ const GridType& grid = faultFactory.grid(); @@ -211,7 +229,7 @@ int main(int argc, char** argv) { try std::vector<std::shared_ptr<LocalInterfaceAssembler>> coarseLocalInterfaceAssemblers(coarseLevelInterfaceNetwork.size()); for (size_t i=0; i<coarseLocalInterfaceAssemblers.size(); i++) { const int k = coarseResolution + coarseLevelInterfaceNetwork.getIntersectionsLevels().at(i); - const double pen = (1+1.0/c) * std::pow(1.0+c, k) * std::pow(2.0, k); + const double pen = (1+1.0/c) * std::pow(1.0+c, k) * Ck(k); coarseLocalInterfaceAssemblers[i] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen); } @@ -228,7 +246,7 @@ int main(int argc, char** argv) { try coarseSol.resize(coarseBasis->size()); std::stringstream solFilename; - solFilename << resultPath << "solutions/level_" << std::to_string(fineLevelIdx-1); + solFilename << outPath << "solutions/level_" << std::to_string(fineLevelIdx-1); std::ifstream file(solFilename.str().c_str(), std::ios::in|std::ios::binary); if (not(file)) @@ -242,27 +260,35 @@ int main(int argc, char** argv) { try // ---------------------------------------------------------------- // --- compute exact solution --- // ---------------------------------------------------------------- - std::cout << std::endl; - std::cout << "Computing exact solution!" << std::endl; + VectorType exactSol; + std::shared_ptr<DGBasis> exactBasis = nullptr; + EnergyNorm<MatrixType, VectorType> exactEnergyNorm; const LevelInterfaceNetwork<GV>& exactLevelInterfaceNetwork = interfaceNetwork.levelInterfaceNetwork(exactLevelIdx); LevelInterfaceNetworkProblem<DGBasis, MatrixType, VectorType> exactGlobalAssembler(exactLevelInterfaceNetwork); - LocalOscAssembler exactLocalAssembler(oscDataHierarchy[exactLevelIdx]->data(), oscDataHierarchy[exactLevelIdx]->mapper()); - std::vector<std::shared_ptr<LocalInterfaceAssembler>> exactLocalInterfaceAssemblers(exactLevelInterfaceNetwork.size()); - for (size_t i=0; i<exactLocalInterfaceAssemblers.size(); i++) { - const int k = coarseResolution + exactLevelInterfaceNetwork.getIntersectionsLevels().at(i); - const double pen = (1+1.0/c) * std::pow(1.0+c, k) * std::pow(2.0, k); - exactLocalInterfaceAssemblers[i] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen); - } + if (useExact) { + std::cout << std::endl; + std::cout << "Computing exact solution!" << std::endl; + + + LocalOscAssembler exactLocalAssembler(oscDataHierarchy[exactLevelIdx]->data(), oscDataHierarchy[exactLevelIdx]->mapper()); + + std::vector<std::shared_ptr<LocalInterfaceAssembler>> exactLocalInterfaceAssemblers(exactLevelInterfaceNetwork.size()); + for (size_t i=0; i<exactLocalInterfaceAssemblers.size(); i++) { + const int k = coarseResolution + exactLevelInterfaceNetwork.getIntersectionsLevels().at(i); + const double pen = (1+1.0/c) * std::pow(1.0+c, k) * Ck(k); + exactLocalInterfaceAssemblers[i] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen); + } - exactGlobalAssembler.assemble(exactLocalAssembler, exactLocalInterfaceAssemblers, l2FunctionalAssembler); - exactGlobalAssembler.solve(); + exactGlobalAssembler.assemble(exactLocalAssembler, exactLocalInterfaceAssemblers, l2FunctionalAssembler); + exactGlobalAssembler.solve(); - const VectorType& exactSol = exactGlobalAssembler.getSol(); - std::shared_ptr<DGBasis> exactBasis = exactGlobalAssembler.basis(); + exactSol = exactGlobalAssembler.getSol(); + exactBasis = exactGlobalAssembler.basis(); - EnergyNorm<MatrixType, VectorType> exactEnergyNorm(exactGlobalAssembler.matrix()); + exactEnergyNorm.setMatrix(&exactGlobalAssembler.matrix()); + } // ---------------------------------------------------------------- // --- set up cg solver --- @@ -280,7 +306,7 @@ int main(int argc, char** argv) { try std::vector<std::shared_ptr<LocalInterfaceAssembler>> fineLocalInterfaceAssemblers(fineLevelInterfaceNetwork.size()); for (size_t i=0; i<fineLocalInterfaceAssemblers.size(); i++) { const int k = coarseResolution + fineLevelInterfaceNetwork.getIntersectionsLevels().at(i); - const double pen = (1+1.0/c) * std::pow(1.0+c, k) * std::pow(2.0, k); + const double pen = (1+1.0/c) * std::pow(1.0+c, k) * Ck(k); fineLocalInterfaceAssemblers[i] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen); } @@ -295,9 +321,9 @@ int main(int argc, char** argv) { try std::cout << "Setting up the preconditioner!" << std::endl; - BitVector activeLevels(fineLevelIdx+1, false); - activeLevels[0] = true; - activeLevels[fineLevelIdx] = true; + BitVector activeLevels(fineLevelIdx+1, true); + //activeLevels[0] = true; + //activeLevels[fineLevelIdx] = true; std::vector<std::shared_ptr<LocalOscAssembler>> localAssemblers; std::vector<std::vector<std::shared_ptr<LocalInterfaceAssembler>>> localIntersectionAssemblers; @@ -313,15 +339,15 @@ int main(int argc, char** argv) { try levelLocalIntersectionAssemblers.resize(levelInterfaceNetwork.size()); for (size_t j=0; j<levelLocalIntersectionAssemblers.size(); j++) { const int k = coarseResolution + levelInterfaceNetwork.getIntersectionsLevels().at(j); - const double pen = (1+1.0/c) * std::pow(1.0+c, k) * std::pow(2.0, k); + const double pen = (1+1.0/c) * std::pow(1.0+c, k) * Ck(k); levelLocalIntersectionAssemblers[j] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen); } } - const auto& preconditionerParset = parameterSet.sub("preconditioner"); - using MultilevelPatchPreconditioner = MultilevelPatchPreconditioner<DGBasis, LocalOscAssembler, LocalInterfaceAssembler, MatrixType, VectorType>; + const auto& preconditionerParset = parset.sub("preconditioner"); + using Preconditioner = MultilevelPatchPreconditioner<DGBasis, LocalOscAssembler, LocalInterfaceAssembler, MatrixType, VectorType>; - MultilevelPatchPreconditioner preconditioner(interfaceNetwork, activeLevels, localAssemblers, localIntersectionAssemblers, preconditionerParset); + Preconditioner preconditioner(interfaceNetwork, activeLevels, localAssemblers, localIntersectionAssemblers, preconditionerParset); /*for (size_t i=0; i<preconditioner.size()-1; i++) { auto& levelFaultPreconditioner = *preconditioner.levelPatchPreconditioner(i); @@ -335,40 +361,92 @@ int main(int argc, char** argv) { try // set initial iterate VectorType initialX; - DGMGTransfer<DGBasis> coarseFineTransfer(*coarseBasis, *fineBasis); + DGMGTransfer<DGBasis, VectorType, BitVector, MatrixType> coarseFineTransfer(*coarseBasis, *fineBasis); coarseFineTransfer.prolong(coarseSol, initialX); VectorType rhs(fineRhs); - DGMGTransfer<DGBasis> discMGTransfer(*fineBasis, *exactBasis); - - // solve - /*OscCGSolver<MatrixType, VectorType, DGBasis> - solver(*fineBasis, &fineGlobalAssembler.matrix(), &initialX, &rhs, &preconditioner, - maxIterations, solverTolerance, &exactEnergyNorm, - Solver::FULL, &exactSol, &fineSol, discMGTransfer, 1.0, true); //((oscGridN+0.0)/fineGridN) */ VectorType x = initialX; preconditioner.setProblem(fineGlobalAssembler.matrix(), x, rhs); + + const auto& gridView = fineLevelInterfaceNetwork.levelGridView(); + + // set boundary conditions + Dune::BitSetVector<1> boundaryDofs; + BoundaryPatch<GV> boundaryPatch(gridView, true); + constructBoundaryDofs(boundaryPatch, *fineBasis, boundaryDofs); + + typedef MultigridStep<MatrixType, VectorType, BitVector> MGStep; + typedef MultigridTransfer<VectorType, BitVector, MatrixType> Transfer; + typedef DGMGTransfer<DGBasis, VectorType, BitVector, MatrixType> TransferImplementation; + + std::vector<std::shared_ptr<TransferImplementation>> transfer(0); + /* for (size_t i=0; i<activeLevels.size(); i++) { + if (activeLevels[i][0] && i+1<activeLevels.size()) { + // init global problem on coarsest level + const LevelInterfaceNetwork<GV>& coarseLevelNetwork = interfaceNetwork.levelInterfaceNetwork(i); + const LevelInterfaceNetwork<GV>& fineLevelNetwork = interfaceNetwork.levelInterfaceNetwork(i+1); + + DGBasis coarseBasis(coarseLevelNetwork); + DGBasis fineBasis(fineLevelNetwork); + auto t = std::make_shared<DGMGTransfer<DGBasis, VectorType, BitVector, MatrixType>>(coarseBasis, fineBasis); + transfer.push_back(t); + } + }*/ + + // set up smoothers and basesolver + auto presmoother = + Dune::Solvers::BlockGSStepFactory<MatrixType, VectorType>::create( + Dune::Solvers::BlockGS::LocalSolvers::gs(), Dune::Solvers::BlockGS::Direction::BACKWARD); + auto postsmoother = + Dune::Solvers::BlockGSStepFactory<MatrixType, VectorType>::create( + Dune::Solvers::BlockGS::LocalSolvers::gs(), Dune::Solvers::BlockGS::Direction::FORWARD); + + auto basesolver_step = + Dune::Solvers::BlockGSStepFactory<MatrixType, VectorType>::create( + Dune::Solvers::BlockGS::LocalSolvers::gs()); + + EnergyNorm<MatrixType, VectorType> basenorm(basesolver_step); + + Dune::Solvers::LoopSolver<VectorType> basesolver(basesolver_step, 1000, 1e-10, basenorm, Solver::QUIET); VectorType x1 = x; - auto richardsonStep = RichardsonStep<MatrixType, VectorType>(preconditioner, 1.0); - richardsonStep.setProblem(fineGlobalAssembler.matrix(), x1, rhs); + + // create iteration step and set all kinds of data + MGStep mgStep; + mgStep.setTransferOperators(transfer); + mgStep.setSmoother(presmoother, postsmoother); + mgStep.setMGType(1,3,3); + mgStep.setIgnore(boundaryDofs); + mgStep.setBaseSolver(basesolver); + mgStep.setProblem(fineGlobalAssembler.matrix(), x1, rhs); + + auto cgStep = CGStep<MatrixType, VectorType>(); //RichardsonStep<MatrixType, VectorType>(preconditioner, 1.0); //CGStep<MatrixType, VectorType>(); + cgStep.setProblem(fineGlobalAssembler.matrix(), x1, rhs); + cgStep.setPreconditioner(preconditioner); // setup loop solver Dune::Solvers::LoopSolver<VectorType> - solver (richardsonStep, maxIterations, solverTolerance, + solver (cgStep, maxIterations, solverTolerance, fineEnergyNorm, Solver::FULL, true, &fineSol); - auto discAccuracyCriterion = discretizationAccuracy(richardsonStep, fineSol, exactSol, fineEnergyNorm, exactEnergyNorm, discMGTransfer); - solver.addCriterion(discAccuracyCriterion); + if (useExact) { + DGMGTransfer<DGBasis, VectorType, BitVector, MatrixType> discMGTransfer(*fineBasis, *exactBasis); + auto discAccuracyCriterion = discretizationAccuracy(cgStep, fineSol, exactSol, fineEnergyNorm, exactEnergyNorm, discMGTransfer); + solver.addCriterion(discAccuracyCriterion); + } + //solver.historyBuffer_ = resultPath + "solutions/level_" + std::to_string(fineCantorLevel) + ".vec"; solver.check(); solver.preprocess(); solver.solve(); + auto& res = solver.getResult(); + std::cout << "convRate: " << std::setprecision(3) << std::pow(res.conv_rate, 1.0/(res.iterations-1.0)) << std::endl; + std::cout << std::endl; std::cout << "Total time: " << totalTimer.elapsed() << " seconds." << std::endl; diff --git a/src/geofaultnetworks/rockfaultnetwork.parset b/src/geofaultnetworks/rockfaultnetwork.parset index 725219196e6c366ab6eb9fa6bafbfa7753b62b62..834e6a609fc83f19d689c17ce061dc397555fc33 100644 --- a/src/geofaultnetworks/rockfaultnetwork.parset +++ b/src/geofaultnetworks/rockfaultnetwork.parset @@ -1,14 +1,14 @@ -path = ../data/ -resultPath = ../geofaultnetworks/results/rock/ +path = /home/joscha/software/dune/dune-faultnetworks/src/data/ +resultPath = /home/joscha/software/dune/dune-faultnetworks/src/geofaultnetworks/results/rock/cell/mult/ [network] -randomSeed = 15 +randomSeed = 14 maxAngle = 2.0 [preconditioner] patch = CELL # CELL , SUPPORT mode = MULTIPLICATIVE # ADDITIVE , MULTIPLICATIVE -multDirection = BACKWARD # SYMMETRIC , FORWARD , BACKWARD +multDirection = SYMMETRIC # SYMMETRIC , FORWARD , BACKWARD patchDepth = 1 ########################################### @@ -30,70 +30,3 @@ maxIterations = 10 ########################################### -[problem1] -oscDataFile = oscDataLaplace4.mat - -# level resolution in 2^(-...) -coarseResolution = 4 -fineResolution = 8 -exactResolution = 9 - -penaltyFactor = 1 - -# cg solver -nestedIteration = 0 -tol = 1e-12 -maxIterations = 10 - -########################################### - -[problem2] -oscDataFile = oscDataLaplace4.mat - -# level resolution in 2^(-...) -coarseResolution = 4 -fineResolution = 7 -exactResolution = 8 - -penaltyFactor = 1 - -# cg solver -nestedIteration = 0 -tol = 1e-12 -maxIterations = 10 - -########################################### - -[problem3] -oscDataFile = oscDataLaplace4.mat - -# level resolution in 2^(-...) -coarseResolution = 4 -fineResolution = 6 -exactResolution = 7 - -penaltyFactor = 1 - -# cg solver -nestedIteration = 0 -tol = 1e-12 -maxIterations = 10 - -########################################### - -[problem4] -oscDataFile = oscDataLaplace4.mat - -# level resolution in 2^(-...) -coarseResolution = 4 -fineResolution = 5 -exactResolution = 6 - -penaltyFactor = 1 - -# cg solver -nestedIteration = 0 -tol = 1e-12 -maxIterations = 10 - -###########################################