diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc index 8877245a776297b6fb5394fb1f6f3d5f2d6b290d..c1bbbba4a7bc30395d8ed85e3e33476978c2e7f9 100644 --- a/src/one-body-sample.cc +++ b/src/one-body-sample.cc @@ -269,6 +269,54 @@ int main(int argc, char *argv[]) { assemble_nonlinearity(grid.size(grid.maxLevel(), dim), parset, myGlobalNonlinearity, nodalIntegrals); + // {{{ Set up TNNMG solver + auto linearBaseSolverStep = + new TruncatedBlockGSStep<OperatorType, VectorType>; + auto const baseEnergyNorm = + new EnergyNorm<OperatorType, VectorType>(*linearBaseSolverStep); + auto linearBaseSolver = new LoopSolver<VectorType>( + linearBaseSolverStep, solver_maxIterations, solver_tolerance, + baseEnergyNorm, Solver::QUIET); + + // {{{ Smoothers + auto linearPresmoother = new TruncatedBlockGSStep<OperatorType, VectorType>; + auto linearPostsmoother = + new TruncatedBlockGSStep<OperatorType, VectorType>; + // }}} + + auto linearIterationStep = new MultigridStep<OperatorType, VectorType>; + linearIterationStep->setNumberOfLevels(levels); + linearIterationStep->setMGType(parset.get<int>("solver.tnnmg.linear.mu"), + parset.get<int>("solver.tnnmg.linear.nu1"), + parset.get<int>("solver.tnnmg.linear.nu2")); + linearIterationStep->basesolver_ = linearBaseSolver; + linearIterationStep->setSmoother(linearPresmoother, linearPostsmoother); + + // {{{ Transfer operators + std::vector<CompressedMultigridTransfer<VectorType> *> transferOperators; + transferOperators.resize(refinements); + for (auto &x : transferOperators) + x = new CompressedMultigridTransfer<VectorType>; + + TransferOperatorAssembler<GridType> transferOperatorAssembler(grid); + transferOperatorAssembler.assembleOperatorPointerHierarchy( + transferOperators); + linearIterationStep->setTransferOperators(transferOperators); + // }}} + + typedef GenericNonlinearGS<MyBlockProblemType> NonlinearSmootherType; + typedef TruncatedNonsmoothNewtonMultigrid< + MyBlockProblemType, NonlinearSmootherType> TNNMGStepType; + + auto nonlinearSmoother = new NonlinearSmootherType; + auto multigridStep = + new TNNMGStepType(*linearIterationStep, *nonlinearSmoother); + multigridStep->setSmoothingSteps(parset.get<int>("solver.tnnmg.main.nu1"), + parset.get<int>("solver.tnnmg.main.mu"), + parset.get<int>("solver.tnnmg.main.nu2")); + multigridStep->ignoreNodes_ = &ignoreNodes; + // }}} + for (size_t run = 1; run <= parset.get<size_t>("timesteps"); ++run) { if (parset.get<bool>("printProgress")) { std::cout << "."; @@ -303,81 +351,16 @@ int main(int argc, char *argv[]) { u1 += u1_diff; if (parset.get<bool>("useTNNMG")) { - // {{{ Linear Solver; - auto linearBaseSolverStep = - new TruncatedBlockGSStep<OperatorType, VectorType>; - auto const baseEnergyNorm = - new EnergyNorm<OperatorType, VectorType>(*linearBaseSolverStep); - auto linearBaseSolver = new LoopSolver<VectorType>( - linearBaseSolverStep, solver_maxIterations, solver_tolerance, - baseEnergyNorm, Solver::QUIET); - // }}} - - // {{{ Smoothers - auto linearPresmoother = - new TruncatedBlockGSStep<OperatorType, VectorType>; - auto linearPostsmoother = - new TruncatedBlockGSStep<OperatorType, VectorType>; - // }}} - - auto linearIterationStep = new MultigridStep<OperatorType, VectorType>; - linearIterationStep->setNumberOfLevels(levels); - linearIterationStep->setMGType( - parset.get<int>("solver.tnnmg.linear.mu"), - parset.get<int>("solver.tnnmg.linear.nu1"), - parset.get<int>("solver.tnnmg.linear.nu2")); - linearIterationStep->basesolver_ = linearBaseSolver; - linearIterationStep->setSmoother(linearPresmoother, linearPostsmoother); - - // {{{ Transfer operators - std::vector<CompressedMultigridTransfer<VectorType> *> - transferOperators; - transferOperators.resize(refinements); - for (auto &x : transferOperators) - x = new CompressedMultigridTransfer<VectorType>; - - TransferOperatorAssembler<GridType> transferOperatorAssembler(grid); - transferOperatorAssembler.assembleOperatorPointerHierarchy( - transferOperators); - linearIterationStep->setTransferOperators(transferOperators); - // }}} - - typedef GenericNonlinearGS<MyBlockProblemType> NonlinearSmootherType; - typedef TruncatedNonsmoothNewtonMultigrid< - MyBlockProblemType, NonlinearSmootherType> TNNMGStepType; - MyConvexProblemType myConvexProblem(stiffnessMatrix, *myGlobalNonlinearity, b4, u4_diff); auto myBlockProblem = new MyBlockProblemType(parset, myConvexProblem); - auto nonlinearSmoother = new NonlinearSmootherType; - auto multigridStep = - new TNNMGStepType(*linearIterationStep, *nonlinearSmoother); multigridStep->setProblem(u4_diff, *myBlockProblem); - multigridStep->setSmoothingSteps( - parset.get<int>("solver.tnnmg.main.nu1"), - parset.get<int>("solver.tnnmg.main.mu"), - parset.get<int>("solver.tnnmg.main.nu2")); - multigridStep->ignoreNodes_ = &ignoreNodes; - - auto const energyNorm = - new EnergyNorm<OperatorType, VectorType>(stiffnessMatrix); LoopSolver<VectorType> overallSolver( - multigridStep, solver_maxIterations, solver_tolerance, energyNorm, + multigridStep, solver_maxIterations, solver_tolerance, &energyNorm, verbosity); overallSolver.solve(); - delete linearBaseSolverStep; - delete baseEnergyNorm; - delete linearBaseSolver; - delete linearPresmoother; - delete linearPostsmoother; - for (auto &x : transferOperators) - delete x; - delete linearIterationStep; delete myBlockProblem; - delete nonlinearSmoother; - delete multigridStep; - delete energyNorm; } u4 += u4_diff; @@ -434,6 +417,19 @@ int main(int argc, char *argv[]) { } std::cout << std::endl; + // {{{ Clean up after the TNNMG solver + delete linearBaseSolverStep; + delete baseEnergyNorm; + delete linearBaseSolver; + delete linearPresmoother; + delete linearPostsmoother; + for (auto &x : transferOperators) + delete x; + delete linearIterationStep; + delete nonlinearSmoother; + delete multigridStep; + // }}} + if (parset.get<bool>("printDifference")) { VectorType diff2 = u2; diff2 -= u1;