Using special functions (bessel, gamma, etc) in parsed expressions

[richmondodufisan]

Check these boxes if you have followed the posting rules.

• Q&A General is the most appropriate section for my question
• I have consulted the posting Guidelines on the Discussions front page
• I have searched the Discussions forum and my question has not been asked before
• I have searched the MOOSE website and the documentation does not answer my question
• I have formatted my post following the posting guidelines (screenshots as a last resort, triple back quotes around pasted text)

Question

Hi,

I was wondering if there was a way to use special functions in MOOSE. I looked at the function parser documentation:

http://warp.povusers.org/FunctionParser/fparser.html#literals

and saw that it only really covers basic ones like sin/tan/etc. I was wondering if there is a way to implement gamma, bessel functions, etc without needing to use an approximation.

[GiudGiud]

Hello

I don't think they can be used in parsed expressions.
I think the solution here is to code your own object instead of relying on the parsed capabilities

[richmondodufisan]

if I go this route, say, create an object file. Do I just need to include cmath and then I can use bessel functions, etc? And then also use _q_point to get the xyz quadrature points?

And then I'm guessing the material file would need to be block restricted (or surface restricted) since I'm using the function to define boundary conditions.

[richmondodufisan]

If I go this route, what is the default value for "h"? Is it 1.0? Or is h the "value" parameter with the default set to 0?

https://mooseframework.inl.gov/source/bcs/MatNeumannBC.html

That is, if I only set the required parameters, will it fix the flux as the value of the constant material property?

[richmondodufisan]

Does this option work, and/or is there a better way to do this than coding a material property?

[GiudGiud]

I'm not sure why I did not get notifications on this thread. Sorry about that.

. Do I just need to include cmath and then I can use bessel functions, etc? And then also use _q_point to get the xyz quadrature points?

yes that's all you need to do

And then I'm guessing the material file would need to be block restricted (or surface restricted) since I'm using the function to define boundary conditions.

If you need to define a BC, why are you creating a Material? I think you might want to create a BC object instead

If I go this route, what is the default value for "h"? Is it 1.0? Or is h the "value" parameter with the default set to 0?

hM is defined by the material. So effectively I d say h=1 and M is the property

Does this option work, and/or is there a better way to do this than coding a material property?

I would say the better way is to create a NeumannBC, somewhat similar to FunctionNeumannBC, but uses the hard coded function you need instead of a parameter

[richmondodufisan]

Okay I think I'm going to go this route, this makes sense. I should be able to do Dirichlet BCs too with this then. Thanks!

[dschwen]

The function parser is quite easy to extend. This would have to be done in libmesh to get the derivatives implemented as well. If there is interest I could add bessel functions etc. However the standard library only supports a subset of those functions. What exactly do you need?

[richmondodufisan]

Okay I found it, the search function didn't find it for some reason but it was in the directory you shared. I added the functions in the source file:

////////////////////////////////// Added by Richmond Odufisan //////////////////////////////////////

template <>
std::complex<double> FunctionParserADBase<std::complex<double>>::fp_i0(const std::complex<double>* params) 
{
    double order = 0.0;
    double argument = params[0].real(); // Extract the real part
    return std::complex<double>(std::cyl_bessel_i(order, argument), 0.0); // Result is real
}

template <>
std::complex<double> FunctionParserADBase<std::complex<double>>::fp_j0(const std::complex<double>* params) 
{
    double order = 0.0;
    double argument = params[0].real(); // Extract the real part
    return std::complex<double>(std::cyl_bessel_j(order, argument), 0.0); // Result is real
}

template <>
std::complex<double> FunctionParserADBase<std::complex<double>>::fp_gamma(const std::complex<double>* params) 
{
    double argument = params[0].real(); // Extract the real part
    return std::complex<double>(std::tgamma(argument), 0.0); // Result is real
}

template <>
std::complex<double> FunctionParserADBase<std::complex<double>>::fp_erfc(const std::complex<double>* params) 
{
    double argument = params[0].real(); // Extract the real part
    return std::complex<double>(std::erfc(argument), 0.0); // Result is real
}



template <>
std::complex<long double> FunctionParserADBase<std::complex<long double>>::fp_i0(const std::complex<long double>* params) 
{
    long double order = 0.0L;
    long double argument = params[0].real(); // Extract the real part
    return std::complex<long double>(std::cyl_bessel_i(order, argument), 0.0L); // Result is real
}

template <>
std::complex<long double> FunctionParserADBase<std::complex<long double>>::fp_j0(const std::complex<long double>* params) 
{
    long double order = 0.0L;
    long double argument = params[0].real(); // Extract the real part
    return std::complex<long double>(std::cyl_bessel_j(order, argument), 0.0L); // Result is real
}

template <>
std::complex<long double> FunctionParserADBase<std::complex<long double>>::fp_gamma(const std::complex<long double>* params) 
{
    long double argument = params[0].real(); // Extract the real part
    return std::complex<long double>(std::tgamma(argument), 0.0L); // Result is real
}

template <>
std::complex<long double> FunctionParserADBase<std::complex<long double>>::fp_erfc(const std::complex<long double>* params) 
{
    long double argument = params[0].real(); // Extract the real part
    return std::complex<long double>(std::erfc(argument), 0.0L); // Result is real
}


////////////////////////////////////////////////////////////////////////////////////////////////


template<typename Value_t>
FunctionParserADBase<Value_t>::FunctionParserADBase() :
    FunctionParserBase<Value_t>(),
    compiledFunction(NULL),
    mFPlog(this->mData->mFuncPtrs.size()),
    mADFlags(ADJITCache),
    mRegisteredDerivatives(),
    ad(new ADImplementation<Value_t>(this))
{
  this->AddFunction("plog", fp_plog, 2);
  mFErf = this->mData->mFuncPtrs.size();
  this->AddFunction("erf", fp_erf, 1);
  
  //////////////////// Added by Richmond Odufisan /////////////////////////////////
  
  this->AddFunction("i0", fp_i0, 1);
  this->AddFunction("j0", fp_j0, 1);
  this->AddFunction("gamma", fp_gamma, 1);
  this->AddFunction("erfc", fp_erfc, 1);
  
  /////////////////////////////////////////////////////////////////////////////////
}

And built libmesh with the command you shared, and got this error at the end:

---------------------------------------------
-- Done configuring core library features ---
---------------------------------------------
configure: error: You must run "git submodule update --init --recursive" before configuring libmesh
Running make -j 1...
make: *** No targets specified and no makefile found.  Stop.

Is it okay to continue even with this error? Since the command specifically tried to skip submodule updates. The last step to follow would be to do:

cd ~/projects/moose/test
make -j 6

to build MOOSE, and then to create my app and try to run something with the newly added functions

[GiudGiud]

go in moose/libmesh and run git submodule update --init --recursive

then try that script again. --skip-submodule-update only works well if you got all the submodules properly first :)

[richmondodufisan]

okay I got it to work now. built libmesh (no errors), was asked to build wasp, and then tried to compile moose. but then I got this error, clearly caused by the functions I added. Did I miss a step somewhere? What does this mean?

/home/richmond98/mambaforge3/envs/moose/bin/../lib/gcc/x86_64-conda-linux-gnu/12.3.0/../../../../x86_64-conda-linux-gnu/bin/ld: /home/richmond98/projects/moose/scripts/../libmesh/installed/lib/libmesh_opt.so: undefined reference to `FunctionParserADBase<double>::fp_j0(double const*)'
/home/richmond98/mambaforge3/envs/moose/bin/../lib/gcc/x86_64-conda-linux-gnu/12.3.0/../../../../x86_64-conda-linux-gnu/bin/ld: /home/richmond98/projects/moose/scripts/../libmesh/installed/lib/libmesh_opt.so: undefined reference to `FunctionParserADBase<float>::fp_gamma(float const*)'
/home/richmond98/mambaforge3/envs/moose/bin/../lib/gcc/x86_64-conda-linux-gnu/12.3.0/../../../../x86_64-conda-linux-gnu/bin/ld: /home/richmond98/projects/moose/scripts/../libmesh/installed/lib/libmesh_opt.so: undefined reference to `FunctionParserADBase<long double>::fp_gamma(long double const*)'
/home/richmond98/mambaforge3/envs/moose/bin/../lib/gcc/x86_64-conda-linux-gnu/12.3.0/../../../../x86_64-conda-linux-gnu/bin/ld: /home/richmond98/projects/moose/scripts/../libmesh/installed/lib/libmesh_opt.so: undefined reference to `FunctionParserADBase<long double>::fp_erfc(long double const*)'
/home/richmond98/mambaforge3/envs/moose/bin/../lib/gcc/x86_64-conda-linux-gnu/12.3.0/../../../../x86_64-conda-linux-gnu/bin/ld: /home/richmond98/projects/moose/scripts/../libmesh/installed/lib/libmesh_opt.so: undefined reference to `FunctionParserADBase<float>::fp_erfc(float const*)'
/home/richmond98/mambaforge3/envs/moose/bin/../lib/gcc/x86_64-conda-linux-gnu/12.3.0/../../../../x86_64-conda-linux-gnu/bin/ld: /home/richmond98/projects/moose/scripts/../libmesh/installed/lib/libmesh_opt.so: undefined reference to `FunctionParserADBase<std::complex<float> >::fp_j0(std::complex<float> const*)'
/home/richmond98/mambaforge3/envs/moose/bin/../lib/gcc/x86_64-conda-linux-gnu/12.3.0/../../../../x86_64-conda-linux-gnu/bin/ld: /home/richmond98/projects/moose/scripts/../libmesh/installed/lib/libmesh_opt.so: undefined reference to `FunctionParserADBase<float>::fp_j0(float const*)'
/home/richmond98/mambaforge3/envs/moose/bin/../lib/gcc/x86_64-conda-linux-gnu/12.3.0/../../../../x86_64-conda-linux-gnu/bin/ld: /home/richmond98/projects/moose/scripts/../libmesh/installed/lib/libmesh_opt.so: undefined reference to `FunctionParserADBase<long double>::fp_j0(long double const*)'
/home/richmond98/mambaforge3/envs/moose/bin/../lib/gcc/x86_64-conda-linux-gnu/12.3.0/../../../../x86_64-conda-linux-gnu/bin/ld: /home/richmond98/projects/moose/scripts/../libmesh/installed/lib/libmesh_opt.so: undefined reference to `FunctionParserADBase<double>::fp_i0(double const*)'
/home/richmond98/mambaforge3/envs/moose/bin/../lib/gcc/x86_64-conda-linux-gnu/12.3.0/../../../../x86_64-conda-linux-gnu/bin/ld: /home/richmond98/projects/moose/scripts/../libmesh/installed/lib/libmesh_opt.so: undefined reference to `FunctionParserADBase<double>::fp_gamma(double const*)'
/home/richmond98/mambaforge3/envs/moose/bin/../lib/gcc/x86_64-conda-linux-gnu/12.3.0/../../../../x86_64-conda-linux-gnu/bin/ld: /home/richmond98/projects/moose/scripts/../libmesh/installed/lib/libmesh_opt.so: undefined reference to `FunctionParserADBase<std::complex<float> >::fp_erfc(std::complex<float> const*)'
/home/richmond98/mambaforge3/envs/moose/bin/../lib/gcc/x86_64-conda-linux-gnu/12.3.0/../../../../x86_64-conda-linux-gnu/bin/ld: /home/richmond98/projects/moose/scripts/../libmesh/installed/lib/libmesh_opt.so: undefined reference to `FunctionParserADBase<float>::fp_i0(float const*)'
/home/richmond98/mambaforge3/envs/moose/bin/../lib/gcc/x86_64-conda-linux-gnu/12.3.0/../../../../x86_64-conda-linux-gnu/bin/ld: /home/richmond98/projects/moose/scripts/../libmesh/installed/lib/libmesh_opt.so: undefined reference to `FunctionParserADBase<std::complex<float> >::fp_gamma(std::complex<float> const*)'
/home/richmond98/mambaforge3/envs/moose/bin/../lib/gcc/x86_64-conda-linux-gnu/12.3.0/../../../../x86_64-conda-linux-gnu/bin/ld: /home/richmond98/projects/moose/scripts/../libmesh/installed/lib/libmesh_opt.so: undefined reference to `FunctionParserADBase<double>::fp_erfc(double const*)'
/home/richmond98/mambaforge3/envs/moose/bin/../lib/gcc/x86_64-conda-linux-gnu/12.3.0/../../../../x86_64-conda-linux-gnu/bin/ld: /home/richmond98/projects/moose/scripts/../libmesh/installed/lib/libmesh_opt.so: undefined reference to `FunctionParserADBase<std::complex<float> >::fp_i0(std::complex<float> const*)'
/home/richmond98/mambaforge3/envs/moose/bin/../lib/gcc/x86_64-conda-linux-gnu/12.3.0/../../../../x86_64-conda-linux-gnu/bin/ld: /home/richmond98/projects/moose/scripts/../libmesh/installed/lib/libmesh_opt.so: undefined reference to `FunctionParserADBase<long double>::fp_i0(long double const*)'
collect2: error: ld returned 1 exit status
make: *** [/home/richmond98/projects/moose/framework/moose.mk:475: /home/richmond98/projects/moose/framework/contrib/exodiff/exodiff] Error 1
make: *** Waiting for unfinished jobs....

I don't know if it has anything to do with the header file:

template<typename Value_t>
class FunctionParserADBase : public FunctionParserBase<Value_t>
{
public:
  FunctionParserADBase();
  FunctionParserADBase(const FunctionParserADBase& cpy);
  virtual ~FunctionParserADBase();

  /**
   * This class manages its own memory, so the compiler-generated copy
   * assignment, move assignment, and move constructor implementations
   * are not safe to use.  We therefore explicitly delete them so they
   * can't be called accidentally.
   */
  FunctionParserADBase (FunctionParserADBase &&) = delete;
  FunctionParserADBase & operator= (const FunctionParserADBase &) = delete;
  FunctionParserADBase & operator= (FunctionParserADBase &&) = delete;

  ...

  // user function plog
  static Value_t fp_plog(const Value_t * params);

  // user function erf
  static Value_t fp_erf(const Value_t * params);
  
  /////////////////// Added by Richmond Odufisan ///////////////////////////
  
  static Value_t fp_i0(const Value_t* params);
  
  static Value_t fp_j0(const Value_t* params);
  
  static Value_t fp_gamma(const Value_t* params);
  
  static Value_t fp_erfc(const Value_t* params);

/////////////////////////////////////////////////////////////////////////

  // function ID for the plog function
  unsigned int mFPlog;

  // function ID for the erf function
  unsigned int mFErf;

};

Where there is a function ID for the plog and erf functions. Looking in the source code for what those "function ID" attributes do does not look like a simple fix/implementation.

[GiudGiud]

did you define these functions in a source file?

It looks like you did some long double and doubles. Did you implement everything that was declared?
Ideally there would be a default template implementation, and you only use template specification for the few types that require it

[richmondodufisan]

Yes I did, I tried several different ways to get it to compile. The existing example in the source file had three definitions:

template <typename Value_t>
Value_t FunctionParserADBase<Value_t>::fp_erf(const Value_t * params)
{
  return Faddeeva::erf(params[0]);
}

template<>
std::complex<float>
FunctionParserADBase<std::complex<float> >::fp_erf(const std::complex<float> * params)
{
  std::complex<double> result =
    Faddeeva::erf(std::complex<double>(params[0].real(),
                                       params[0].imag()));
  return std::complex<float>(result.real(), result.imag());
}

template<>
std::complex<long double>
FunctionParserADBase<std::complex<long double> >::fp_erf(const std::complex<long double> * params)
{
  std::complex<double> result =
    Faddeeva::erf(std::complex<double>(params[0].real(),
                                       params[0].imag()));
  return std::complex<long double>(result.real(), result.imag());
}

And after a lot of tinkering I was able to get mine to compile similarly (with all three). I eventually ended up leaving out the last two (complex long doubles and doubles) when I got this to compile:

template <typename Value_t>
static Value_t fp_i0(const Value_t* params) {
    float argument;
    if constexpr (std::is_scalar<Value_t>::value) {
        argument = static_cast<float>(params[0]);
    } else {
        argument = static_cast<float>(params[0].real());
    }
    return static_cast<Value_t>(std::cyl_bessel_i(0.0, argument));
}

template <typename Value_t>
static Value_t fp_j0(const Value_t* params) {
    float argument;
    if constexpr (std::is_scalar<Value_t>::value) {
        argument = static_cast<float>(params[0]);
    } else {
        argument = static_cast<float>(params[0].real());
    }
    return static_cast<Value_t>(std::cyl_bessel_j(0.0, argument));
}

template <typename Value_t>
static Value_t fp_gamma(const Value_t* params) {
    float argument;
    if constexpr (std::is_scalar<Value_t>::value) {
        argument = static_cast<float>(params[0]);
    } else {
        argument = static_cast<float>(params[0].real());
    }
    return static_cast<Value_t>(std::tgamma(argument));
}

template <typename Value_t>
static Value_t fp_erfc(const Value_t* params) {
    float argument;
    if constexpr (std::is_scalar<Value_t>::value) {
        argument = static_cast<float>(params[0]);
    } else {
        argument = static_cast<float>(params[0].real());
    }
    return static_cast<Value_t>(std::erfc(argument));
}

But it still gave me the same error as before when compiling MOOSE.