Newton.h 3.17 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28
// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
// ==  http://www.amdis-fem.org                                              ==
// ==                                                                        ==
// ============================================================================
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.



/** \file Newton.h */

#ifndef AMDIS_NEWTON_H
#define AMDIS_NEWTON_H

#include "CreatorInterface.h"
#include "NonLinSolver.h"
#include "OEMSolver.h"
29
#include "io/VtkWriter.h"
30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48

namespace AMDiS {

  /**
   * \ingroup Solver
   * 
   * \Brief
   * Implements the newton method for solving a non linear system. Sub class of
   * NonLinSolver.
   */
  class Newton : public NonLinSolver
  {
  public:
    /// Creator class used in the NonLinSolverMap.
    class Creator : public NonLinSolverCreator
    {
    public:
      virtual ~Creator() {}

49
      /// Returns a new Newton object.
50 51
      NonLinSolver* create() 
      { 
52
	return new Newton(this->name, this->linearSolver); 
53 54 55
      }
    };

56 57 58
    /// Calls constructor of base class NonLinSolver
    Newton(const std::string& name, OEMSolver *linSolver)
      : NonLinSolver(name, linSolver),
59 60 61 62
	b(NULL)
    {}

  private:
63 64
    /// Realisation of NonLinSolver::init
    void init() {}
65

66
    /// realisation of NonLinSolver::nlsolve
67
    int nlsolve(SolverMatrix<Matrix<DOFMatrix*> >& mat,
68 69 70
		SystemVector& x, SystemVector& rhs, 
		AdaptInfo *adaptInfo,
		ProblemStat *prob)
71
    {
72
      FUNCNAME("Newton::nlsolve()");
73

74 75
      if (b == NULL)
	b = new SystemVector(x);
76

77 78
      double err = 0.0, errOld = -1.0;
      int iter, n;
79

80
      MSG("iter. |     this->residual |     red. |    n |\n");
81 82

      for (iter = 1; iter <= this->maxIter; iter++) {
83 84 85 86 87 88 89
	// Assemble DF(x) and F(x)
	prob->buildAfterCoarsen(adaptInfo, 0, true, true);

	// Initial guess is zero
	b->set(0.0);

	// Solve linear system
90
	n = solveLinearSystem(mat, *b, rhs);
91 92 93

	// x = x + d
	x += *b;
94 95

	if (this->usedNorm == NO_NORM || this->usedNorm == L2_NORM)
96
	  err = L2Norm(b);
97
	else
98
	  err = H1Norm(b);
99 100
    

101 102
	if (iter == 1)  
	  this->initialResidual = err;
103

104 105 106 107 108 109 110 111 112 113 114
	if (errOld <= 0)
	  MSG("%5d | %12.5e | -------- | %4d |\n", iter, err, n);
	else
	  MSG("%5d | %12.5e | %8.2e | %4d |\n", iter, err, err/errOld, n);	

	residual = err;
 	if (err < this->tolerance) {
 	  MSG("Finished successfully!\n");
 	  return iter;
 	}
	errOld = err;
115 116
      }

117 118
      MSG("iter. %d, residual: %12.5e\n", iter, err);
      MSG("tolerance %e not reached\n", this->tolerance);
119 120 121 122 123 124

      this->residual = err;

      return iter;
    }

125
    /// Realisation of NonLinSolver::exit
126 127
    void exit()
    {
128
      if (b != NULL) {
129
	delete b;
130 131
	b = NULL;
      }
132 133 134
    }

  private:
135
    /// Internal used data
136 137 138 139 140 141
    SystemVector *b;
  };

}

#endif // AMDIS_NEWTON_H