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
#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/cantorconvergencefaultfactory.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::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("/home/mi/podlesny/software/dune/faultnetworks/src/cantorfaultnetworks/cantorconvergence.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 int minResolution = problemParameters.get<int>("minResolution");
const int maxResolution = problemParameters.get<int>("maxResolution");
const int minCantorResolution = problemParameters.get<int>("minCantorResolution");
std::ofstream out(resultPath + oscDataFile + "_" + std::to_string(minResolution) + "_" + std::to_string(maxResolution) + "_" + std::to_string(exactResolution) + "_" + std::to_string(minCantorResolution) + "_" + std::to_string((int) c) + ".log");
118
119
120
121
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
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)) - minResolution;
if (oscGridN>maxGridN)
DUNE_THROW(Dune::Exception, "Provided oscData too large!");
std::cout << "OscData file read successfully!" << std::endl << std::endl;
const int maxLevelIdx = maxResolution - minResolution;
const int exactLevelIdx = exactResolution - minResolution;
std::vector<size_t> maxCantorLevels(exactLevelIdx+1);
for (int i=0; i<=exactLevelIdx; i++) {
maxCantorLevels[i] = i+minCantorResolution;
CantorConvergenceFaultFactory<GridType> faultFactory(minResolution, exactLevelIdx, minCantorResolution, maxCantorLevels);
const std::vector<std::shared_ptr<InterfaceNetwork<GridType>>>& interfaceNetworks = faultFactory.interfaceNetworks();
/*
for (size_t i=0; i<interfaceNetworks[1]->size(); i++) {
LevelInterfaceNetworkWriter networkWriter(resultPath + "levelinterfacenetwork_" + std::to_string(i) + ".tikz");
networkWriter.write(interfaceNetworks[1]->levelInterfaceNetwork(i), true);
}*/
std::cout << "Interface networks built successfully!" << std::endl << std::endl;
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 direct solutions ---
// ----------------------------------------------------------------
std::cout << std::endl;
std::cout << "Computing direct solutions!" << std::endl;
std::vector<VectorType> solutions;
solutions.resize(exactLevelIdx+1);
std::vector<double> approxErrors;
approxErrors.resize(maxLevelIdx+1);
double discError = NAN;
OscRhs<WorldVectorType, Dune::FieldVector<double,1>> f;
L2FunctionalAssembler<GridType, DGFE> l2FunctionalAssembler(f);
const LevelInterfaceNetwork<GV>& exactLevelInterfaceNetwork = interfaceNetworks.back()->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 j=0; j<exactLocalInterfaceAssemblers.size(); j++) {
const int k = minCantorResolution + exactLevelInterfaceNetwork.getIntersectionsLevels().at(j);
const double pen = (1+1.0/c) * std::pow(1.0+c, k) * std::pow(2.0, k-1);
exactLocalInterfaceAssemblers[j] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
exactGlobalAssembler.assemble(exactLocalAssembler, exactLocalInterfaceAssemblers, l2FunctionalAssembler);
exactGlobalAssembler.solve();
solutions[exactLevelIdx] = exactGlobalAssembler.getSol();
approxErrors[maxLevelIdx] = 0;
EnergyNorm<MatrixType, VectorType> energyNorm(exactGlobalAssembler.matrix());
for (int i=maxLevelIdx; i>=0; i--) {
const LevelInterfaceNetwork<GV>& levelInterfaceNetwork = interfaceNetworks[i]->levelInterfaceNetwork(maxLevelIdx);
LevelInterfaceNetworkProblem<DGBasis, MatrixType, VectorType> globalAssembler(levelInterfaceNetwork);
LocalOscAssembler localAssembler(oscDataHierarchy[maxLevelIdx]->data(), oscDataHierarchy[maxLevelIdx]->mapper());
std::vector<std::shared_ptr<LocalInterfaceAssembler>> localInterfaceAssemblers(levelInterfaceNetwork.size());
for (size_t j=0; j<localInterfaceAssemblers.size(); j++) {
//const double pen = penaltyFactor;
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);
localInterfaceAssemblers[j] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
}
globalAssembler.assemble(localAssembler, localInterfaceAssemblers, l2FunctionalAssembler);
globalAssembler.solve();
const VectorType& sol = globalAssembler.getSol();
const DGBasis& basis = globalAssembler.basis();
if (i==maxLevelIdx) {
VectorType discErrorV(solutions[exactLevelIdx]);
discErrorV -= solutions[maxLevelIdx];
//approxErrors[i] = approxError.two_norm();
discError= energyNorm.operator ()(discErrorV);
} else {
VectorType approxError(solutions[maxLevelIdx]);
approxError -= solutions[i];
//approxErrors[i] = approxError.two_norm();
approxErrors[i] = energyNorm.operator ()(approxError);
}
}
// print results
std::cout << std::endl << std::endl;
std::cout << "cantorLevel h approxError"<< std::endl;
std::cout << "------------------------------" << std::endl;
for (size_t i=0; i<approxErrors.size(); i++) {
std::cout << (minCantorResolution+i) << " 2^-" << maxResolution << " " << approxErrors[i] << std::endl;
std::cout << std::endl << std::endl;
std::cout << "exactResolution: 2^-" << exactResolution << std::endl;
std::cout << "discError: " << discError << 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) {