From 3cf78d553206f78eb60d271efb1133527207dbc0 Mon Sep 17 00:00:00 2001
From: podlesny <podlesny@zedat.fu-berlin.de>
Date: Tue, 2 Jun 2020 14:39:18 +0200
Subject: [PATCH] .

---
 .../multilevelpatchpreconditioner.hh          | 45 +++++++------------
 .../sparsecantorfaultnetwork.cc               | 15 ++++---
 .../sparsecantorfaultnetwork.parset           |  4 +-
 src/data/oscDataLaplace16.mat                 | 18 ++++++++
 src/data/oscDataLaplace8.mat                  |  9 ++++
 5 files changed, 53 insertions(+), 38 deletions(-)
 create mode 100644 src/data/oscDataLaplace16.mat
 create mode 100644 src/data/oscDataLaplace8.mat

diff --git a/dune/faultnetworks/preconditioners/multilevelpatchpreconditioner.hh b/dune/faultnetworks/preconditioners/multilevelpatchpreconditioner.hh
index 106edb1..d238009 100644
--- a/dune/faultnetworks/preconditioners/multilevelpatchpreconditioner.hh
+++ b/dune/faultnetworks/preconditioners/multilevelpatchpreconditioner.hh
@@ -178,8 +178,6 @@ private:
     }
 
     void iterateAdd() {
-        //*(this->x_) = 0;
-
         VectorType x;
 
         for (size_t i=0; i<levelPatchPreconditioners_.size(); i++) {
@@ -198,31 +196,6 @@ private:
 
         mgTransfer_[levelPatchPreconditioners_.size()]->prolong(it, x);
         *(this->x_) += x;
-
-
-        /*
-            levelPatchPreconditioners_[itCounter_]->iterate();
-            const VectorType& it = levelPatchPreconditioners_[itCounter_]->getSol();
-
-            mgTransfer_[itCounter_]->prolong(it, x);
-
-            //writeToVTK(*allFaultBasis_, x, "/home/mi/podlesny/data/faultnetworks/sol/", "preconditioner_step_"+std::to_string(i));
-
-            *(this->x_) += x;
-
-        coarseGlobalPreconditioner_->iterate();
-        const VectorType& it2 = coarseGlobalPreconditioner_->getSol();
-
-        mgTransfer_[levelPatchPreconditioners_.size()]->prolong(it2, x);
-        *(this->x_) += x;
-
-        itCounter_++;
-        itCounter_ = itCounter_ % (levelPatchPreconditioners_.size());
-
-        if (itCounter_==0)
-            std::cout << "Multilevel pass complete!" << std::endl;
-
-        */
     }
 
 
@@ -238,7 +211,13 @@ private:
         //print(levelRhs_[j], "localCoarseRhs: ");
         //writeToVTK(coarseGlobalPreconditioner_->basis(), levelRhs_[j], "/storage/mi/podlesny/data/faultnetworks/rhs/coarse", "exactvertexdata_step_"+std::to_string(itCounter_));
 
-        coarseGlobalPreconditioner_->setProblem(*(this->mat_), levelX_[j], levelRhs_[j]);
+        VectorType residual = *(this->rhs_);
+        Dune::MatrixVector::subtractProduct(residual, *(this->mat_), *(this->x_));
+
+        VectorType localResidual;
+        mgTransfer_[j]->restrict(residual, localResidual);
+
+        coarseGlobalPreconditioner_->setProblem(*(this->mat_), levelX_[j], localResidual);
         coarseGlobalPreconditioner_->iterate();
         const VectorType& it = coarseGlobalPreconditioner_->getSol();
 
@@ -256,12 +235,18 @@ private:
     }
 
     void iterateStep(const size_t i, const MultDirection dir) {
-        mgTransfer_[i]->restrict(*(this->x_), levelX_[i]);
+        //mgTransfer_[i]->restrict(*(this->x_), levelX_[i]);
 
+        VectorType residual = *(this->rhs_);
+        Dune::MatrixVector::subtractProduct(residual, *(this->mat_), *(this->x_));
+
+        VectorType localResidual;
+        mgTransfer_[i]->restrict(residual, localResidual);
         //print(levelRhs_[i], "levelLocalCoarseRhs: ");
         //writeToVTK(levelPatchPreconditioners_[i]->basis(), levelRhs_[i], "/storage/mi/podlesny/data/faultnetworks/rhs/fine", "exactvertexdata_step_"+std::to_string(itCounter_));
 
-        levelPatchPreconditioners_[i]->setProblem(*(this->mat_), levelX_[i], levelRhs_[i]);
+
+        levelPatchPreconditioners_[i]->setProblem(*(this->mat_), levelX_[i], localResidual);
 
         levelPatchPreconditioners_[i]->setDirection(dir);
         levelPatchPreconditioners_[i]->iterate();
diff --git a/src/cantorfaultnetworks/sparsecantorfaultnetwork.cc b/src/cantorfaultnetworks/sparsecantorfaultnetwork.cc
index 4ecac4b..b359981 100644
--- a/src/cantorfaultnetworks/sparsecantorfaultnetwork.cc
+++ b/src/cantorfaultnetworks/sparsecantorfaultnetwork.cc
@@ -153,7 +153,7 @@ int main(int argc, char** argv) { try
     const auto& interfaceNetwork = faultFactory.interfaceNetwork();
 
     if (problemCount==0) {
-        std::vector<int> writeLevels {0, 2, 4, 6};
+        std::vector<int> writeLevels {0, 2};
         for (size_t i=0; i<writeLevels.size(); i++) {
             int writeLevel = writeLevels[i];
             bool writeGrid = !(writeLevel==8);
@@ -194,7 +194,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);
-            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) * std::pow(2.0, k);
             coarseLocalInterfaceAssemblers[i] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
         }
 
@@ -235,7 +235,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);
-        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) * std::pow(2.0, k);
         exactLocalInterfaceAssemblers[i] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
     }
 
@@ -263,7 +263,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);
-        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) * std::pow(2.0, k);
         fineLocalInterfaceAssemblers[i] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
     }
 
@@ -278,7 +278,10 @@ int main(int argc, char** argv) { try
 
     std::cout << "Setting up the preconditioner!" << std::endl;
 
-    BitVector activeLevels(fineLevelIdx+1, true);
+    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;
@@ -294,7 +297,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);
-            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) * std::pow(2.0, k);
             levelLocalIntersectionAssemblers[j] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
         }
     }
diff --git a/src/cantorfaultnetworks/sparsecantorfaultnetwork.parset b/src/cantorfaultnetworks/sparsecantorfaultnetwork.parset
index 1f8e671..f5370c0 100644
--- a/src/cantorfaultnetworks/sparsecantorfaultnetwork.parset
+++ b/src/cantorfaultnetworks/sparsecantorfaultnetwork.parset
@@ -2,7 +2,7 @@ path        = ../data/
 resultPath  = ../cantorfaultnetworks/results/sparse/
 
 [preconditioner]
-patch = SUPPORT            # CELL , SUPPORT
+patch = CELL            # CELL , SUPPORT
 mode = ADDITIVE         # ADDITIVE , MULTIPLICATIVE
 multDirection = SYMMETRIC    # SYMMETRIC , FORWARD , BACKWARD
 patchDepth = 1
@@ -10,7 +10,7 @@ patchDepth = 1
 ###########################################
 
 [problem0]
-oscDataFile         = oscDataLaplace32.mat
+oscDataFile         = oscDataLaplace4.mat
 
 # level resolution in 2^(-...)
 coarseCantorLevel = 1
diff --git a/src/data/oscDataLaplace16.mat b/src/data/oscDataLaplace16.mat
new file mode 100644
index 0000000..e8d945a
--- /dev/null
+++ b/src/data/oscDataLaplace16.mat
@@ -0,0 +1,18 @@
+16 16
+1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1 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/data/oscDataLaplace8.mat b/src/data/oscDataLaplace8.mat
new file mode 100644
index 0000000..a30360f
--- /dev/null
+++ b/src/data/oscDataLaplace8.mat
@@ -0,0 +1,9 @@
+8 8 
+1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1
-- 
GitLab