diff --git a/dune/faultnetworks/levelinterfacenetworkproblem.hh b/dune/faultnetworks/levelinterfacenetworkproblem.hh
index 0ff14d8438729b0826b8d7d6333c4d9a3ebe5628..ffa55382f9082822e43d5416d89ea58411771a5d 100644
--- a/dune/faultnetworks/levelinterfacenetworkproblem.hh
+++ b/dune/faultnetworks/levelinterfacenetworkproblem.hh
@@ -10,7 +10,6 @@
#include <dune/faultnetworks/faultinterface.hh>
#include <dune/faultnetworks/oscdata.hh>
-//#include <dune/istl/superlu.hh>
#include <dune/istl/umfpack.hh>
#include <dune/faultnetworks/utils/debugutils.hh>
@@ -92,14 +91,14 @@ class LevelInterfaceNetworkProblem
timer.reset();
timer.start();
- #if HAVE_SUPERLU
+ #if HAVE_UMFPACK
Dune::InverseOperatorResult res;
VectorType rhsCopy(rhs_);
Dune::UMFPack<MatrixType> solver(matrix_);
solver.apply(x_, rhsCopy, res);
#else
- #error No SuperLU!
+ #error No UMFPack!
#endif
std::cout << "Solving global problem took: " << timer.elapsed() << " seconds" << std::endl;
diff --git a/dune/faultnetworks/localproblem.hh b/dune/faultnetworks/localproblem.hh
index 9e409f4b239a1eefe7171de52ed10ade9494ede7..3da70633f19152350c2a445195a2ad1881f75f2d 100644
--- a/dune/faultnetworks/localproblem.hh
+++ b/dune/faultnetworks/localproblem.hh
@@ -7,7 +7,7 @@
#include <dune/common/timer.hh>
#include <dune/istl/matrixindexset.hh>
-//#include <dune/istl/superlu.hh>
+
#include <dune/istl/umfpack.hh>
#include <dune/fufem/assemblers/localoperatorassembler.hh>
@@ -146,7 +146,7 @@ public:
}
void solve(DomainType& x){
- #if HAVE_SUPERLU
+ #if HAVE_UMFPACK
RangeType localRhsCopy(localRhs);
Dune::InverseOperatorResult res;
@@ -155,7 +155,7 @@ public:
Dune::UMFPack<MatrixType> directSolver(localMat);
directSolver.apply(x, localRhsCopy, res);
#else
- #error No SuperLU!
+ #error No UMFPack!
#endif
}
diff --git a/dune/faultnetworks/solvers/osccgsolver.cc b/dune/faultnetworks/solvers/osccgsolver.cc
index 67527edf233c5787f44129a04092a174095c6586..25dac2d5f575268ebdbe838e65dcd8da2c0caced 100644
--- a/dune/faultnetworks/solvers/osccgsolver.cc
+++ b/dune/faultnetworks/solvers/osccgsolver.cc
@@ -162,8 +162,6 @@ void OscCGSolver<MatrixType, VectorType, DGBasisType>::solve()
VectorType iterationError(*x_);
iterationError -= *fineSol_;
- writeToVTK(fineBasis_, iterationError, "/storage/mi/podlesny/data/faultnetworks/iterationError/multilevel", "exactvertexdata_step_"+std::to_string(i));
-
double normOfIterationError = fineErrorNorm.operator()(iterationError);
double convRate = normOfIterationError / normOfOldIterationError;
@@ -199,14 +197,14 @@ void OscCGSolver<MatrixType, VectorType, DGBasisType>::solve()
std::cout << std::setiosflags(std::ios::scientific);
std::cout << std::setw(15) << std::setprecision(7) << history;
std::cout << std::resetiosflags(std::ios::scientific);
-/*
- std::cout << std::setiosflags(std::ios::scientific);
+
+ std::cout << std::setiosflags(std::ios::scientific);
std::cout << std::setw(15) << std::setprecision(7) << normOfIterationError;
std::cout << std::resetiosflags(std::ios::scientific);
- std::cout << std::setiosflags(std::ios::scientific);
+ std::cout << std::setiosflags(std::ios::scientific);
std::cout << std::setw(15) << std::setprecision(7) << (terminationFactor_*normOfFinalDiscretizationError);
- std::cout << std::resetiosflags(std::ios::scientific); */
+ std::cout << std::resetiosflags(std::ios::scientific);
std::cout << std::endl;
}
diff --git a/dune/faultnetworks/utils/levelinterfacenetworkwriter.hh b/dune/faultnetworks/utils/levelinterfacenetworkwriter.hh
index 0c76116b925222bb8ea07aba99a88b3459fd422d..b4f43e15c5417d324c286ba05535f9c2858b24b6 100644
--- a/dune/faultnetworks/utils/levelinterfacenetworkwriter.hh
+++ b/dune/faultnetworks/utils/levelinterfacenetworkwriter.hh
@@ -107,7 +107,7 @@ class LevelInterfaceNetworkWriter {
const std::vector<int>& intersectionsLevels = levelInterfaceNetwork.getIntersectionsLevels();
for (size_t i=0; i<intersections.size(); i++) {
if (highlightNewLevelInterfaces) {
- writeIntersection(intersections[i], file, 2*(intersectionsLevels[i]==levelInterfaceNetwork.level()));
+ writeIntersection(intersections[i], file, 3*(intersectionsLevels[i]==levelInterfaceNetwork.level()));
} else {
writeIntersection(intersections[i], file, intersectionsLevels[i]);
}
diff --git a/src/cantorfaultnetworks/cantorconvergence.cc b/src/cantorfaultnetworks/cantorconvergence.cc
index f72393746222628a5ff29b68ffb4b19bbaf8a8c6..271b9fbdedd418317da988b66a473c40c3a645e6 100644
--- a/src/cantorfaultnetworks/cantorconvergence.cc
+++ b/src/cantorfaultnetworks/cantorconvergence.cc
@@ -79,7 +79,6 @@ int main(int argc, char** argv) { try
#error No UG!
#endif
-
typedef typename GridType::LevelGridView GV;
typedef FaultP1NodalBasis<GV, double> DGBasis;
typedef typename DGBasis::LocalFiniteElement DGFE;
@@ -92,7 +91,7 @@ int main(int argc, char** argv) { try
if (argc==2)
ParameterTreeParser::readINITree(argv[1], parameterSet);
else
- ParameterTreeParser::readINITree("/home/mi/podlesny/software/dune/faultnetworks/src/cantorfaultnetworks/cantorconvergence.parset", parameterSet);
+ ParameterTreeParser::readINITree("cantorconvergence.parset", parameterSet);
const std::string path = parameterSet.get<std::string>("path");
const std::string resultPath = parameterSet.get<std::string>("resultPath");
@@ -154,12 +153,6 @@ int main(int argc, char** argv) { try
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();
@@ -196,7 +189,6 @@ int main(int argc, char** argv) { try
std::vector<std::shared_ptr<LocalInterfaceAssembler>> exactLocalInterfaceAssemblers(exactLevelInterfaceNetwork.size());
for (size_t j=0; j<exactLocalInterfaceAssemblers.size(); j++) {
- //const double pen = penaltyFactor;
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);
@@ -219,7 +211,6 @@ int main(int argc, char** argv) { try
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);
@@ -237,12 +228,10 @@ int main(int argc, char** argv) { try
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);
}
}
diff --git a/src/cantorfaultnetworks/cantorconvergence.parset b/src/cantorfaultnetworks/cantorconvergence.parset
index 1be006a4b6e9615f16c77637e7695bea024bbcff..1fc5ed4f5761c61adb46f9189b565410932ae13f 100644
--- a/src/cantorfaultnetworks/cantorconvergence.parset
+++ b/src/cantorfaultnetworks/cantorconvergence.parset
@@ -1,5 +1,5 @@
-path = /storage/mi/podlesny/data/osccoeff/
-resultPath = /home/mi/podlesny/results/faultnetworks/cantorfaultnetworks/cantorconvergence/
+path = ../data/
+resultPath = ../results/cantorfaultnetworks/cantorconvergence
[problem0]
diff --git a/src/cantorfaultnetworks/cantorfaultnetwork.cc b/src/cantorfaultnetworks/cantorfaultnetwork.cc
index a21885f3da57861c9b8d3ac0eb9baa730c8581de..695f5f9dd3deddda83a4a6d6ec7c98f18649fee6 100644
--- a/src/cantorfaultnetworks/cantorfaultnetwork.cc
+++ b/src/cantorfaultnetworks/cantorfaultnetwork.cc
@@ -93,7 +93,7 @@ int main(int argc, char** argv) { try
if (argc==2)
ParameterTreeParser::readINITree(argv[1], parameterSet);
else
- ParameterTreeParser::readINITree("/home/mi/podlesny/software/dune/faultnetworks/src/cantorfaultnetworks/cantorfaultnetwork.parset", parameterSet);
+ ParameterTreeParser::readINITree("cantorfaultnetwork.parset", parameterSet);
const std::string path = parameterSet.get<std::string>("path");
const std::string resultPath = parameterSet.get<std::string>("resultPath");
@@ -118,7 +118,6 @@ int main(int argc, char** argv) { try
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::streambuf *coutbuf = std::cout.rdbuf(); //save old buf
std::cout.rdbuf(out.rdbuf()); //redirect std::cout to out.txt
@@ -160,9 +159,11 @@ 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();
- for (size_t i=0; i<interfaceNetwork.size(); i++) {
- LevelInterfaceNetworkWriter networkWriter(resultPath + "levelinterfacenetwork_" + std::to_string(i) + ".tikz");
- networkWriter.write(interfaceNetwork.levelInterfaceNetwork(i));
+
+ if (problemCount==0) {
+ int writeLevel = 1;
+ LevelInterfaceNetworkWriter networkWriter(resultPath + "levelinterfacenetwork_" + std::to_string(writeLevel) + ".tikz");
+ networkWriter.write(interfaceNetwork.levelInterfaceNetwork(writeLevel), false);
}
const GridType& grid = faultFactory.grid();
@@ -191,8 +192,6 @@ 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 double pen = penaltyFactor;
- //const double pen = penaltyFactor*std::pow(2.0, coarseResolution);
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);
coarseLocalInterfaceAssemblers[i] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
@@ -216,8 +215,6 @@ 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 double pen = penaltyFactor;
- //const double pen = penaltyFactor*std::pow(2.0, coarseResolution + exactLevelInterfaceNetwork.getIntersectionsLevels().at(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);
exactLocalInterfaceAssemblers[i] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
@@ -231,7 +228,6 @@ int main(int argc, char** argv) { try
EnergyNorm<MatrixType, VectorType> exactEnergyNorm(exactGlobalAssembler.matrix());
-
// ----------------------------------------------------------------
// --- set up cg solver ---
// ----------------------------------------------------------------
@@ -247,8 +243,6 @@ 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 double pen = penaltyFactor;
- //const double pen = penaltyFactor*std::pow(2.0, coarseResolution + fineLevelInterfaceNetwork.getIntersectionsLevels().at(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);
fineLocalInterfaceAssemblers[i] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
@@ -263,18 +257,12 @@ int main(int argc, char** argv) { try
EnergyNorm<MatrixType, VectorType> fineEnergyNorm(fineGlobalAssembler.matrix());
- // ---------------------------------------------------------------
// multilevel transfer
- // ---------------------------------------------------------------
DGMGTransfer<DGBasis> coarseFineTransfer(coarseBasis, fineBasis);
std::cout << "Setting up the preconditioner!" << std::endl;
BitVector activeLevels(fineLevelIdx+1, true);
- //activeLevels[0][0] = true;
- //activeLevels[1][0] = false;
- //activeLevels[2][0] = false;
- //activeLevels[fineLevelIdx][0] = true;
std::vector<std::shared_ptr<LocalOscAssembler>> localAssemblers;
std::vector<std::vector<std::shared_ptr<LocalInterfaceAssembler>>> localIntersectionAssemblers;
@@ -289,8 +277,6 @@ int main(int argc, char** argv) { try
std::vector<std::shared_ptr<LocalInterfaceAssembler>>& levelLocalIntersectionAssemblers = localIntersectionAssemblers[i];
levelLocalIntersectionAssemblers.resize(levelInterfaceNetwork.size());
for (size_t j=0; j<levelLocalIntersectionAssemblers.size(); j++) {
- //const double pen = penaltyFactor;
- //const double pen = penaltyFactor*std::pow(2.0, coarseResolution + levelInterfaceNetwork.getIntersectionsLevels().at(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);
levelLocalIntersectionAssemblers[j] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
@@ -299,7 +285,6 @@ int main(int argc, char** argv) { try
typedef MultilevelPatchPreconditioner<DGBasis, LocalOscAssembler, LocalInterfaceAssembler, MatrixType, VectorType> MultilevelPatchPreconditioner;
typedef typename MultilevelPatchPreconditioner::LevelPatchPreconditionerType LevelPatchPreconditioner;
- //typedef typename MultilevelPatchPreconditioner::LevelGlobalPreconditionerType LevelGlobalPreconditioner;
MultilevelPatchPreconditioner::Mode preconditionerMode = MultilevelPatchPreconditioner::Mode::additive;
MultilevelPatchPreconditioner preconditioner(interfaceNetwork, activeLevels, localAssemblers, localIntersectionAssemblers, preconditionerMode);
@@ -313,33 +298,12 @@ int main(int argc, char** argv) { try
std::cout << "Setup complete, starting preconditioned cg iteration!" << std::endl;
std::cout << std::endl << std::endl;
-
- /*Dune::BlockVector<Dune::FieldVector<double,1>> initialError(fineX.size());
- initialError = initialX;
- initialError -= fineX;
- std::cout << "Initial error in energy norm: " << energyNorm.operator()(initialError) << std::endl;
- std::cout << "Initial error in two norm: " << initialError.two_norm() << std::endl << std::endl;*/
-
VectorType initialX;
coarseFineTransfer.prolong(coarseSol, initialX);
- VectorType initialXCopy(initialX);
- VectorType fineRhsCopy(fineRhs);
VectorType rhs(fineRhs);
- //print(initialX, "initialX");
- //print(fineRhs, "fineRhs");
- //print(exactSol, "exactSol");
- //print(fineSol, "fineSol");
-
- /* Basis fineContBasis(fineGridView);
- std::cout << "fineContBasis.size(): " << fineContBasis.size() << std::endl;
- DGMGTransfer<Basis> contMGTransfer(coarseBasis, fineContBasis);
-
- Basis exactContBasis(exactGridView);
- std::cout << "exactContBasis.size(): " << exactContBasis.size() << std::endl; */
-
DGMGTransfer<DGBasis> discMGTransfer(fineBasis, exactBasis);
// solve
@@ -351,31 +315,9 @@ int main(int argc, char** argv) { try
solver.preprocess();
solver.solve();
-
- /*
- VectorType finalError(fineSol.size());
- finalError = initialX;
- finalError -= fineSol;
- std::cout << "Final error in energy norm: " << fineEnergyNorm.operator()(finalError) << std::endl;
- std::cout << "Final error in two norm: " << finalError.two_norm() << std::endl << std::endl;*/
-
std::cout << std::endl;
std::cout << "Total time: " << totalTimer.elapsed() << " seconds." << std::endl;
- /*
- // regular unpreconditioned cg solver for comparisson
- std::cout << "Setup complete, starting cg iteration!" << std::endl;
- std::cout << std::endl;
-
- // solve
- OscCGSolver<MatrixType, VectorType, DGBasis >
- regSolver(fineBasis, &fineGlobalAssembler.matrix(), &initialXCopy, &fineRhsCopy, nullptr,
- maxIterations, solverTolerance, &exactEnergyNorm,
- Solver::FULL, &exactSol, &fineSol, discMGTransfer, 1.0, true);
- regSolver.check();
- regSolver.preprocess();
- regSolver.solve();*/
-
std::cout.rdbuf(coutbuf); //reset to standard output again
std::cout << "Problem " << problemCount << " done!" << std::endl;
diff --git a/src/cantorfaultnetworks/cantorfaultnetwork.parset b/src/cantorfaultnetworks/cantorfaultnetwork.parset
index 4215d68442f32e2ab89ad697fcc182faa6a9a31e..6aa5bb6e122b40c435892505464dd9ce4582568a 100644
--- a/src/cantorfaultnetworks/cantorfaultnetwork.parset
+++ b/src/cantorfaultnetworks/cantorfaultnetwork.parset
@@ -1,5 +1,5 @@
-path = /storage/mi/podlesny/data/osccoeff/
-resultPath = /home/mi/podlesny/results/faultnetworks/cantorfaultnetworks/cantorfaultnetwork/
+path = ../data/
+resultPath = ../results/cantorfaultnetworks/
@@ -10,7 +10,7 @@ oscDataFile = oscDataLaplace32.mat
coarseResolution = 4
fineResolution = 5
exactResolution = 10
-minCantorResolution = 1
+minCantorResolution = 4
penaltyFactor = 1
patchDepth = 1
@@ -29,7 +29,7 @@ oscDataFile = oscDataLaplace32.mat
coarseResolution = 4
fineResolution = 6
exactResolution = 10
-minCantorResolution = 1
+minCantorResolution = 4
penaltyFactor = 1
patchDepth = 1
diff --git a/src/data/oscDataLaplace32.mat b/src/data/oscDataLaplace32.mat
new file mode 100644
index 0000000000000000000000000000000000000000..e422df8bbe17487fb5afcf7c258f6aa081772d3d
--- /dev/null
+++ b/src/data/oscDataLaplace32.mat
@@ -0,0 +1,33 @@
+32 32
+1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
diff --git a/src/geofaultnetworks/geofaultnetwork.cc b/src/geofaultnetworks/geofaultnetwork.cc
index 1572c52fcdad3990703c8a48bd17748e9b275a82..1da03affcfa8491b56bb9ec001924f408dfbb091 100644
--- a/src/geofaultnetworks/geofaultnetwork.cc
+++ b/src/geofaultnetworks/geofaultnetwork.cc
@@ -69,8 +69,6 @@ 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;
@@ -98,7 +96,7 @@ int main(int argc, char** argv) { try
if (argc==2)
ParameterTreeParser::readINITree(argv[1], parameterSet);
else
- ParameterTreeParser::readINITree("/home/mi/podlesny/software/dune/faultnetworks/src/geofaultnetworks/geofaultnetwork.parset", parameterSet);
+ ParameterTreeParser::readINITree("geofaultnetwork.parset", parameterSet);
const std::string path = parameterSet.get<std::string>("path");
const std::string resultPath = parameterSet.get<std::string>("resultPath");
@@ -121,12 +119,9 @@ int main(int argc, char** argv) { try
const int maxIterations = problemParameters.get<int>("maxIterations");
const double solverTolerance = problemParameters.get<double>("tol");
- //const bool useExactSol = problemParameters.get<bool>("useExactSol");
- //const int coarseGridN = std::pow(2, coarseResolution);
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!
@@ -157,8 +152,6 @@ int main(int argc, char** argv) { try
std::cout << "OscData file read successfully!" << std::endl << std::endl;
const int exactLevelIdx = exactResolution - coarseResolution;
- // TODO: reverse after HomTap computations
- //const int fineLevelIdx = fineResolution - coarseResolution;
int fineLevelIdx = exactLevelIdx - 1;
/*
@@ -189,7 +182,6 @@ int main(int argc, char** argv) { try
}
}
- // das auch
fineLevelIdx = fineResolution - coarseResolution;
print(faultToLevel, "faultToLevel: ");
@@ -199,15 +191,12 @@ int main(int argc, char** argv) { try
const GridType& grid = faultFactory.grid();
const InterfaceNetwork<GridType>& interfaceNetwork = faultFactory.interfaceNetwork();
- for (size_t i=0; i<interfaceNetwork.size(); i++) {
- LevelInterfaceNetworkWriter networkWriter(resultPath + "levelinterfacenetwork_" + std::to_string(i) + ".tikz");
- networkWriter.write(interfaceNetwork.levelInterfaceNetwork(i));
+ if (problemCount==0) {
+ int writeLevel = fineLevelIdx;
+ LevelInterfaceNetworkWriter networkWriter(resultPath + "levelinterfacenetwork_" + std::to_string(writeLevel) + ".tikz");
+ networkWriter.write(interfaceNetwork.levelInterfaceNetwork(writeLevel), false);
}
- /*int writeLevel = fineLevelIdx;
- LevelInterfaceNetworkWriter networkWriter(resultPath + "levelinterfacenetwork_" + std::to_string(writeLevel) + ".tikz");
- networkWriter.write(interfaceNetwork.levelInterfaceNetwork(writeLevel), true);*/
-
OscData<GridType, OscMatrixType> oscData(grid, oscGridLevelIdx);
oscData.set(matrix);
OscDataHierarchy<GridType, OscMatrixType> oscDataHierarchy(oscData);
@@ -232,8 +221,6 @@ 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 double pen = penaltyFactor;
- //const double pen = penaltyFactor*std::pow(2.0, coarseResolution);
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);
coarseLocalInterfaceAssemblers[i] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
@@ -257,21 +244,16 @@ 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 double pen = penaltyFactor;
- //const double pen = penaltyFactor*std::pow(2.0, coarseResolution + exactLevelInterfaceNetwork.getIntersectionsLevels().at(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);
exactLocalInterfaceAssemblers[i] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
}
exactGlobalAssembler.assemble(exactLocalAssembler, exactLocalInterfaceAssemblers, l2FunctionalAssembler);
- // TODO: turn on exact sol? needed for computation of discretization error
- //exactGlobalAssembler.solve();
-
- //const VectorType& exactSol = exactGlobalAssembler.getSol();
+ exactGlobalAssembler.solve();
+ const VectorType& exactSol = exactGlobalAssembler.getSol();
const DGBasis& exactBasis = exactGlobalAssembler.basis();
- VectorType exactSol(exactBasis.size(), 0);
EnergyNorm<MatrixType, VectorType> exactEnergyNorm(exactGlobalAssembler.matrix());
@@ -290,8 +272,6 @@ 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 double pen = penaltyFactor;
- //const double pen = penaltyFactor*std::pow(2.0, coarseResolution + fineLevelInterfaceNetwork.getIntersectionsLevels().at(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);
fineLocalInterfaceAssemblers[i] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
@@ -306,18 +286,12 @@ int main(int argc, char** argv) { try
EnergyNorm<MatrixType, VectorType> fineEnergyNorm(fineGlobalAssembler.matrix());
- // ---------------------------------------------------------------
// multilevel transfer
- // ---------------------------------------------------------------
DGMGTransfer<DGBasis> coarseFineTransfer(coarseBasis, fineBasis);
std::cout << "Setting up the preconditioner!" << std::endl;
BitVector activeLevels(fineLevelIdx+1, true);
- //activeLevels[0][0] = true;
- //activeLevels[1][0] = false;
- //activeLevels[2][0] = false;
- //activeLevels[fineLevelIdx][0] = true;
std::vector<std::shared_ptr<LocalOscAssembler>> localAssemblers;
std::vector<std::vector<std::shared_ptr<LocalInterfaceAssembler>>> localIntersectionAssemblers;
@@ -332,8 +306,6 @@ int main(int argc, char** argv) { try
std::vector<std::shared_ptr<LocalInterfaceAssembler>>& levelLocalIntersectionAssemblers = localIntersectionAssemblers[i];
levelLocalIntersectionAssemblers.resize(levelInterfaceNetwork.size());
for (size_t j=0; j<levelLocalIntersectionAssemblers.size(); j++) {
- //const double pen = penaltyFactor;
- //const double pen = penaltyFactor*std::pow(2.0, coarseResolution + levelInterfaceNetwork.getIntersectionsLevels().at(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);
levelLocalIntersectionAssemblers[j] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
@@ -342,7 +314,6 @@ int main(int argc, char** argv) { try
typedef MultilevelPatchPreconditioner<DGBasis, LocalOscAssembler, LocalInterfaceAssembler, MatrixType, VectorType> MultilevelPatchPreconditioner;
typedef typename MultilevelPatchPreconditioner::LevelPatchPreconditionerType LevelPatchPreconditioner;
- //typedef typename MultilevelPatchPreconditioner::LevelGlobalPreconditionerType LevelGlobalPreconditioner;
MultilevelPatchPreconditioner::Mode preconditionerMode = MultilevelPatchPreconditioner::Mode::additive;
MultilevelPatchPreconditioner preconditioner(interfaceNetwork, activeLevels, localAssemblers, localIntersectionAssemblers, preconditionerMode);
@@ -357,32 +328,11 @@ int main(int argc, char** argv) { try
std::cout << "Setup complete, starting preconditioned cg iteration!" << std::endl;
std::cout << std::endl << std::endl;
- /*Dune::BlockVector<Dune::FieldVector<double,1>> initialError(fineX.size());
- initialError = initialX;
- initialError -= fineX;
- std::cout << "Initial error in energy norm: " << energyNorm.operator()(initialError) << std::endl;
- std::cout << "Initial error in two norm: " << initialError.two_norm() << std::endl << std::endl;*/
-
-
VectorType initialX;
coarseFineTransfer.prolong(coarseSol, initialX);
- VectorType initialXCopy(initialX);
- VectorType fineRhsCopy(fineRhs);
VectorType rhs(fineRhs);
- //print(initialX, "initialX");
- //print(fineRhs, "fineRhs");
- //print(exactSol, "exactSol");
- //print(fineSol, "fineSol");
-
- /* Basis fineContBasis(fineGridView);
- std::cout << "fineContBasis.size(): " << fineContBasis.size() << std::endl;
- DGMGTransfer<Basis> contMGTransfer(coarseBasis, fineContBasis);
-
- Basis exactContBasis(exactGridView);
- std::cout << "exactContBasis.size(): " << exactContBasis.size() << std::endl; */
-
DGMGTransfer<DGBasis> discMGTransfer(fineBasis, exactBasis);
// solve
@@ -394,32 +344,9 @@ int main(int argc, char** argv) { try
solver.preprocess();
solver.solve();
-
- /*
- VectorType finalError(fineSol.size());
- finalError = initialX;
- finalError -= fineSol;
- std::cout << "Final error in energy norm: " << fineEnergyNorm.operator()(finalError) << std::endl;
- std::cout << "Final error in two norm: " << finalError.two_norm() << std::endl << std::endl;*/
-
std::cout << std::endl;
std::cout << "Total time: " << totalTimer.elapsed() << " seconds." << std::endl;
- // regular unpreconditioned cg solver for comparisson
-/*
- std::cout << "Setup complete, starting cg iteration!" << std::endl;
- std::cout << std::endl;
-
- // solve
- OscCGSolver<MatrixType, VectorType, DGBasis >
- regSolver(fineBasis, &fineGlobalAssembler.matrix(), &initialXCopy, &fineRhsCopy, nullptr,
- maxIterations, solverTolerance, &exactEnergyNorm,
- Solver::FULL, &exactSol, &fineSol, discMGTransfer, 1.0, true);
- regSolver.check();
- regSolver.preprocess();
- regSolver.solve();
-*/
- // TODO: use once debugging finished
std::cout.rdbuf(coutbuf); //reset to standard output again
std::cout << "Problem " << problemCount << " done!" << std::endl;
diff --git a/src/geofaultnetworks/geofaultnetwork.parset b/src/geofaultnetworks/geofaultnetwork.parset
index 5deed292facda5f71d2c06a4141d265e6921b4f0..3bb72548ad3d5ef3c6a5254a82caa0e458600e59 100644
--- a/src/geofaultnetworks/geofaultnetwork.parset
+++ b/src/geofaultnetworks/geofaultnetwork.parset
@@ -1,5 +1,5 @@
-path = /storage/mi/podlesny/data/osccoeff/
-resultPath = /home/mi/podlesny/results/faultnetworks/geofaultnetworks/testing/fractalhomogenization/
+path = ../data/
+resultPath = ../results/geofaultnetworks/
[problem0]
@@ -17,6 +17,77 @@ patchDepth = 1
# cg solver
tol = 1e-12
maxIterations = 8
-useExactSol = 0
+
+###########################################
+
+[problem1]
+oscDataFile = oscDataLaplace32.mat
+
+# level resolution in 2^(-...)
+coarseResolution = 4
+fineResolution = 6
+exactResolution = 10
+faultCount = 30
+
+penaltyFactor = 1
+patchDepth = 1
+
+# cg solver
+tol = 1e-12
+maxIterations = 8
+
+###########################################
+
+[problem2]
+oscDataFile = oscDataLaplace32.mat
+
+# level resolution in 2^(-...)
+coarseResolution = 4
+fineResolution = 7
+exactResolution = 10
+faultCount = 30
+
+penaltyFactor = 1
+patchDepth = 1
+
+# cg solver
+tol = 1e-12
+maxIterations = 8
+
+###########################################
+
+[problem3]
+oscDataFile = oscDataLaplace32.mat
+
+# level resolution in 2^(-...)
+coarseResolution = 4
+fineResolution = 8
+exactResolution = 10
+faultCount = 30
+
+penaltyFactor = 1
+patchDepth = 1
+
+# cg solver
+tol = 1e-12
+maxIterations = 8
+
+###########################################
+
+[problem4]
+oscDataFile = oscDataLaplace32.mat
+
+# level resolution in 2^(-...)
+coarseResolution = 4
+fineResolution = 9
+exactResolution = 10
+faultCount = 30
+
+penaltyFactor = 1
+patchDepth = 1
+
+# cg solver
+tol = 1e-12
+maxIterations = 8
###########################################
diff --git a/src/results/cantorfaultnetworks/plot.tex b/src/results/cantorfaultnetworks/plot.tex
new file mode 100644
index 0000000000000000000000000000000000000000..5f0576d36dfcf10200cfa18a40fe4f7450244d54
--- /dev/null
+++ b/src/results/cantorfaultnetworks/plot.tex
@@ -0,0 +1,9 @@
+\documentclass[tikz]{standalone}
+
+\usepackage{tikz}
+
+\begin{document}
+
+\def\scale{3}
+\input{levelinterfacenetwork_1.tikz}
+\end{document}
\ No newline at end of file
diff --git a/src/results/geofaultnetworks/plot.tex b/src/results/geofaultnetworks/plot.tex
new file mode 100644
index 0000000000000000000000000000000000000000..b9761d1b2e3028c434b385d429265eb0356586b2
--- /dev/null
+++ b/src/results/geofaultnetworks/plot.tex
@@ -0,0 +1,9 @@
+\documentclass[tikz]{standalone}
+
+\usepackage{tikz}
+
+\begin{document}
+
+\def\scale{3}
+\input{levelinterfacenetwork_5.tikz}
+\end{document}
\ No newline at end of file