Code owners
Assign users and groups as approvers for specific file changes. Learn more.
sparsecantorfaultnetwork.cc 14.79 KiB
#include <config.h>
#include <cmath>
#include <iostream>
#include <sstream>
#include <fstream>
#include <string>
/*
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wunused-parameter"
#include <dune/grid/uggrid/ugwrapper.hh>
#pragma GCC diagnostic pop
*/
// dune-common includes
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/shared_ptr.hh>
#include <dune/common/parametertree.hh>
#include <dune/common/parametertreeparser.hh>
#include <dune/common/timer.hh>
// dune-istl includes
#include <dune/istl/matrix.hh>
#include <dune/istl/bcrsmatrix.hh>
// dune-fufem includes
#include <dune/fufem/assemblers/localassemblers/l2functionalassembler.hh>
// dune-grid includes
#include <dune/grid/uggrid.hh>
#include <dune/grid/io/file/vtk/vtkwriter.hh>
// dune-solver includes
#include <dune/solvers/norms/energynorm.hh>
#include <dune/solvers/norms/twonorm.hh>
// dune-faultnetworks includes
#include <dune/faultnetworks/utils/debugutils.hh>
#include <dune/faultnetworks/utils/matrixreader.hh>
#include <dune/faultnetworks/utils/vectorreader.hh>
#include <dune/faultnetworks/utils/levelinterfacenetworkwriter.hh>
#include <dune/faultnetworks/preconditioners/multilevelpatchpreconditioner.hh>
#include <dune/faultnetworks/solvers/osccgsolver.hh>
#include <dune/faultnetworks/assemblers/intersectionjumpmassassembler.hh>
#include <dune/faultnetworks/assemblers/osclocalassembler.hh>
#include <dune/faultnetworks/oscrhs.hh>
#include <dune/faultnetworks/oscdata.hh>
#include <dune/faultnetworks/oscdatahierarchy.hh>
#include <dune/faultnetworks/faultp1nodalbasis.hh>
#include <dune/faultnetworks/interfacenetwork.hh>
#include <dune/faultnetworks/levelinterfacenetwork.hh>
#include <dune/faultnetworks/levelinterfacenetworkproblem.hh>
#include <dune/faultnetworks/dgmgtransfer.hh>
#include <dune/faultnetworks/faultfactories/sparsecantorfaultfactory.hh>
int main(int argc, char** argv) { try
{
MPIHelper::instance(argc, argv);
const int dim = 2;
//typedef double ctype;
typedef typename Dune::FieldVector<double, dim> WorldVectorType;
typedef typename Dune::BlockVector<Dune::FieldVector<double,1>> VectorType;
typedef typename Dune::BitSetVector<1> BitVector;
typedef typename Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>> MatrixType;
#if HAVE_UG
typedef typename Dune::UGGrid<dim> GridType;
#else
#error No UG!
#endif
typedef typename GridType::LevelGridView GV;
typedef FaultP1NodalBasis<GV, double> DGBasis;
typedef typename DGBasis::LocalFiniteElement DGFE;
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("sparsecantorfaultnetwork.parset", parameterSet);
const std::string path = parameterSet.get<std::string>("path");
const std::string resultPath = parameterSet.get<std::string>("resultPath");
int problemCount = 0;
while (parameterSet.hasSub("problem" + std::to_string(problemCount))) {
ParameterTree& problemParameters = parameterSet.sub("problem" + std::to_string(problemCount));
const std::string oscDataFile = problemParameters.get<std::string>("oscDataFile");
const size_t minCantorLevel = problemParameters.get<size_t>("coarseCantorLevel");
const size_t fineCantorLevel = problemParameters.get<size_t>("fineCantorLevel");
const size_t maxCantorLevel = problemParameters.get<size_t>("maxCantorLevel");
const double c = problemParameters.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 fineGridN = std::pow(4, fineCantorLevel);
std::ofstream out(resultPath + oscDataFile + "_" + std::to_string(minCantorLevel) + "_" + std::to_string(fineCantorLevel) + "_" + std::to_string(maxCantorLevel) + ".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;
typedef Dune::Matrix<Dune::FieldMatrix<double,1,1>, std::allocator<Dune::FieldMatrix<double,1,1>>> OscMatrixType;
OscMatrixType matrix;
// interface integral jump penalty, B:
OscMatrixType B(2,2);
B[0][0] = 1;
B[1][1] = 1;
B[0][1] = -1;
B[1][0] = -1;
MatrixReader<OscMatrixType> matrixReader(path + oscDataFile);
matrixReader.read(matrix);
const int oscGridN = matrix.N(); const int oscGridLevelIdx = std::round(std::log2(oscGridN)) - minCantorLevel;
if (oscGridN>fineGridN)
DUNE_THROW(Dune::Exception, "Provided oscData too large!");
std::cout << "OscData file read successfully!" << std::endl << std::endl;
const int fineLevelIdx = 2*(fineCantorLevel - minCantorLevel);
const int exactLevelIdx = 2*(maxCantorLevel - minCantorLevel);
// build multilevel cantor fault network
SparseCantorFaultFactory<GridType> faultFactory(minCantorLevel, maxCantorLevel);
const auto& interfaceNetwork = faultFactory.interfaceNetwork();
if (problemCount==0) {
std::vector<int> writeLevels {0, 2};
for (size_t i=0; i<writeLevels.size(); i++) {
int writeLevel = writeLevels[i];
bool writeGrid = !(writeLevel==8);
LevelInterfaceNetworkWriter networkWriter(resultPath + "levelinterfacenetwork_" + std::to_string(writeLevel) + ".tikz");
networkWriter.write(interfaceNetwork.levelInterfaceNetwork(writeLevel), writeGrid);
}
}
const GridType& grid = faultFactory.grid();
OscData<GridType, OscMatrixType> oscData(grid, oscGridLevelIdx);
oscData.set(matrix);
OscDataHierarchy<GridType, OscMatrixType> oscDataHierarchy(oscData);
Timer totalTimer;
totalTimer.reset();
totalTimer.start();
// ----------------------------------------------------------------
// --- compute initial iterate ---
// ----------------------------------------------------------------
VectorType coarseSol;
std::shared_ptr<DGBasis> coarseBasis;
OscRhs<WorldVectorType, Dune::FieldVector<double,1>> f;
L2FunctionalAssembler<GridType, DGFE> l2FunctionalAssembler(f);
std::cout << std::endl;
std::cout << "Computing initial iterate!" << std::endl;
if (!nestedIteration) {
const LevelInterfaceNetwork<GV>& coarseLevelInterfaceNetwork = interfaceNetwork.levelInterfaceNetwork(0);
LevelInterfaceNetworkProblem<DGBasis, MatrixType, VectorType> coarseGlobalAssembler(coarseLevelInterfaceNetwork);
LocalOscAssembler coarseLocalAssembler(oscDataHierarchy[0]->data(), oscDataHierarchy[0]->mapper());
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);
const double pen = (1+1.0/c) * std::pow(1.0+c, k) * std::pow(2.0, k);
coarseLocalInterfaceAssemblers[i] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
}
coarseGlobalAssembler.assemble(coarseLocalAssembler, coarseLocalInterfaceAssemblers, l2FunctionalAssembler);
coarseGlobalAssembler.solve();
coarseSol = coarseGlobalAssembler.getSol();
coarseBasis = coarseGlobalAssembler.basis();
} else {
//TODO
const LevelInterfaceNetwork<GV>& nLevelInterfaceNetwork = interfaceNetwork.levelInterfaceNetwork(fineLevelIdx-1);
coarseBasis = std::make_shared<DGBasis>(nLevelInterfaceNetwork);
coarseSol.resize(coarseBasis->size());
std::stringstream solFilename;
solFilename << resultPath << "solutions/level_" << std::to_string(fineCantorLevel-1);
std::ifstream file(solFilename.str().c_str(), std::ios::in|std::ios::binary);
if (not(file))
DUNE_THROW(SolverError, "Couldn't open file " << solFilename.str() << " for reading");
Dune::MatrixVector::Generic::readBinary(file, coarseSol);
file.close();
}
// ----------------------------------------------------------------
// --- compute exact solution ---
// ----------------------------------------------------------------
std::cout << std::endl;
std::cout << "Computing exact solution!" << std::endl;
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 = minCantorLevel + 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);
}
exactGlobalAssembler.assemble(exactLocalAssembler, exactLocalInterfaceAssemblers, l2FunctionalAssembler);
exactGlobalAssembler.solve();
const VectorType& exactSol = exactGlobalAssembler.getSol();
std::shared_ptr<DGBasis> exactBasis = exactGlobalAssembler.basis();
EnergyNorm<MatrixType, VectorType> exactEnergyNorm(exactGlobalAssembler.matrix());
// ----------------------------------------------------------------
// --- set up cg solver ---
// ----------------------------------------------------------------
std::cout << std::endl;
std::cout << "---------------------------------------------" << std::endl;
std::cout << "Setting up the fine level CG solver:" << std::endl;
std::cout << "---------------------------------------------" << std::endl << std::endl;
const LevelInterfaceNetwork<GV>& fineLevelInterfaceNetwork = interfaceNetwork.levelInterfaceNetwork(fineLevelIdx);
LevelInterfaceNetworkProblem<DGBasis, MatrixType, VectorType> fineGlobalAssembler(fineLevelInterfaceNetwork);
LocalOscAssembler fineLocalAssembler(oscDataHierarchy[fineLevelIdx]->data(), oscDataHierarchy[fineLevelIdx]->mapper());
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);
const double pen = (1+1.0/c) * std::pow(1.0+c, k) * std::pow(2.0, k);
fineLocalInterfaceAssemblers[i] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
}
fineGlobalAssembler.assemble(fineLocalAssembler, fineLocalInterfaceAssemblers, l2FunctionalAssembler);
fineGlobalAssembler.solve();
const VectorType& fineSol = fineGlobalAssembler.getSol();
const VectorType& fineRhs = fineGlobalAssembler.rhs();
std::shared_ptr<DGBasis> fineBasis = fineGlobalAssembler.basis();
EnergyNorm<MatrixType, VectorType> fineEnergyNorm(fineGlobalAssembler.matrix());
std::cout << "Setting up the preconditioner!" << std::endl;
BitVector activeLevels(fineLevelIdx+1, false);
for (size_t i=0; i<activeLevels.size(); i++, i++) {
activeLevels[i][0] = true;
}
std::vector<std::shared_ptr<LocalOscAssembler>> localAssemblers;
std::vector<std::vector<std::shared_ptr<LocalInterfaceAssembler>>> localIntersectionAssemblers;
localAssemblers.resize(activeLevels.size());
localIntersectionAssemblers.resize(activeLevels.size());
for (size_t i=0; i<activeLevels.size(); i++) {
localAssemblers[i] = std::make_shared<LocalOscAssembler>(oscDataHierarchy[i]->data(), oscDataHierarchy[i]->mapper());
const LevelInterfaceNetwork<GV>& levelInterfaceNetwork = interfaceNetwork.levelInterfaceNetwork(i);
std::vector<std::shared_ptr<LocalInterfaceAssembler>>& levelLocalIntersectionAssemblers = localIntersectionAssemblers[i];
levelLocalIntersectionAssemblers.resize(levelInterfaceNetwork.size());
for (size_t j=0; j<levelLocalIntersectionAssemblers.size(); j++) {
const int k = minCantorLevel + levelInterfaceNetwork.getIntersectionsLevels().at(j);
const double pen = (1+1.0/c) * std::pow(1.0+c, k) * std::pow(2.0, k);
levelLocalIntersectionAssemblers[j] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
}
}
const auto& preconditionerParset = parameterSet.sub("preconditioner");
using MultilevelPatchPreconditioner = MultilevelPatchPreconditioner<DGBasis, LocalOscAssembler, LocalInterfaceAssembler, MatrixType, VectorType>;
MultilevelPatchPreconditioner preconditioner(interfaceNetwork, activeLevels, localAssemblers, localIntersectionAssemblers, preconditionerParset);
/*for (size_t i=0; i<preconditioner.size()-1; i++) {
auto& levelFaultPreconditioner = *preconditioner.levelPatchPreconditioner(i);
levelFaultPreconditioner.setPatchDepth(patchDepth);
levelFaultPreconditioner.setBoundaryMode(LevelPatchPreconditioner::BoundaryMode::homogeneous);
}*/
preconditioner.build();
std::cout << "Setup complete, starting preconditioned cg iteration!" << std::endl;
std::cout << std::endl << std::endl;
// set initial iterate
VectorType initialX;
DGMGTransfer<DGBasis> 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)
//solver.historyBuffer_ = resultPath + "solutions/level_" + std::to_string(fineCantorLevel) + ".vec";
solver.check();
solver.preprocess();
solver.solve();
std::cout << std::endl;
std::cout << "Total time: " << totalTimer.elapsed() << " seconds." << std::endl;
std::cout.rdbuf(coutbuf); //reset to standard output again
std::cout << "Problem " << problemCount << " done!" << std::endl;
problemCount++;
}
} catch (Dune::Exception e) {
std::cout << e << std::endl;
}
}