A templated object for real/complex vector system of unsteady equations. More...

#include <PDE_IBVP.h>

Inheritance diagram for CppNoddy::PDE_IBVP< _Type >:

Public Member Functions
	PDE_IBVP (Equation_2matrix< _Type > equation_ptr, const DenseVector< double > &nodes, Residual_with_coords< _Type > ptr_to_left_residual, Residual_with_coords< _Type > *ptr_to_right_residual)
	The class is defined by a vector function for the system. More...

	~PDE_IBVP ()
	Destructor. More...

void	step2 (const double &dt)
	A Crank-Nicolson 'time' stepper. More...

void	assemble_matrix_problem (BandedMatrix< _Type > &a, DenseVector< _Type > &b, const double &dt)
	Assembles the matrix problem for a BVP solve at the current time level. More...

double &	coord ()
	Return a reference to the current value of the 'timelike/parabolic' coordinate. More...

OneD_Node_Mesh< _Type > &	solution ()

OneD_Node_Mesh< _Type > &	previous_solution ()

double &	tolerance ()
	Access method to the tolerance. More...

int &	max_itns ()
	Access method to the maximum number of iterations. More...

Equation_2matrix< _Type > *	get_p_equation ()

Residual_with_coords< _Type > *	get_p_left ()

Residual_with_coords< _Type > *	get_p_right ()

Detailed Description

template<typename _Type>
class CppNoddy::PDE_IBVP< _Type >

A templated object for real/complex vector system of unsteady equations.

Definition at line 37 of file PDE_IBVP.h.

Constructor & Destructor Documentation

◆ PDE_IBVP()

template<typename _Type >

CppNoddy::PDE_IBVP< _Type >::PDE_IBVP	(	Equation_2matrix< _Type > *	equation_ptr,
		const DenseVector< double > &	nodes,
		Residual_with_coords< _Type > *	ptr_to_left_residual,
		Residual_with_coords< _Type > *	ptr_to_right_residual
	)

The class is defined by a vector function for the system.

Parameters

equation_ptr	A pointer to an inherited Equation object.
nodes	A vector that defines the nodal positions.
ptr_to_left_residual	A pointer to a residual object that defines the LHS boundary conditions.
ptr_to_right_residual	A pointer to a residual object that defines the RHS boundary conditions.

Definition at line 34 of file PDE_IBVP.cpp.

                                                                                :
    TOL(1.e-8),
    T(0.0),
    MAX_ITERATIONS(12),
    p_EQUATION(ptr_to_equation),
    p_LEFT_RESIDUAL(ptr_to_left_residual),
    p_RIGHT_RESIDUAL(ptr_to_right_residual) {
    SOLN = OneD_Node_Mesh<_Type>(nodes, p_EQUATION -> get_order());
    PREV_SOLN = SOLN;
    //
    p_EQUATION -> coord(1) = T;
    if(p_EQUATION -> get_order() - p_LEFT_RESIDUAL -> get_order() - p_RIGHT_RESIDUAL -> get_order() != 0) {
      std::string problem("\n The PDE_IBVP class has been constructed, but the order of the \n");
      problem += "system does not match the number of boundary conditions.\n";
      throw ExceptionRuntime(problem);
    }
#ifdef TIME
    // timers
    T_ASSEMBLE = Timer("Assembling of the matrix (incl. equation updates):");
    T_SOLVE = Timer("Solving of the matrix:");
#endif
  }

References CppNoddy::PDE_IBVP< _Type >::coord().

◆ ~PDE_IBVP()

template<typename _Type >

CppNoddy::PDE_IBVP< _Type >::~PDE_IBVP

Destructor.

Definition at line 61 of file PDE_IBVP.cpp.

                             {
#ifdef TIME
    std::cout << "\n";
    T_ASSEMBLE.stop();
    T_ASSEMBLE.print();
    T_SOLVE.stop();
    T_SOLVE.print();
#endif
  }

Member Function Documentation

◆ assemble_matrix_problem()

template<typename _Type >

void CppNoddy::PDE_IBVP< _Type >::assemble_matrix_problem	(	BandedMatrix< _Type > &	a,
		DenseVector< _Type > &	b,
		const double &	dt
	)

Assembles the matrix problem for a BVP solve at the current time level.

Parameters

a	The LHS (banded) matrix.
b	The RHS (dense) vector.
dt	The 'time step' to be taken.

Definition at line 143 of file PDE_IBVP.cpp.

                                                                                                               {
    // clear the Jacobian matrix
    a.assign(0.0);
    // inverse of the time step
    const double inv_dt(1. / dt);
    // the order of the problem
    const unsigned order(p_EQUATION -> get_order());
    // number of spatial nodes
    const unsigned ny(SOLN.get_nnodes());
    // row counter
    std::size_t row(0);
    // a matrix that is used in the Jacobian of the mass matrix terms
    DenseMatrix<_Type> h0(order, order, 0.0);
    DenseMatrix<_Type> h1(order, order, 0.0);
    // local state variable and functions
    DenseVector<_Type> F_midpt(order, 0.0);
    DenseVector<_Type> O_midpt(order, 0.0);
    DenseVector<_Type> state(order, 0.0);
    DenseVector<_Type> state_dt(order, 0.0);
    DenseVector<_Type> state_dy(order, 0.0);
    // BCn equation is evaluated at the next time step
    p_LEFT_RESIDUAL -> coord(0) = T + dt;
    // update the BC residuals for the current iteration
    p_LEFT_RESIDUAL -> update(SOLN.get_nodes_vars(0));
    // add the (linearised) LHS BCs to the matrix problem
    for(unsigned i = 0; i < p_LEFT_RESIDUAL -> get_order(); ++i) {
      // loop thru variables at LHS of the domain
      for(unsigned var = 0; var < order; ++var) {
        a(row, var) = p_LEFT_RESIDUAL -> jacobian()(i, var);
      }
      b[ row ] = - p_LEFT_RESIDUAL -> residual()[ i ];
      ++row;
    }
    // inner nodes of the mesh, node = 0,1,2,...,N-2
    for(std::size_t node = 0; node <= ny - 2; ++node) {
      const std::size_t l_node = node;
      const std::size_t r_node = node + 1;
      // inverse of step size
      const double inv_dy = 1. / (SOLN.coord(r_node) - SOLN.coord(l_node));
      // set the current solution at this node by 2nd order evaluation at mid point
      for(unsigned var = 0; var < order; ++var) {
        const _Type F_midpt = (SOLN(l_node, var) + SOLN(r_node, var)) / 2.;
        const _Type O_midpt = (PREV_SOLN(l_node, var) + PREV_SOLN(r_node, var)) / 2.;
        state_dy[ var ] = (SOLN(r_node, var) - SOLN(l_node, var)
                           + PREV_SOLN(r_node, var) - PREV_SOLN(l_node, var)) * inv_dy / 2.;
        state[ var ] = (F_midpt + O_midpt) / 2.;
        state_dt[ var ] = (F_midpt - O_midpt) * inv_dt;
      }
      // set the equation's y & t values to be mid points
      // y
      p_EQUATION -> coord(0) = 0.5 * (SOLN.coord(l_node) + SOLN.coord(r_node));
      // t
      p_EQUATION -> coord(1) = T + dt / 2;
      // Update the equation to the mid point position
      p_EQUATION -> update(state);
      // evaluate the Jacobian of mass contribution multiplied by state_dy
      p_EQUATION -> get_jacobian_of_matrix0_mult_vector(state, state_dy, h0);
      // evaluate the Jacobian of mass contribution multiplied by state_dt
      p_EQUATION -> get_jacobian_of_matrix1_mult_vector(state, state_dt, h1);
      // loop over all the variables
      //
      // to avoid repeated mapping arithmetic within operator() of the
      // BandedMatrix class we'll access the matrix with the iterator.
      //
      // il = position for (0, lnode*order)
      typename BandedMatrix<_Type>::elt_iter l_iter(a.get_elt_iter(0, l_node * order));
      // ir = position for (0, rnode*order)
      typename BandedMatrix<_Type>::elt_iter r_iter(a.get_elt_iter(0, r_node * order));
      //
      for(unsigned var = 0; var < order; ++var) {
        // add the matrix mult terms to the linearised problem
        for(unsigned i = 0; i < order; ++i) {   // dummy index
          // add the Jacobian terms
          // indirect access: a( row, order * l_node + i ) -= jac_midpt( var, i ) * 0.5;
          const std::size_t offset(i * a.noffdiag() * 3 + row);
          const typename BandedMatrix<_Type>::elt_iter left(l_iter + offset);
          const typename BandedMatrix<_Type>::elt_iter right(r_iter + offset);
          //
          *left -= p_EQUATION -> jacobian()(var, i) * 0.5;
          // indirect access: a( row, order * r_node + i ) -= jac_midpt( var, i ) * 0.5;
          *right -= p_EQUATION -> jacobian()(var, i) * 0.5;
          // add the Jacobian of mass terms for dt terms
          // indirect access: a( row, order * l_node + i ) += h1( var, i ) * 0.5;
          *left += h0(var, i) * 0.5;
          // indirect access: a( row, order * r_node + i ) += h1( var, i ) * 0.5;
          *right += h0(var, i) * 0.5;
          // add the Jacobian of mass terms for dt terms
          // indirect access: a( row, order * l_node + i ) += h2( var, i ) * 0.5;
          *left += h1(var, i) * 0.5;
          // indirect access: a( row, order * r_node + i ) += h2( var, i ) * 0.5;
          *right += h1(var, i) * 0.5;
          // add the mass matrix terms
          // indirect access: a( row, order * l_node + i ) -= mass_midpt( var, i ) * inv_dy;
          *left -= p_EQUATION -> matrix0()(var, i) * inv_dy;
          // indirect access: a( row, order * r_node + i ) += mass_midpt( var, i ) * inv_dy;
          *right += p_EQUATION -> matrix0()(var, i) * inv_dy;
          // indirect access: a( row, order * l_node + i ) += mass_midpt( var, i ) * inv_dt;
          *left += p_EQUATION -> matrix1()(var, i) * inv_dt;
          // indirect access: a( row, order * r_node + i ) += mass_midpt( var, i ) * inv_dt;
          *right += p_EQUATION -> matrix1()(var, i) * inv_dt;
        }
        // RHS
        b[ row ] = p_EQUATION -> residual()[ var ];
        b[ row ] -= Utility::dot(p_EQUATION -> matrix0()[ var ], state_dy);
        b[ row ] -= Utility::dot(p_EQUATION -> matrix1()[ var ], state_dt);
        b[ row ] *= 2;
        // increment the row
        row += 1;
      }
    }
    // BCn equation is evaluated at the next step time point as
    // they cannot depend on d/dt terms
    p_RIGHT_RESIDUAL -> coord(0) = T + dt;
    // update the BC residuals for the current iteration
    p_RIGHT_RESIDUAL -> update(SOLN.get_nodes_vars(ny - 1));
    // add the (linearised) RHS BCs to the matrix problem
    for(unsigned i = 0; i < p_RIGHT_RESIDUAL -> get_order(); ++i) {
      // loop thru variables at RHS of the domain
      for(unsigned var = 0; var < order; ++var) {
        a(row, order * (ny - 1) + var) = p_RIGHT_RESIDUAL -> jacobian()(i, var);
      }
      b[ row ] = - p_RIGHT_RESIDUAL -> residual()[ i ];
      ++row;
    }
#ifdef PARANOID
    if(row != ny * order) {
      std::string problem("\n The ODE_BVP has an incorrect number of boundary conditions. \n");
      throw ExceptionRuntime(problem);
    }
#endif
  }

References CppNoddy::BandedMatrix< _Type >::assign(), CppNoddy::Utility::dot(), CppNoddy::BandedMatrix< _Type >::get_elt_iter(), and CppNoddy::BandedMatrix< _Type >::noffdiag().

◆ coord()

template<typename _Type >

double & CppNoddy::PDE_IBVP< _Type >::coord ( )

inline

Return a reference to the current value of the 'timelike/parabolic' coordinate.

Returns: A handle to the current time stored in the object

Definition at line 65 of file PDE_IBVP.h.

                    {
      return T;
    }

Referenced by main(), and CppNoddy::PDE_IBVP< _Type >::PDE_IBVP().

◆ get_p_equation()

template<typename _Type >

Equation_2matrix< _Type > * CppNoddy::PDE_IBVP< _Type >::get_p_equation ( )

inline

Definition at line 87 of file PDE_IBVP.h.

                                              {
      return p_EQUATION;
    }

◆ get_p_left()

template<typename _Type >

Residual_with_coords< _Type > * CppNoddy::PDE_IBVP< _Type >::get_p_left ( )

inline

Definition at line 91 of file PDE_IBVP.h.

                                              {
      return p_LEFT_RESIDUAL;
    }

◆ get_p_right()

template<typename _Type >

Residual_with_coords< _Type > * CppNoddy::PDE_IBVP< _Type >::get_p_right ( )

inline

Definition at line 95 of file PDE_IBVP.h.

                                               {
      return p_RIGHT_RESIDUAL;
    }

◆ max_itns()

template<typename _Type >

int & CppNoddy::PDE_IBVP< _Type >::max_itns ( )

inline

Access method to the maximum number of iterations.

Returns: A handle to the private member data MAX_ITERATIONS

Definition at line 83 of file PDE_IBVP.h.

                    {
      return MAX_ITERATIONS;
    }

◆ previous_solution()

template<typename _Type >

OneD_Node_Mesh< _Type > & CppNoddy::PDE_IBVP< _Type >::previous_solution

inline

Returns: A handle to the previous step's solution mesh

Definition at line 132 of file PDE_IBVP.h.

                                                                   {
    return PREV_SOLN;
  }

◆ solution()

template<typename _Type >

OneD_Node_Mesh< _Type > & CppNoddy::PDE_IBVP< _Type >::solution

inline

Returns: A handle to the solution mesh

Definition at line 127 of file PDE_IBVP.h.

                                                          {
    return SOLN;
  }

Referenced by main().

◆ step2()

template<typename _Type >

void CppNoddy::PDE_IBVP< _Type >::step2 ( const double & dt )

A Crank-Nicolson 'time' stepper.

Definition at line 72 of file PDE_IBVP.cpp.

                                              {
    // the order of the problem
    unsigned order(p_EQUATION -> get_order());
    // get the number of nodes in the mesh
    // -- this may have been refined by the user since the last call.
    unsigned ny(SOLN.get_nnodes());
    // measure of maximum residual
    double max_residual(1.0);
    // iteration counter
    int counter = 0;
    // store this soln as the 'previous SOLUTION'
    PREV_SOLN = SOLN;
    // ANY LARGE STORAGE USED IN THE MAIN LOOP IS
    // DEFINED HERE TO AVOID REPEATED CONSTRUCTION.
    // Note we blank the A matrix after every iteration.
    //
    // Banded LHS matrix - max obove diagonal band width is
    // from first variable at node i to last variable at node i+1
    BandedMatrix<_Type> a(ny * order, 2 * order - 1, 0.0);
    // RHS
    DenseVector<_Type> b(ny * order, 0.0);
    // linear solver definition
#ifdef LAPACK
    BandedLinearSystem<_Type> system(&a, &b, "lapack");
#else
    BandedLinearSystem<_Type> system(&a, &b, "native");
#endif
    // loop until converged or too many iterations
    do {
      // iteration counter
      ++counter;
#ifdef TIME
      T_ASSEMBLE.start();
#endif
      assemble_matrix_problem(a, b, dt);
      max_residual = b.inf_norm();
#ifdef DEBUG
      std::cout << " PDE_IBVP.solve : Residual_max = " << max_residual << " tol = " << TOL << "\n";
#endif
#ifdef TIME
      T_ASSEMBLE.stop();
      T_SOLVE.start();
#endif
      // linear solver
      system.solve();
      // keep the solution in a OneD_GenMesh object
      for(std::size_t var = 0; var < order; ++var) {
        for(std::size_t i = 0; i < ny; ++i) {
          SOLN(i, var) += b[ i * order + var ];
        }
      }
#ifdef TIME
      T_SOLVE.stop();
#endif
    } while((max_residual > TOL) && (counter < MAX_ITERATIONS));
    if(max_residual > TOL) {
      // restore to previous state because this step failed
      SOLN = PREV_SOLN;
      std::string problem("\n The PDE_IBVP.step2 method took too many iterations.\n");
      problem += "Solution has been restored to the previous accurate state.\n";
      throw ExceptionItn(problem, counter, max_residual);
    }
#ifdef DEBUG
    std::cout << "[DEBUG] time-like variable = " << T << "\n";
#endif
    // set the time to the updated level
    T += dt;
  }

References CppNoddy::DenseVector< _Type >::inf_norm(), and CppNoddy::BandedLinearSystem< _Type >::solve().

Referenced by main().

◆ tolerance()

template<typename _Type >

double & CppNoddy::PDE_IBVP< _Type >::tolerance ( )

inline

Access method to the tolerance.

Returns: A handle to the private member data TOL

Definition at line 77 of file PDE_IBVP.h.

                        {
      return TOL;
    }

The documentation for this class was generated from the following files:

include/PDE_IBVP.h
src/PDE_IBVP.cpp

Public Member Functions

Detailed Description

Constructor & Destructor Documentation

◆ PDE_IBVP()

◆ ~PDE_IBVP()

Member Function Documentation

◆ assemble_matrix_problem()

◆ coord()

◆ get_p_equation()

◆ get_p_left()

◆ get_p_right()

◆ max_itns()

◆ previous_solution()

◆ solution()

◆ step2()

◆ tolerance()