test/dictionary/testCalcNewton: Example of #calc with #{...#} to find the root of a function

using Newton-Raphson iteration.
This commit is contained in:
Henry Weller
2023-07-18 15:26:46 +01:00
parent bebe584842
commit cc9991e9a3
3 changed files with 58 additions and 9 deletions

View File

@ -1,9 +0,0 @@
#!/bin/sh
cd ${0%/*} || exit 1 # Run from this directory
# Source tutorial run functions
. $WM_PROJECT_DIR/bin/tools/CleanFunctions
cleanCase
#------------------------------------------------------------------------------

View File

@ -6,6 +6,7 @@ cd ${0%/*} || exit 1 # Run from this directory
runApplication -s testDict foamDictionary -expand testDict
runApplication -s testCalc foamDictionary -expand testCalc
runApplication -s testCalcNewton foamDictionary -expand testCalcNewton
runApplication -s simpleBlockMeshDict foamDictionary -expand simpleBlockMeshDict
#------------------------------------------------------------------------------

View File

@ -0,0 +1,57 @@
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: dev
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
object testCalcNewton;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Find the root of f(x) = x + B*x/sqrt(1 + sqr(x)) - A
// using Newton-Raphson iteration
A 1;
B 1;
x #calc
#{
// Lookup the coefficients
const scalar A = $A;
const scalar B = $B;
// Initial guess for x
scalar x = 0;
scalar x0 = x;
do
{
// Store the previous iteration x for the convergence check
x0 = x;
// Temporary sub-function evaluations
const scalar f1 = 1 + sqr(x);
const scalar f2 = sqrt(f1);
// Evaluate the function
const scalar f = x + B*x/f2 - A;
// Evaluate the derivative
const scalar df = 1 + B/(f1*f2);
// Update x
x = x0 - f/df;
// Test for convergence
} while (mag(x - x0) > 1e-6);
os << x;
#};
// ************************************************************************* //