NrIterator.h

#ifndef R1D_NR_ITERATOR_H
#define R1D_NR_ITERATOR_H
 
/*
* PSCF - Polymer Self-Consistent Field
*
* Copyright 2015 - 2025, The Regents of the University of Minnesota
* Distributed under the terms of the GNU General Public License.
*/
 
#include "Iterator.h"
#include <r1d/solvers/Mixture.h>
#include <pscf/math/LuSolver.h>
#include <util/containers/Array.h>
#include <util/containers/DArray.h>
#include <util/containers/DMatrix.h>
 
namespace Pscf {
namespace R1d
{
 
   using namespace Util;

   class NrIterator : public Iterator
   {
 
   public:

      typedef Mixture::FieldT FieldT;

      NrIterator();

      NrIterator(System& system);

      virtual ~NrIterator();

      void readParameters(std::istream& in);

      int solve(bool isContinuation = false);
 
   private:

      LuSolver solver_;

      DArray<FieldT> wFieldsNew_;

      DArray<FieldT> cFieldsNew_;

      DArray<double> cArray_;

      DArray<double> wArray_;

      DArray<double> residual_;

      DMatrix<double> jacobian_;

      DArray<double> residualNew_;

      DArray<double> dOmega_;

      double epsilon_;

      int maxItr_;

      bool isAllocated_;

      bool newJacobian_;

      bool needsJacobian_;

      bool isCanonical_;
 
      // Private member functions

      void setup();

      void computeResidual(Array<FieldT> const & wFields, 
                           Array<FieldT> const & cFields, 
                           Array<double>& residual);

      double residualNorm(Array<double> const & residual) const;

      void computeJacobian();

      void incrementWFields(Array<FieldT> const & wOld,
                            Array<double> const & dW,
                            Array<FieldT>& wNew);
 
   };
 
} // namespace R1d
} // namespace Pscf
#endif