Skip to content
Snippets Groups Projects
Commit 15e97530 authored by podlesny's avatar podlesny
Browse files

multiplicative working

parent dbc43ca3
No related branches found
No related tags found
No related merge requests found
...@@ -98,6 +98,8 @@ public: ...@@ -98,6 +98,8 @@ public:
} }
virtual void iterate() { virtual void iterate() {
*(this->x_) = 0;
if (mode_ == ADDITIVE) if (mode_ == ADDITIVE)
iterateAdd(); iterateAdd();
else else
...@@ -151,8 +153,6 @@ private: ...@@ -151,8 +153,6 @@ private:
} }
void iterateAdd() { void iterateAdd() {
//*(this->x_) = 0;
VectorType it, x; VectorType it, x;
for (size_t i=0; i<localProblems_.size(); i++) { for (size_t i=0; i<localProblems_.size(); i++) {
localProblems_[i]->solve(it); localProblems_[i]->solve(it);
...@@ -182,44 +182,15 @@ private: ...@@ -182,44 +182,15 @@ private:
auto& localProblem = *localProblems_[i]; auto& localProblem = *localProblems_[i];
VectorType v, v1, v2, v3;
localProblem.restrict(*(this->x_), v);
v1.resize(localProblem.getLocalMat().N());
localProblem.getLocalMat().mv(v, v1);
v3.resize(matrix_.N());
matrix_.mv(*(this->x_), v3);
localProblem.restrict(v3, v2);
//print(v1, "v1");
//print(v2, "v2");
VectorType r1;
localProblem.restrict(*(this->rhs_), r1);
Dune::MatrixVector::subtractProduct(r1, localProblem.getLocalMat(), v);
/*localProblem.updateLocalRhs(r);
VectorType x;
localProblem.solve(x);
v += x;
localProblem.updateVector(v, *(this->x_));*/
VectorType localR; VectorType localR;
VectorType r = *(this->rhs_); VectorType r = *(this->rhs_);
Dune::MatrixVector::subtractProduct(r, matrix_, *(this->x_)); Dune::MatrixVector::subtractProduct(r, matrix_, *(this->x_));
localProblem.restrict(r, localR); localProblem.restrict(r, localR);
// print(r1, "r1");
// print(localR, "localR");
localProblem.updateLocalRhs(localR); localProblem.updateLocalRhs(localR);
VectorType x; VectorType x, v;
localProblem.solve(x); localProblem.solve(x);
//VectorType v;
localProblem.prolong(x, v); localProblem.prolong(x, v);
*(this->x_) += v; *(this->x_) += v;
} }
......
...@@ -141,7 +141,7 @@ int main(int argc, char** argv) { try ...@@ -141,7 +141,7 @@ int main(int argc, char** argv) { try
matrixReader.read(matrix); matrixReader.read(matrix);
const int oscGridN = matrix.N(); const int oscGridN = matrix.N();
const int oscGridLevelIdx = std::round(std::log2(oscGridN)) - coarseResolution; const int oscGridLevelIdx = std::max(std::round(std::log2(oscGridN)) - coarseResolution, 0.0);
if (oscGridN>fineGridN) if (oscGridN>fineGridN)
DUNE_THROW(Dune::Exception, "Provided oscData too large!"); DUNE_THROW(Dune::Exception, "Provided oscData too large!");
...@@ -161,7 +161,7 @@ int main(int argc, char** argv) { try ...@@ -161,7 +161,7 @@ int main(int argc, char** argv) { try
CantorFaultFactory<GridType> faultFactory(levelToCantorLevel, coarseResolution, exactLevelIdx, maxCantorLevel); CantorFaultFactory<GridType> faultFactory(levelToCantorLevel, coarseResolution, exactLevelIdx, maxCantorLevel);
const InterfaceNetwork<GridType>& interfaceNetwork = faultFactory.interfaceNetwork(); const InterfaceNetwork<GridType>& interfaceNetwork = faultFactory.interfaceNetwork();
if (problemCount==0) { /*if (problemCount==0) {
std::vector<int> writeLevels {1, 2, 8}; std::vector<int> writeLevels {1, 2, 8};
for (size_t i=0; i<writeLevels.size(); i++) { for (size_t i=0; i<writeLevels.size(); i++) {
int writeLevel = writeLevels[i]; int writeLevel = writeLevels[i];
...@@ -170,7 +170,7 @@ int main(int argc, char** argv) { try ...@@ -170,7 +170,7 @@ int main(int argc, char** argv) { try
LevelInterfaceNetworkWriter networkWriter(resultPath + "levelinterfacenetwork_" + std::to_string(writeLevel) + ".tikz"); LevelInterfaceNetworkWriter networkWriter(resultPath + "levelinterfacenetwork_" + std::to_string(writeLevel) + ".tikz");
networkWriter.write(interfaceNetwork.levelInterfaceNetwork(writeLevel), writeGrid); networkWriter.write(interfaceNetwork.levelInterfaceNetwork(writeLevel), writeGrid);
} }
} }*/
const GridType& grid = faultFactory.grid(); const GridType& grid = faultFactory.grid();
...@@ -332,23 +332,23 @@ int main(int argc, char** argv) { try ...@@ -332,23 +332,23 @@ int main(int argc, char** argv) { try
DGMGTransfer<DGBasis> discMGTransfer(*fineBasis, *exactBasis); DGMGTransfer<DGBasis> discMGTransfer(*fineBasis, *exactBasis);
// solve // solve
/*OscCGSolver<MatrixType, VectorType, DGBasis> OscCGSolver<MatrixType, VectorType, DGBasis>
solver(*fineBasis, &fineGlobalAssembler.matrix(), &initialX, &rhs, &preconditioner, solver(*fineBasis, &fineGlobalAssembler.matrix(), &initialX, &rhs, &preconditioner,
maxIterations, solverTolerance, &exactEnergyNorm, maxIterations, solverTolerance, &exactEnergyNorm,
Solver::FULL, &exactSol, &fineSol, discMGTransfer, 1.0, true); //((oscGridN+0.0)/fineGridN) */ Solver::FULL, &exactSol, &fineSol, discMGTransfer, 1.0, true); //((oscGridN+0.0)/fineGridN)
//solver.historyBuffer_ = resultPath + "solutions/level_" + std::to_string(fineResolution) + ".vec"; //solver.historyBuffer_ = resultPath + "solutions/level_" + std::to_string(fineResolution) + ".vec";
VectorType x = initialX; /*VectorType x = initialX;
preconditioner.setProblem(fineGlobalAssembler.matrix(), x, rhs); preconditioner.setProblem(fineGlobalAssembler.matrix(), x, rhs);
Dune::Solvers::LoopSolver<VectorType> Dune::Solvers::LoopSolver<VectorType>
solver (preconditioner, maxIterations, solverTolerance, solver (preconditioner, maxIterations, solverTolerance,
fineEnergyNorm, fineEnergyNorm,
Solver::FULL, Solver::FULL,
true, &fineSol); true, &fineSol);*/
solver.check(); solver.check();
solver.preprocess(); solver.preprocess();
......
...@@ -2,7 +2,7 @@ path = ../data/ ...@@ -2,7 +2,7 @@ path = ../data/
resultPath = ../cantorfaultnetworks/results/ resultPath = ../cantorfaultnetworks/results/
[preconditioner] [preconditioner]
patch = SUPPORT # CELL , SUPPORT patch = CELL # CELL , SUPPORT
mode = MULTIPLICATIVE # ADDITIVE , MULTIPLICATIVE mode = MULTIPLICATIVE # ADDITIVE , MULTIPLICATIVE
multDirection = SYMMETRIC # SYMMETRIC , FORWARD , BACKWARD multDirection = SYMMETRIC # SYMMETRIC , FORWARD , BACKWARD
patchDepth = 1 patchDepth = 1
...@@ -10,12 +10,12 @@ patchDepth = 1 ...@@ -10,12 +10,12 @@ patchDepth = 1
########################################### ###########################################
[problem0] [problem0]
oscDataFile = oscDataLaplace32.mat oscDataFile = oscDataLaplace4.mat
# level resolution in 2^(-...) # level resolution in 2^(-...)
coarseResolution = 0 coarseResolution = 0
fineResolution = 5 fineResolution = 3
exactResolution = 8 exactResolution = 4
minCantorResolution = 0 minCantorResolution = 0
penaltyFactor = 1 penaltyFactor = 1
......
...@@ -14,8 +14,8 @@ oscDataFile = oscDataLaplace4.mat ...@@ -14,8 +14,8 @@ oscDataFile = oscDataLaplace4.mat
# level resolution in 2^(-...) # level resolution in 2^(-...)
coarseCantorLevel = 1 coarseCantorLevel = 1
fineCantorLevel = 2 fineCantorLevel = 4
maxCantorLevel = 3 maxCantorLevel = 5
penaltyFactor = 1 penaltyFactor = 1
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment