Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
#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
*/
#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/bcrsmatrix.hh>
#include <dune/istl/matrix.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/assemblers/interfaceoperatorassembler.hh>
#include <dune/faultnetworks/assemblers/intersectionjumpmassassembler.hh>
#include <dune/faultnetworks/assemblers/osclocalassembler.hh>
#include <dune/faultnetworks/solvers/osccgsolver.hh>
#include <dune/faultnetworks/preconditioners/multilevelpatchpreconditioner.hh>
#include <dune/faultnetworks/faultp1nodalbasis.hh>
#include <dune/faultnetworks/levelinterfacenetwork.hh>
#include <dune/faultnetworks/faultinterface.hh>
#include <dune/faultnetworks/levelinterfacenetworkproblem.hh>
#include <dune/faultnetworks/compressedmultigridtransfer.hh>
#include <dune/faultnetworks/dgmgtransfer.hh>
#include <dune/faultnetworks/oscdata.hh>
#include <dune/faultnetworks/oscdatahierarchy.hh>
#include <dune/faultnetworks/oscrhs.hh>
#include <dune/faultnetworks/faultfactories/geofaultfactory.hh>
int main(int argc, char** argv) { try
{
MPIHelper::instance(argc, argv);
const int dim = 2;
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
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 int coarseResolution = problemParameters.get<int>("coarseResolution");
const int fineResolution = problemParameters.get<int>("fineResolution");
const int exactResolution = problemParameters.get<int>("exactResolution");
const int faultCount = problemParameters.get<int>("faultCount");
const double c = problemParameters.get<double>("penaltyFactor");
const int patchDepth = problemParameters.get<int>("patchDepth");
const int maxIterations = problemParameters.get<int>("maxIterations");
const double solverTolerance = problemParameters.get<double>("tol");
const bool nestedIteration = problemParameters.get<bool>("nestedIteration");
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
const int fineGridN = std::pow(2, fineResolution);
std::ofstream out(resultPath + oscDataFile + "_" + std::to_string(coarseResolution) + "_" + std::to_string(fineResolution) + "_" + std::to_string(exactResolution) + "_" + std::to_string(faultCount) + ".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)) - coarseResolution;
if (oscGridN>fineGridN)
DUNE_THROW(Dune::Exception, "Provided oscData too large!");
std::cout << "OscData file read successfully!" << std::endl << std::endl;
const int exactLevelIdx = exactResolution - coarseResolution;
int fineLevelIdx = exactLevelIdx - 1;
/*
// exponential distribution of faults onto levels
std::vector<int> faultToLevel(1);
faultToLevel[0] = 1;
//faultCountdown = faultCount-1;
for (int i=2; i<=fineLevelIdx; i++) {
std::vector<int> help(faultToLevel.size() + (int) std::pow(2.0, i-1), i);
for (size_t j=0; j<faultToLevel.size(); j++) {
help[1+2*j] = faultToLevel[j];
}
faultToLevel = help;
}*/
// uniform distribution of faults onto levels
std::vector<int> faultToLevel(faultCount);
for (size_t i=0; i<faultToLevel.size(); ) {
for (int j=fineLevelIdx; j>0 && i<faultToLevel.size(); j--) {
faultToLevel[i] = j;
i++;
}
for (int j=0; j<fineLevelIdx && i<faultToLevel.size(); j++) {
faultToLevel[i] = j;
i++;
}
}
fineLevelIdx = fineResolution - coarseResolution;
print(faultToLevel, "faultToLevel: ");
GeoFaultFactory<GridType> faultFactory(faultToLevel, coarseResolution, exactLevelIdx, false);
const GridType& grid = faultFactory.grid();
const InterfaceNetwork<GridType>& interfaceNetwork = faultFactory.interfaceNetwork();
std::vector<int> writeLevels {0, 1, fineLevelIdx};
for (size_t i=0; i<writeLevels.size(); i++) {
int writeLevel = writeLevels[i];
bool writeGrid = !(writeLevel==fineLevelIdx);
LevelInterfaceNetworkWriter networkWriter(resultPath + "levelinterfacenetwork_" + std::to_string(writeLevel) + ".tikz");
networkWriter.write(interfaceNetwork.levelInterfaceNetwork(writeLevel), writeGrid);
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 = coarseResolution + coarseLevelInterfaceNetwork.getIntersectionsLevels().at(i);
const double pen = (1+1.0/c) * std::pow(1.0+c, k) * (std::pow(2.0, k) - 1.0);
coarseLocalInterfaceAssemblers[i] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
}
coarseGlobalAssembler.assemble(coarseLocalAssembler, coarseLocalInterfaceAssemblers, l2FunctionalAssembler);
coarseGlobalAssembler.solve();
coarseSol = coarseGlobalAssembler.getSol();
coarseBasis = coarseGlobalAssembler.basis();
} else {
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(fineResolution-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 = coarseResolution + exactLevelInterfaceNetwork.getIntersectionsLevels().at(i);
const double pen = (1+1.0/c) * std::pow(1.0+c, k) * (std::pow(2.0, k) - 1.0);
exactLocalInterfaceAssemblers[i] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
}
exactGlobalAssembler.assemble(exactLocalAssembler, exactLocalInterfaceAssemblers, l2FunctionalAssembler);
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 = coarseResolution + fineLevelInterfaceNetwork.getIntersectionsLevels().at(i);
const double pen = (1+1.0/c) * std::pow(1.0+c, k) * (std::pow(2.0, k) - 1.0);
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, 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 = coarseResolution + levelInterfaceNetwork.getIntersectionsLevels().at(j);
const double pen = (1+1.0/c) * std::pow(1.0+c, k) * (std::pow(2.0, k) - 1.0);
levelLocalIntersectionAssemblers[j] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
}
}
typedef MultilevelPatchPreconditioner<DGBasis, LocalOscAssembler, LocalInterfaceAssembler, MatrixType, VectorType> MultilevelPatchPreconditioner;
typedef typename MultilevelPatchPreconditioner::LevelPatchPreconditionerType LevelPatchPreconditioner;
MultilevelPatchPreconditioner::Mode preconditionerMode = MultilevelPatchPreconditioner::Mode::additive;
MultilevelPatchPreconditioner preconditioner(interfaceNetwork, activeLevels, localAssemblers, localIntersectionAssemblers, preconditionerMode);
for (size_t i=0; i<preconditioner.size()-1; i++) {
LevelPatchPreconditioner& 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;
DGMGTransfer<DGBasis> coarseFineTransfer(*coarseBasis, *fineBasis);
coarseFineTransfer.prolong(coarseSol, initialX);
VectorType rhs(fineRhs);
DGMGTransfer<DGBasis> discMGTransfer(*fineBasis, *exactBasis);
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(fineResolution);
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;
}
}