From: Joscha Gedicke Date: Mon, 2 May 2016 20:45:32 +0000 (+0200) Subject: fix overlapping block Jacobi in relaxation preconditioner X-Git-Tag: v8.5.0-rc1~1052^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=3ed9da792305f037c1f9969a1b359ffd937169e2;p=dealii.git fix overlapping block Jacobi in relaxation preconditioner --- diff --git a/doc/news/changes.h b/doc/news/changes.h index adbfeb3c43..241540efc4 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -233,6 +233,14 @@ inconvenience this causes.

Specific improvements

    +
  1. Fixed: Bug in the RelaxationBlock class function do_step. Before, the + corrections were not added together, which leads to a wrong update whenever the + Jacobi blocks are overlapping. For SOR, SSOR and non-overlapping Jacobi this was + not an issue. +
    + (Joscha Gedicke, 2016/05/07) +
    2. +
  2. New: Added function GridOut::write_mesh_per_processor_as_vtu. This allows the visualization of a parallel finite element mesh that can be separated into each processor's owned and ghost cells. It also allows for the visualization of each level diff --git a/include/deal.II/lac/relaxation_block.templates.h b/include/deal.II/lac/relaxation_block.templates.h index ce4174f01d..27c7254427 100644 --- a/include/deal.II/lac/relaxation_block.templates.h +++ b/include/deal.II/lac/relaxation_block.templates.h @@ -210,7 +210,7 @@ RelaxationBlock::do_step (Vector &dst, // Store in result vector row=additional_data->block_list.begin(block); for (size_type row_cell=0; row_cellcolumn()) = prev(row->column()) + additional_data->relaxation * x_cell(row_cell); + dst(row->column()) += additional_data->relaxation * x_cell(row_cell); } } } diff --git a/tests/lac/solver_relaxation_04.cc b/tests/lac/solver_relaxation_04.cc new file mode 100644 index 0000000000..cbe98890a5 --- /dev/null +++ b/tests/lac/solver_relaxation_04.cc @@ -0,0 +1,142 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 1998 - 2015 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Compare overlapping block Jacobi relaxation with different +// permutations of the blocks. All output diffs should be zero. + +#include +#include +#include +#include +#include +#include +#include + +#include "../tests.h" +#include "testmatrix.h" +#include +#include +#include + + +template +double +check_solve (SolverType &solver, + const MatrixType &A, + VectorType &u, + VectorType &f, + const PRECONDITION &P) +{ + double result = 0.; + u = 0.; + f = 1.; + try + { + solver.solve(A,u,f,P); + } + catch (SolverControl::NoConvergence &e) + { + result = e.last_residual; + } + return result; +} + +int main() +{ + const std::string logname = "output"; + std::ofstream logfile(logname.c_str()); + // logfile.setf(std::ios::fixed); + deallog << std::setprecision(4); + deallog.attach(logfile); + deallog.threshold_double(1.e-10); + + SolverControl control(10, 1.e-3); + SolverRelaxation<> relax(control); + + // Solve non-symmetric laplace with five-point FD + for (unsigned int size=33; size <= 33; size *= 3) + { + unsigned int dim = (size-1)*(size-1); + + deallog << "Size " << size << " Unknowns " << dim << std::endl; + + // Make matrix + FDMatrix testproblem(size, size); + SparsityPattern structure(dim, dim, 5); + testproblem.five_point_structure(structure); + structure.compress(); + SparseMatrix A(structure); + testproblem.five_point(A,true); + + for (unsigned int blocksize = 4; blocksize < 32; blocksize <<= 1) + { + deallog << "Block size " << blocksize << std::endl; + + const unsigned int n_blocks = dim/blocksize; + RelaxationBlock,double>::AdditionalData relax_data(0.7); + RelaxationBlock,double>::AdditionalData relax_data_reorder(0.7); + + relax_data.block_list.reinit(n_blocks, dim, blocksize+2); + relax_data_reorder.block_list.reinit(n_blocks, dim, blocksize+2); + for (unsigned int block=0; block-1 && (i+block*blocksize),double> relax_jacobi; + relax_jacobi.initialize(A, relax_data); + + // reverse the order of the blocks + relax_data_reorder.order.resize(1); + relax_data_reorder.order[0].resize(n_blocks); + for (unsigned int i=0; i,double> relax_jacobi_reorder; + relax_jacobi_reorder.initialize(A, relax_data_reorder); + + Vector f(dim); + Vector u(dim); + Vector res(dim); + + f = 1.; + u = 1.; + + try + { + double r1, r2; + + deallog.push("Jacobi"); + r1 = check_solve(relax,A,u,f,relax_jacobi); + r2 = check_solve(relax,A,u,f,relax_jacobi_reorder); + deallog << "diff " << std::fabs(r1-r2)/r1 << std::endl; + deallog.pop(); + + } + catch (std::exception &e) + { + std::cerr << "Exception: " << e.what() << std::endl; + } + } + } +} \ No newline at end of file diff --git a/tests/lac/solver_relaxation_04.output b/tests/lac/solver_relaxation_04.output new file mode 100644 index 0000000000..0ab871937d --- /dev/null +++ b/tests/lac/solver_relaxation_04.output @@ -0,0 +1,20 @@ + +DEAL::Size 33 Unknowns 1024 +DEAL::Block size 4 +DEAL:Jacobi:Relaxation::Starting value 32.00 +DEAL:Jacobi:Relaxation::Failure step 10 value 28.14 +DEAL:Jacobi:Relaxation::Starting value 32.00 +DEAL:Jacobi:Relaxation::Failure step 10 value 28.14 +DEAL:Jacobi::diff 0 +DEAL::Block size 8 +DEAL:Jacobi:Relaxation::Starting value 32.00 +DEAL:Jacobi:Relaxation::Failure step 10 value 26.83 +DEAL:Jacobi:Relaxation::Starting value 32.00 +DEAL:Jacobi:Relaxation::Failure step 10 value 26.83 +DEAL:Jacobi::diff 0 +DEAL::Block size 16 +DEAL:Jacobi:Relaxation::Starting value 32.00 +DEAL:Jacobi:Relaxation::Failure step 10 value 26.12 +DEAL:Jacobi:Relaxation::Starting value 32.00 +DEAL:Jacobi:Relaxation::Failure step 10 value 26.12 +DEAL:Jacobi::diff 0