From 4e37d37094260e0e6cecd537e44074600973d67d Mon Sep 17 00:00:00 2001 From: "Balos, Cody, J" Date: Tue, 18 Jun 2024 11:35:33 -0700 Subject: [PATCH] fix fixedpoint sunnonlinear solver fortran test --- .../test_fsunnonlinsol_fixedpoint_mod.f90 | 155 +++++++++++------- 1 file changed, 99 insertions(+), 56 deletions(-) diff --git a/examples/sunnonlinsol/fixedpoint/test_fsunnonlinsol_fixedpoint_mod.f90 b/examples/sunnonlinsol/fixedpoint/test_fsunnonlinsol_fixedpoint_mod.f90 index f084b816bf..7b6cfd2322 100644 --- a/examples/sunnonlinsol/fixedpoint/test_fsunnonlinsol_fixedpoint_mod.f90 +++ b/examples/sunnonlinsol/fixedpoint/test_fsunnonlinsol_fixedpoint_mod.f90 @@ -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 @@ -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 @@ -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) @@ -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 @@ -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 @@ -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)