Skip to content

Commit

Permalink
fix fixedpoint sunnonlinear solver fortran test
Browse files Browse the repository at this point in the history
  • Loading branch information
@balos1
balos1 committed Jun 18, 2024
1 parent 45243ac commit 4e37d37
Showing 1 changed file with 99 additions and 56 deletions.
155 changes: 99 additions & 56 deletions examples/sunnonlinsol/fixedpoint/test_fsunnonlinsol_fixedpoint_mod.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,15 +22,17 @@ module test_fsunnonlinsol_fixedpoint
implicit none

integer(kind=myindextype), parameter :: NEQ = 3 ! number of equations
integer(C_INT), parameter :: MAXIT = 50 ! max nonlinear iters.
integer(C_INT), parameter :: MAXIT = 20 ! max nonlinear iters.
real(C_DOUBLE), parameter :: TOL = 1.0e-4 ! nonlinear solver tolerance

real(C_DOUBLE), parameter :: PI = 3.1415926535898

! approximate solution
real(C_DOUBLE) :: Y1 = 0.5d0
real(C_DOUBLE) :: Y2 = 0.d0
real(C_DOUBLE) :: Y3 = -PI/6.0d0
real(C_DOUBLE) :: XTRUE = 0.5d0
real(C_DOUBLE) :: YTRUE = 1.0d0
real(C_DOUBLE) :: ZTRUE = -PI/6.0d0

type(N_Vector), pointer :: y0

contains

Expand All @@ -43,26 +45,30 @@ integer(C_INT) function unit_tests() result(retval)
implicit none

type(SUNNonlinearSolver), pointer :: NLS ! test nonlinear solver
type(N_Vector), pointer :: x, y0, y, w ! test vectors
real(C_DOUBLE), pointer :: ydata(:)
type(N_Vector), pointer :: ycur, ycor, w ! test vectors
real(C_DOUBLE), pointer :: data(:)
integer(C_LONG) :: niters(1)
integer(C_INT) :: tmp

x => FN_VNew_Serial(NEQ, sunctx)
y0 => FN_VClone(x)
y => FN_VClone(x)
w => FN_VClone(x)
y0 => FN_VNew_Serial(NEQ, sunctx)
ycor => FN_VClone(y0)
ycur => FN_VClone(y0)
w => FN_VClone(y0)

! set weights
ydata => FN_VGetArrayPointer(y0)
ydata(1) = 0.1d0
ydata(2) = 0.1d0
ydata(3) = -0.1d0
! set initial guess
data => FN_VGetArrayPointer(y0)
data(1) = 0.1d0
data(2) = 0.1d0
data(3) = -0.1d0

! set initial correction
call FN_VConst(0.0d0, ycor)

! set weights
call FN_VConst(1.0d0, w)

! create and test NLS
NLS => FSUNNonlinsol_FixedPoint(y, 0, sunctx)
NLS => FSUNNonlinsol_FixedPoint(y0, 0, sunctx)

retval = FSUNNonlinSolSetSysFn(NLS, c_funloc(FPFunction))
if (retval /= 0) then
Expand All @@ -82,37 +88,35 @@ integer(C_INT) function unit_tests() result(retval)
return
end if

retval = FSUNNonlinSolSolve(NLS, y0, y, w, TOL, 1, c_loc(x))
retval = FSUNNonlinSolSolve(NLS, y0, ycor, w, TOL, 1, c_loc(y0))
if (retval /= 0) then
write(*,'(A,I0)') ' >>> FAIL: FSUNNonlinSolSolve returned ', retval
return
end if

! extract and print solution
ydata => FN_VGetArrayPointer(y)

write(*,*) 'Solution:'
write(*,'(A,E14.7)') 'y1 = ', ydata(1)
write(*,'(A,E14.7)') 'y2 = ', ydata(2)
write(*,'(A,E14.7)') 'y3 = ', ydata(3)

write(*,*) 'Solution Error:'
write(*,'(A,E14.7)') 'e1 = ', ydata(1) - Y1
write(*,'(A,E14.7)') 'e2 = ', ydata(2) - Y2
write(*,'(A,E14.7)') 'e3 = ', ydata(3) - Y3
! update the initial guess with the final correction
call FN_VLinearSum(1.0d0, y0, 1.0d0, ycor, ycur);

! print number of iterations
retval = FSUNNonlinSolGetNumIters(NLS, niters)
if (retval /= 0) then
write(*,'(A,I0)') ' >>> FAIL: FSUNNonlinSolGetNumIters returned ', retval
return
end if

write(*,'(A,I0)') 'Number of nonlinear iterations:', niters(1)
write(*,'(A,I0)') 'Number of nonlinear iterations: ', niters(1)

! check answer
retval = check_ans(ycur, TOL)
if (retval /= 0) then
write(*,'(A,I0)') ' >>> FAIL: check_ans failed'
return
end if

! cleanup
call FN_VDestroy(x)
call FN_VDestroy(y0)
call FN_VDestroy(y)
call FN_VDestroy(ycor)
call FN_VDestroy(ycur)
call FN_VDestroy(w)
tmp = FSUNNonlinSolFree(NLS)

Expand All @@ -122,8 +126,6 @@ integer(C_INT) function ConvTest(NLS, y, del, tol, ewt, mem) &
result(retval) bind(C)
use, intrinsic :: iso_c_binding



implicit none

type(SUNNonlinearSolver) :: NLS
Expand All @@ -144,46 +146,87 @@ integer(C_INT) function ConvTest(NLS, y, del, tol, ewt, mem) &
end function

! -----------------------------------------------------------------------------
! Nonlinear system
! Nonlinear system F(x,y,z):
!
! 3x - cos(yz) - 1/2 = 0
! x^2 - 81(y+0.1)^2 + sin(z) + 1.06 = 0
! exp(-xy) + 20z + (10 pi - 3)/3 = 0
! 3x - cos((y-1)z) - 1/2 = 0
! x^2 - 81(y-0.9)^2 + sin(z) + 1.06 = 0
! exp(-x(y-1)) + 20z + (10 pi - 3)/3 = 0
!
! Nonlinear fixed point function
! Nonlinear fixed point function G(x,y,z):
!
! g1(x,y,z) = 1/3 cos(yz) + 1/6
! g2(x,y,z) = 1/9 sqrt(x^2 + sin(z) + 1.06) - 0.1
! g3(x,y,z) = -1/20 exp(-xy) - (10 pi - 3) / 60
! G1(x,y,z) = 1/3 cos((y-1)yz) + 1/6
! G2(x,y,z) = 1/9 sqrt(x^2 + sin(z) + 1.06) + 0.9
! G3(x,y,z) = -1/20 exp(-x(y-1)) - (10 pi - 3) / 60
!
! Corrector form g(x,y,z):
!
! g1(x,y,z) = 1/3 cos((y-1)yz) + 1/6 - x0
! g2(x,y,z) = 1/9 sqrt(x^2 + sin(z) + 1.06) + 0.9 - y0
! g3(x,y,z) = -1/20 exp(-x(y-1)) - (10 pi - 3) / 60 - z0
!
! ----------------------------------------------------------------------------
integer(C_INT) function FPFunction(y, f, mem) &
integer(C_INT) function FPFunction(ycor, f, mem) &
result(retval) bind(C)
use, intrinsic :: iso_c_binding


implicit none

type(N_Vector) :: y, f
type(N_Vector) :: ycor, f
type(C_PTR), value :: mem
real(C_DOUBLE), pointer :: ydata(:), fdata(:)
real(C_DOUBLE) :: y1, y2, y3
real(C_DOUBLE), pointer :: data(:), fdata(:)
real(C_DOUBLE) :: x, y, z

ydata => FN_VGetArrayPointer(y)
data => FN_VGetArrayPointer(ycor)
fdata => FN_VGetArrayPointer(f)

y1 = ydata(1)
y2 = ydata(2)
y3 = ydata(3)
x = data(1)
y = data(2)
z = data(3)

fdata(1) = (1/3.0d0) * cos(y2*y3) + (1/6.0d0)
fdata(2) = (1/9.0d0) * sqrt(y1*y1 + sin(y3) + 1.06d0) - 0.1d0
fdata(3) = -(1/20.d0) * exp(-y1*y2) - (10.d0 * PI - 3.0d0) / 60.d0
fdata(1) = (1.0d0/3.0d0) * cos((y - 1.0d0) * z) + (1.0d0/6.0d0)
fdata(2) = (1.0d0/9.0d0) * sqrt(x*x + sin(z) + 1.06d0) + 0.9d0
fdata(3) = -(1/20.d0) * exp(-x*(y-1.0d0)) - (10.d0 * PI - 3.0d0) / 60.0d0

call FN_VLinearSum(1.0d0, f, -1.0d0, y0, f)

retval = 0

end function

integer(C_INT) function check_ans(ycor, tol) &
result(retval) bind(C)
use, intrinsic :: iso_c_binding
implicit none

type(N_Vector) :: ycor
real(C_DOUBLE), value :: tol
real(C_DOUBLE) :: ex, ey, ez
real(C_DOUBLE), pointer :: data(:)

! extract and print solution
data => FN_VGetArrayPointer(ycor)

write(*,*) 'Solution:'
write(*,'(A,E14.7)') ' x = ', data(1)
write(*,'(A,E14.7)') ' y = ', data(2)
write(*,'(A,E14.7)') ' z = ', data(3)

ex = data(1) - XTRUE
ey = data(2) - YTRUE
ez = data(3) - ZTRUE

write(*,*) 'Solution Error:'
write(*,'(A,E14.7)') ' ex = ', ex
write(*,'(A,E14.7)') ' ey = ', ey
write(*,'(A,E14.7)') ' ez = ', ez

tol = tol * 10.0d0
if (ex > tol .or. ey > tol .or. ez > tol) then
retval = 1
end if

retval = 0
end function

end module

program main
Expand All @@ -198,7 +241,7 @@ program main
integer(C_INT) :: fails = 0

!============== Introduction =============
print *, 'fixedpoint SUNNonlinearSolver Fortran 2003 interface test'
write(*,*) 'SUNNonlinearSolver_FixedPoint Fortran 2003 interface test'

call Test_Init(SUN_COMM_NULL)

Expand Down

0 comments on commit 4e37d37

Please sign in to comment.