Unverified Commit 3ebc783f authored by COLLINSBS email's avatar COLLINSBS email Committed by GitHub
Browse files

Implements Base Solution Acceleration Methods (#4717) (#299) 

* Initial base class implementation

* Adds picard and modified picard types

* Adds unit test for Relaxed and Modified Picard

* Updates to comments and fixes test failures

* Switch Anderson to use base Accelerator type

* Fixing unit test

* Removes unneeded DBC statement

* Addresses review comments

* Adds more DBC checking, updates comments
parent 9615fa63

src/AndersonAcceleration.f90

+79 −83
Original line number Diff line number Diff line
@@ -20,6 +20,7 @@ USE ParallelEnv 
USE VectorTypes

USE MatrixTypes

USE LinearSolverTypes

USE SolutionAccelerationModule



IMPLICIT NONE
@@ -29,21 +30,11 @@ PRIVATE 
PUBLIC :: AndersonAccelerationType



!> @brief the native Anderson Accelerator Type

TYPE :: AndersonAccelerationType

  !> Initialization status

  LOGICAL(SBK) :: isInit=.FALSE.

  !> Size of solution vector/matrix

  INTEGER(SIK) :: n=-1

  !> Current iteration count

  INTEGER(SIK) :: s=0

TYPE,EXTENDS(SolutionAccelerationType) :: AndersonAccelerationType

  !> Iteration depth of Anderson solver

  INTEGER(SIK) :: depth=1

  !> Starting iteration for Anderson

  INTEGER(SIK) :: start=1

  !> Value of mixing parameter

  REAL(SRK) :: beta=0.5_SRK

  !> Initial iterates

  CLASS(VectorType),ALLOCATABLE :: x(:)

  !> Gx vectors:

  CLASS(VectorType),ALLOCATABLE :: Gx(:)

  !> Difference vectors r=Gx-x
@@ -54,8 +45,6 @@ TYPE :: AndersonAccelerationType 
  REAL(SRK),ALLOCATABLE :: alpha(:)

  !> Linear solver

  TYPE(LinearSolverType_Direct) :: LS

  !> Futility computing environment

  TYPE(FutilityComputingEnvironment),POINTER :: ce => NULL()

!

!List of Type Bound Procedures

  CONTAINS
@@ -65,6 +54,9 @@ TYPE :: AndersonAccelerationType 
    !> @copybrief AndersonAccelerationModule::clear_AndersonAccelerationType

    !> @copydetails AndersonAccelerationModule::clear_AndersonAccelerationType

    PROCEDURE,PASS :: clear => clear_AndersonAccelerationType

    !> @copybrief AndersonAccelerationModule::setInitial_AndersonAccelerationType

    !> @copydetails AndersonAccelerationModule::setInitial_AndersonAccelerationType

    PROCEDURE,PASS :: setInitial => setInitial_AndersonAccelerationType

    !> @copybrief AndersonAccelerationModule::step_AndersonAccelerationType

    !> @copydetails AndersonAccelerationModule::step_AndersonAccelerationType

    PROCEDURE,PASS :: step => step_AndersonAccelerationType
@@ -82,8 +74,8 @@ CONTAINS 
!-------------------------------------------------------------------------------

!> @brief Initializes the Anderson Acceleration Type with an input parameter list

!> @param solver The Anderson Acceleration solver to act on

!> @param Params The Anderson options parameter list

!> @param ce The computing environment to use for the calculation

!> @param Params The Anderson options parameter list

!>

SUBROUTINE init_AndersonAccelerationType(solver,ce,Params)

  CLASS(AndersonAccelerationType),INTENT(INOUT) :: solver
@@ -94,27 +86,18 @@ SUBROUTINE init_AndersonAccelerationType(solver,ce,Params) 
  TYPE(ParamType) :: LSparams

  INTEGER(SIK) :: i,j,m



  REQUIRE(Params%has('AndersonAccelerationType->n'))

  REQUIRE(.NOT.solver%isInit)



  CALL solver%init_base(ce,Params)

  !Pull Data from Parameter List

  CALL Params%get('AndersonAccelerationType->n',solver%n)

  IF(Params%has('AndersonAccelerationType->depth')) &

      CALL Params%get('AndersonAccelerationType->depth',solver%depth)

  IF(Params%has('AndersonAccelerationType->beta')) &

      CALL Params%get('AndersonAccelerationType->beta',solver%beta)

  IF(Params%has('AndersonAccelerationType->start')) &

      CALL Params%get('AndersonAccelerationType->start',solver%start)



  IF(solver%n < 1) CALL ce%exceptHandler%raiseError('Incorrect input to '//modName// &

      '::'//myName//' - Number of unkowns (n) must be greater than 0!')

  IF(Params%has('SolutionAcceleration->anderson_depth')) &

      CALL Params%get('SolutionAcceleration->anderson_depth',solver%depth)

  IF(Params%has('SolutionAcceleration->anderson_mixing_parameter')) &

      CALL Params%get('SolutionAcceleration->anderson_mixing_parameter',solver%beta)



  IF(solver%depth < 0) CALL ce%exceptHandler%raiseError('Incorrect input to '//modName// &

      '::'//myName//' - Anderson depth must be greater than or equal to 0!')



  IF(solver%start < 1) CALL ce%exceptHandler%raiseError('Incorrect input to '//modName// &

      '::'//myName//' - Anderson starting point must be greater than 1!')



  IF(solver%beta <= 0.0_SRK .OR. solver%beta > 1.0_SRK) THEN

    CALL ce%exceptHandler%raiseError('Incorrect input to '//modName// &

        '::'//myName//' - Anderson mixing parameter must be between 0.0 and 1.0!')
@@ -150,8 +133,6 @@ SUBROUTINE init_AndersonAccelerationType(solver,ce,Params) 
  m=solver%depth+1

  ALLOCATE(solver%alpha(m))

  solver%alpha(:)=0.0_SRK

  solver%s=0

  solver%ce => ce

  solver%isInit=.TRUE.



ENDSUBROUTINE init_AndersonAccelerationType
@@ -165,32 +146,73 @@ SUBROUTINE clear_AndersonAccelerationType(solver) 


  INTEGER(SIK) :: i



  IF(solver%isInit) THEN

    IF(ALLOCATED(solver%x)) THEN

  CALL solver%clear_base()



  IF(ALLOCATED(solver%tmpvec)) THEN

    DO i=1,solver%depth+1

        CALL solver%x(i)%clear()

      CALL solver%Gx(i)%clear()

      CALL solver%r(i)%clear()

    ENDDO

    CALL solver%tmpvec%clear()

      DEALLOCATE(solver%x)

    DEALLOCATE(solver%Gx)

    DEALLOCATE(solver%r)

    DEALLOCATE(solver%tmpvec)

  ENDIF

  DEALLOCATE(solver%alpha)

  IF(solver%LS%isinit) CALL solver%LS%clear()

    solver%s=0

    solver%n=-1

  solver%depth=1

    solver%start=1

  solver%beta=0.5_SRK

  solver%isInit=.FALSE.

  ENDIF



ENDSUBROUTINE clear_AndersonAccelerationType

!

!-------------------------------------------------------------------------------

!> @brief Set the initial iterate for the Anderson solver

!> @param solver The Anderson solver to act on

!> @param x Initial iterate solve is starting from

!>

SUBROUTINE setInitial_AndersonAccelerationType(solver,x)

  CLASS(AndersonAccelerationType),INTENT(INOUT) :: solver

  CLASS(VectorType),INTENT(INOUT) :: x



#if defined(HAVE_MPI) || defined(FUTILITY_HAVE_Trilinos)

  CHARACTER(LEN=*),PARAMETER :: myName='setInitial_AndersonAccelerationType'

#endif



  REQUIRE(x%n == solver%n)

  REQUIRE(solver%s == 0)



  !If this is the first call to set/reset must actually create vectors of needed type

  IF(.NOT.ALLOCATED(solver%x)) THEN



    !!!TODO Once missing BLAS interfaces have been added for the following vector types

    !do away with this error catch to allow use of these with Anderson Acceleration.



    SELECT TYPE(x)

#ifdef HAVE_MPI

    TYPE IS(NativeDistributedVectorType)

      CALL solver%ce%exceptHandler%raiseError('Incorrect call to '//modName// &

          '::'//myName//' - Input vector type not supported!')

#endif

#ifdef FUTILITY_HAVE_Trilinos

    TYPE IS(TrilinosVectorType)

      CALL solver%ce%exceptHandler%raiseError('Incorrect call to '//modName// &

          '::'//myName//' - Input vector type not supported!')

#endif

    ENDSELECT



    CALL VectorResembleAlloc(solver%x,x,solver%depth+1)

    CALL VectorResembleAlloc(solver%Gx,x,solver%depth+1)

    CALL VectorResembleAlloc(solver%r,x,solver%depth+1)

    CALL VectorResembleAlloc(solver%tmpvec,x)

  ENDIF



  !Grab initial iterate to initiate Anderson from

  CALL BLAS_copy(x,solver%x(1))



ENDSUBROUTINE setInitial_AndersonAccelerationType

!

!-------------------------------------------------------------------------------

!> @brief Performs a single step of Anderson Acceleration acting on the input solution vector.

!>        If depth is set to zero, or if the iteration count is below the Anderson starting

!>        point the behavior is to under-relax the solution using the mixing parameter (beta)
@@ -206,6 +228,10 @@ SUBROUTINE step_AndersonAccelerationType(solver,x_new) 
  REAL(SRK) :: tmpA,tmpb



  REQUIRE(x_new%n == solver%n)

  REQUIRE(ALLOCATED(solver%x))

  REQUIRE(ALLOCATED(solver%Gx))

  REQUIRE(ALLOCATED(solver%r))

  REQUIRE(ALLOCATED(solver%tmpvec))



  !Update iteration counter

  solver%s=solver%s+1
@@ -289,7 +315,7 @@ SUBROUTINE step_AndersonAccelerationType(solver,x_new) 
ENDSUBROUTINE step_AndersonAccelerationType

!

!-------------------------------------------------------------------------------

!> @brief Set or reset the initial iterate and linear solver state for the Anderson solver

!> @brief Reset the initial iterate and linear solver state for the Anderson solver

!> @param solver The Anderson solver to act on

!> @param x Initial iterate solve is starting from

!>
@@ -297,43 +323,13 @@ SUBROUTINE reset_AndersonAccelerationType(solver,x) 
  CLASS(AndersonAccelerationType),INTENT(INOUT) :: solver

  CLASS(VectorType),INTENT(INOUT) :: x



#if defined(HAVE_MPI) || defined(FUTILITY_HAVE_Trilinos)

  CHARACTER(LEN=*),PARAMETER :: myName='reset_AndersonAccelerationType'

#endif

  INTEGER(SIK) :: i,j



  REQUIRE(x%n == solver%n)



  !If this is the first call to set/reset must actually create vectors of needed type

  IF(.NOT.ALLOCATED(solver%x)) THEN



    !!!TODO Once missing BLAS interfaces have been added for the following vector types

    !do away with this error catch to allow use of these with Anderson Acceleration.



    SELECT TYPE(x)

#ifdef HAVE_MPI

    TYPE IS(NativeDistributedVectorType)

      CALL solver%ce%exceptHandler%raiseError('Incorrect call to '//modName// &

          '::'//myName//' - Input vector type not supported!')

#endif

#ifdef FUTILITY_HAVE_Trilinos

    TYPE IS(TrilinosVectorType)

      CALL solver%ce%exceptHandler%raiseError('Incorrect call to '//modName// &

          '::'//myName//' - Input vector type not supported!')

#endif

    ENDSELECT



    CALL VectorResembleAlloc(solver%x,x,solver%depth+1)

    CALL VectorResembleAlloc(solver%Gx,x,solver%depth+1)

    CALL VectorResembleAlloc(solver%r,x,solver%depth+1)

    CALL VectorResembleAlloc(solver%tmpvec,x)

  ENDIF



  !Reset iteration counter

  solver%s=0_SIK



  !Grab initial iterate to initiate Anderson from

  CALL BLAS_copy(x,solver%x(1))

  CALL solver%setInitial(x)



  !Reset Anderson coefficient matrix

  IF(solver%LS%isinit) THEN

src/CMakeLists.txt

+1 −0
Original line number Diff line number Diff line
@@ -92,6 +92,7 @@ TRIBITS_ADD_LIBRARY(Utils 
      EigenvalueSolverTypes.f90

      ODESolverTypes.f90

      NonLinearSolver.f90

      SolutionAcceleration.f90

      AndersonAcceleration.f90

      StochasticSampling.f90

      Sorting.f90

src/SolutionAcceleration.f90

0 → 100644
+390 −0
Original line number Diff line number Diff line 
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!

!                          Futility Development Group                          !

!                             All rights reserved.                             !

!                                                                              !

! Futility is a jointly-maintained, open-source project between the University !

! of Michigan and Oak Ridge National Laboratory.  The copyright and license    !

! can be found in LICENSE.txt in the head directory of this repository.        !

!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!

!> @brief A Fortran 2003 module that encompasses a base class for solution

!>        acceleration.  Module also includes both picard iteration and a

!>        modified picard iteration.  Anderson Acceleration extends this base

!>        class in the AndersonAcceleration module.

!>

!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!

MODULE SolutionAccelerationModule

#include "Futility_DBC.h"

USE Futility_DBC

USE IntrType

USE FutilityComputingEnvironmentModule

USE ParameterLists

USE ParallelEnv

USE VectorTypes



IMPLICIT NONE



PRIVATE



!List of public members

PUBLIC :: SolutionAccelerationType

PUBLIC :: RelaxedPicardType

PUBLIC :: ModifiedPicardType



!> @brief the base Solution Acceleration Type

TYPE,ABSTRACT :: SolutionAccelerationType

  !> Initialization status

  LOGICAL(SBK) :: isInit=.FALSE.

  !> Size of solution vector/matrix

  INTEGER(SIK) :: n=-1

  !> Current iteration count

  INTEGER(SIK) :: s=0

  !> Starting iteration

  INTEGER(SIK) :: start=1

  !> Initial iterates

  CLASS(VectorType),ALLOCATABLE :: x(:)

  !> Futility computing environment

  TYPE(FutilityComputingEnvironment),POINTER :: ce => NULL()

!

!List of Type Bound Procedures

  CONTAINS

    !> @copybrief SolutionAccelerationType::init

    !> @copydetails SolutionAccelerationType::init

    PROCEDURE(init_absintfc),DEFERRED,PASS :: init

    !> @copybrief SolutionAccelerationType::init_base_SolutionAccelerationType

    !> @copydetails SolutionAccelerationType::init_base_SolutionAccelerationType

    PROCEDURE,PASS :: init_base => init_base_SolutionAccelerationType

    !> @copybrief SolutionAccelerationType::clear_SolutionAccelerationType

    !> @copydetails SolutionAccelerationType::clear_SolutionAccelerationType

    PROCEDURE,PASS :: clear => clear_base_SolutionAccelerationType

    PROCEDURE,PASS :: clear_base => clear_base_SolutionAccelerationType

    !> @copybrief SolutionAccelerationType::setInitial

    !> @copydetails SolutionAccelerationType::setInitial

    PROCEDURE(set_absintfc),DEFERRED,PASS :: setInitial

    !> @copybrief SolutionAccelerationType::step

    !> @copydetails SolutionAccelerationType::step

    PROCEDURE(step_absintfc),DEFERRED,PASS :: step

    !> @copybrief SolutionAccelerationType::reset

    !> @copydetails SolutionAccelerationType::reset

    PROCEDURE(set_absintfc),DEFERRED,PASS :: reset

ENDTYPE SolutionAccelerationType



!> Definition of interface for init

ABSTRACT INTERFACE

  SUBROUTINE init_absintfc(solver,ce,Params)

    IMPORT :: SolutionAccelerationType,FutilityComputingEnvironment,ParamType

    CLASS(SolutionAccelerationType),INTENT(INOUT) :: solver

    TYPE(FutilityComputingEnvironment),TARGET,INTENT(IN) :: ce

    TYPE(ParamType),INTENT(IN) :: Params

  ENDSUBROUTINE init_absintfc

ENDINTERFACE



!> Definition of interface for step

ABSTRACT INTERFACE

  SUBROUTINE step_absintfc(solver,x_new)

    IMPORT :: SolutionAccelerationType,VectorType

    CLASS(SolutionAccelerationType),INTENT(INOUT) :: solver

    CLASS(VectorType),INTENT(INOUT) :: x_new

  ENDSUBROUTINE step_absintfc

ENDINTERFACE



!> Definition of interface for setInitial and reset

ABSTRACT INTERFACE

  SUBROUTINE set_absintfc(solver,x)

    IMPORT :: SolutionAccelerationType,VectorType

    CLASS(SolutionAccelerationType),INTENT(INOUT) :: solver

    CLASS(VectorType),INTENT(INOUT) :: x

  ENDSUBROUTINE set_absintfc

ENDINTERFACE



!> @brief the Relaxed Picard Acceleration Type

TYPE,EXTENDS(SolutionAccelerationType) :: RelaxedPicardType

  !> Relaxation coefficient

  REAL(SRK) :: alpha=1.0_SRK

!

!List of Type Bound Procedures

  CONTAINS

    !> @copybrief RelaxedPicardType::init_RelaxedPicardType

    !> @copydetails RelaxedPicardType::init_RelaxedPicardType

    PROCEDURE,PASS :: init => init_RelaxedPicardType

    !> @copybrief RelaxedPicardType::clear_RelaxedPicardType

    !> @copydetails RelaxedPicardType::clear_RelaxedPicardType

    PROCEDURE,PASS :: clear => clear_RelaxedPicardType

    PROCEDURE,PASS :: clear_picardBase => clear_RelaxedPicardType

    !> @copybrief RelaxedPicardType::setInitial_RelaxedPicardType

    !> @copydetails RelaxedPicardType::setInitial_RelaxedPicardType

    PROCEDURE,PASS :: setInitial => setInitial_RelaxedPicardType

    PROCEDURE,PASS :: setInitial_picardBase => setInitial_RelaxedPicardType

    !> @copybrief RelaxedPicardType::step_RelaxedPicardType

    !> @copydetails RelaxedPicardType::step_RelaxedPicardType

    PROCEDURE,PASS :: step => step_RelaxedPicardType

    !> @copybrief RelaxedPicardType::reset_RelaxedPicardType

    !> @copydetails RelaxedPicardType::reset_RelaxedPicardType

    PROCEDURE,PASS :: reset => reset_RelaxedPicardType

ENDTYPE RelaxedPicardType



!> @brief the Relaxed Picard Acceleration Type

TYPE,EXTENDS(RelaxedPicardType) :: ModifiedPicardType

  !> Temporary vector to hold previous solution

  CLASS(VectorType),ALLOCATABLE :: tmpvec

!

!List of Type Bound Procedures

  CONTAINS

    !> @copybrief ModifiedPicardType::clear_ModifiedPicardType

    !> @copydetails ModifiedPicardType::clear_ModifiedPicardType

    PROCEDURE,PASS :: clear => clear_ModifiedPicardType

    !> @copybrief ModifiedPicardType::setInitial_ModifiedPicardType

    !> @copydetails ModifiedPicardType::setInitial_ModifiedPicardType

    PROCEDURE,PASS :: setInitial => setInitial_ModifiedPicardType

    !> @copybrief ModifiedPicardType::step_ModifiedPicardType

    !> @copydetails ModifiedPicardType::step_ModifiedPicardType

    PROCEDURE,PASS :: step => step_ModifiedPicardType

ENDTYPE ModifiedPicardType



!> Name of module

CHARACTER(LEN=*),PARAMETER :: modName='SolutionAccelerationModule'

!

!===============================================================================

CONTAINS

!

!-------------------------------------------------------------------------------

!> @brief Initializes the Solution Acceleration base class

!> @param solver The Solution Acceleration solver to act on

!> @param ce The computing environment to use for the calculation

!> @param Params The options parameter list

!>

SUBROUTINE init_base_SolutionAccelerationType(solver,ce,Params)

  CLASS(SolutionAccelerationType),INTENT(INOUT) :: solver

  TYPE(FutilityComputingEnvironment),TARGET,INTENT(IN) :: ce

  TYPE(ParamType),INTENT(IN) :: Params



  CHARACTER(LEN=*),PARAMETER :: myName='init_base_SolutionAccelerationType'



  REQUIRE(Params%has('SolutionAcceleration->num_unknowns'))

  REQUIRE(.NOT.solver%isInit)



  !Pull Data from Parameter List

  CALL Params%get('SolutionAcceleration->num_unknowns',solver%n)

  IF(Params%has('SolutionAcceleration->starting_iteration')) &

      CALL Params%get('SolutionAcceleration->starting_iteration',solver%start)



  IF(solver%n < 1) CALL ce%exceptHandler%raiseError('Incorrect input to '//modName// &

      '::'//myName//' - Number of unkowns (n) must be greater than 0!')



  IF(solver%start < 1) CALL ce%exceptHandler%raiseError('Incorrect input to '//modName// &

      '::'//myName//' - Starting iteration must be greater than 0!')



  solver%s=0

  solver%ce => ce

ENDSUBROUTINE init_base_SolutionAccelerationType

!

!-------------------------------------------------------------------------------

!> @brief Clears the solution acceleration base class

!> @param solver The solver to act on

!>

SUBROUTINE clear_base_SolutionAccelerationType(solver)

    CLASS(SolutionAccelerationType),INTENT(INOUT) :: solver



    INTEGER(SIK) :: i



    IF(ALLOCATED(solver%x)) THEN

      DO i=1,SIZE(solver%x)

        CALL solver%x(i)%clear()

      ENDDO

      DEALLOCATE(solver%x)

    ENDIF

    solver%s=0

    solver%n=-1

    solver%start=1

    solver%isInit=.FALSE.

  ENDSUBROUTINE clear_base_SolutionAccelerationType

!

!-------------------------------------------------------------------------------

!> @brief Initializes the relaxed picard type

!> @param solver The relaxed picard solver to act on

!> @param ce The computing environment to use for the calculation

!> @param Params The options parameter list

!>

SUBROUTINE init_RelaxedPicardType(solver,ce,Params)

  CLASS(RelaxedPicardType),INTENT(INOUT) :: solver

  TYPE(FutilityComputingEnvironment),TARGET,INTENT(IN) :: ce

  TYPE(ParamType),INTENT(IN) :: Params



  CHARACTER(LEN=*),PARAMETER :: myName='init_RelaxedPicardType'



  REQUIRE(.NOT.solver%isInit)



  CALL solver%init_base(ce,Params)



  IF(Params%has('SolutionAcceleration->picard_relaxation_factor')) &

      CALL Params%get('SolutionAcceleration->picard_relaxation_factor',solver%alpha)



  IF(solver%alpha <= 0.0_SRK) CALL ce%exceptHandler%raiseError('Incorrect input to '//modName// &

      '::'//myName//' - Relaxation coefficient must be greater than 0!')



  solver%isInit=.TRUE.

ENDSUBROUTINE init_RelaxedPicardType

!

!-------------------------------------------------------------------------------

!> @brief Clears the relaxed picard type

!> @param solver The solver to act on

!>

SUBROUTINE clear_RelaxedPicardType(solver)

  CLASS(RelaxedPicardType),INTENT(INOUT) :: solver



  CALL solver%clear_base()



  solver%alpha=1.0_SRK

  solver%isInit=.FALSE.



ENDSUBROUTINE clear_RelaxedPicardType

!

!-------------------------------------------------------------------------------

!> @brief Sets the initial vector for relaxation

!> @param solver Solver to set initial vector for

!> @param x  Vector to use as initial condition and all other vector types have

!>           to match size and type of this one

!>

!> This sets the initial vector into internal storage.  Also initializes solver%x

!> based on the type that is passed in.

!>

SUBROUTINE setInitial_RelaxedPicardType(solver,x)

  CLASS(RelaxedPicardType),INTENT(INOUT) :: solver

  CLASS(VectorType),INTENT(INOUT) :: x



#if defined(HAVE_MPI) || defined(FUTILITY_HAVE_Trilinos)

  CHARACTER(LEN=*),PARAMETER :: myName='setInitial_RelaxedPicardType'

#endif



  REQUIRE(x%n == solver%n)

  REQUIRE(solver%s == 0)



  !If this is the first call to set/reset must actually create vectors of needed type

  IF(.NOT.ALLOCATED(solver%x)) THEN

    !!!TODO Once missing BLAS interfaces have been added for the following vector types

    !do away with this error catch to allow use of these



    SELECT TYPE(x)

#ifdef HAVE_MPI

    TYPE IS(NativeDistributedVectorType)

      CALL solver%ce%exceptHandler%raiseError('Incorrect call to '//modName// &

          '::'//myName//' - Input vector type not supported!')

#endif

#ifdef FUTILITY_HAVE_Trilinos

    TYPE IS(TrilinosVectorType)

      CALL solver%ce%exceptHandler%raiseError('Incorrect call to '//modName// &

          '::'//myName//' - Input vector type not supported!')

#endif

    ENDSELECT



    CALL VectorResembleAlloc(solver%x,x,1)

  ENDIF



  CALL BLAS_copy(x,solver%x(1))



ENDSUBROUTINE setInitial_RelaxedPicardType

!

!-------------------------------------------------------------------------------

!> @brief Performs a single step of Relaxed Picard acting on the input solution vector.

!> @param solver Solver to take step with

!> @param x_new New solution iterate and under-relaxed / accelerated return vector

!>

SUBROUTINE step_RelaxedPicardType(solver,x_new)

  CLASS(RelaxedPicardType),INTENT(INOUT) :: solver

  CLASS(VectorType),INTENT(INOUT) :: x_new



  REQUIRE(x_new%n == solver%n)

  REQUIRE(ALLOCATED(solver%x))



  !Update iteration counter

  solver%s=solver%s+1



  IF(solver%s >= solver%start) THEN

    !Get under-relaxed solution using mixing parameter

    CALL BLAS_scal(x_new,solver%alpha)

    CALL BLAS_axpy(solver%x(1),x_new,1.0_SRK-solver%alpha)

  ENDIF

  CALL BLAS_copy(x_new,solver%x(1))



ENDSUBROUTINE step_RelaxedPicardType

!

!-------------------------------------------------------------------------------

!> @brief Reset the initial iterate for the solver

!> @param solver The solver to act on

!> @param x Initial iterate solve is starting from

!>

SUBROUTINE reset_RelaxedPicardType(solver,x)

  CLASS(RelaxedPicardType),INTENT(INOUT) :: solver

  CLASS(VectorType),INTENT(INOUT) :: x



  REQUIRE(x%n == solver%n)



  !Reset iteration counter

  solver%s=0_SIK



  CALL solver%setInitial(x)



ENDSUBROUTINE reset_RelaxedPicardType

!

!-------------------------------------------------------------------------------

!> @brief Clears the modified picard Type

!> @param solver The solver to act on

!>

SUBROUTINE clear_ModifiedPicardType(solver)

  CLASS(ModifiedPicardType),INTENT(INOUT) :: solver



  CALL solver%clear_picardbase()



  IF(ALLOCATED(solver%tmpvec)) THEN

    CALL solver%tmpvec%clear()

    DEALLOCATE(solver%tmpvec)

  ENDIF

  solver%isInit=.FALSE.

ENDSUBROUTINE clear_ModifiedPicardType

!

!-------------------------------------------------------------------------------

!> @brief Sets the initial vector for relaxation

!> @param solver Solver to set initial vector for

!> @param x  Vector to use as initial condition and all other vector types have

!>           to match size and type of this one

!>

!> This sets the initial vector into internal storage.  Also initializes solver%x

!> and tmpvec based on the type that is passed in.

!>

SUBROUTINE setInitial_ModifiedPicardType(solver,x)

  CLASS(ModifiedPicardType),INTENT(INOUT) :: solver

  CLASS(VectorType),INTENT(INOUT) :: x



  CALL solver%setInitial_picardBase(x)



  IF(.NOT. ALLOCATED(solver%tmpvec)) CALL VectorResembleAlloc(solver%tmpvec,x)



ENDSUBROUTINE setInitial_ModifiedPicardType

!

!-------------------------------------------------------------------------------

!> @brief Performs a single step of Modified Picard acting on the input solution vector.

!> @param solver Solver to take step with

!> @param x_new New solution iterate and under-relaxed / accelerated return vector

!>

SUBROUTINE step_ModifiedPicardType(solver,x_new)

  CLASS(ModifiedPicardType),INTENT(INOUT) :: solver

  CLASS(VectorType),INTENT(INOUT) :: x_new



  REQUIRE(x_new%n == solver%n)

  REQUIRE(ALLOCATED(solver%x))

  REQUIRE(ALLOCATED(solver%tmpvec))



  !Update iteration counter

  solver%s=solver%s+1



  CALL BLAS_copy(x_new,solver%tmpvec)



  IF(solver%s >= solver%start) THEN

    !Get under-relaxed solution using mixing parameter

    CALL BLAS_scal(x_new,solver%alpha)

    CALL BLAS_axpy(solver%x(1),x_new,1.0_SRK-solver%alpha)

  ENDIF

  CALL BLAS_copy(solver%tmpvec,solver%x(1))



ENDSUBROUTINE step_ModifiedPicardType

!

ENDMODULE SolutionAccelerationModule

unit_tests/CMakeLists.txt

+1 −0
Original line number Diff line number Diff line
@@ -66,6 +66,7 @@ SET(UNIT_TEST_NAMES 
    testEigenvalueSolver

    testODESolver

    testNonLinearSolver

    testPicard

    testAndersonAcceleration

    testStochasticSampler

    testPartitionGraph

unit_tests/testAndersonAcceleration/testAndersonAcceleration.f90

+13 −13

File changed.

Preview size limit exceeded, changes collapsed.