How to use modified Bessel function of 2nd kind from scipy within Expression?

Question

I tried this code

   from dolfin import *
   import numpy
   import os
   from scipy import pi, linspace, loadtxt, meshgrid, exp, cos, sqrt, power
   from scipy.special import kn

   # System size definition 
   # system = (x1, x1+Lx) x (y1, y1+Ly)
   Lx = 16.0; Ly = 16.0
   x1 = 0.01; y1 = 0.01

   # Create mesh and define function space
   dx0 = dx1 = 0.1
   nx = int(Lx/dx0) 
   ny = int(Ly/dx1)
   mesh = RectangleMesh(Point(x1, y1), Point(x1+Lx, y1+Ly), nx, ny,"right/left")
   V = FunctionSpace(mesh, 'Lagrange', 1)

   # parameters for the problem
   c0 = Constant(1.0)
   D0 = Constant(0.1)
   x0 = Constant(8.0)
   y0 = Constant(8.0)
   #check 'kn' function - it works
   hhh = kn(0,2)
   print('hhh = %s' % hhh)

  def exsol():
        p = Expression("sqrt( log(x[0]/x0)*log(x[0]/x0) + log(x[1]/y0)*log(x[1]/y0) )/c0", x0=x0, y0=y0, c0=c0) 
        q = Expression("(p*sqrt(1.0+c0*c0*D0*D0))/(sqrt(2.0)*D0)",p=p,c0=c0,D0=D0)
        uext = Expression("kn(0,q)",q=q) #this line does not get compiler in JIT
        ue = interpolate(uext, V)
        return ue

   # Plot the exact solution at t=t_start
   ue = exsol()
   plot(ue, interactive=True,
     title="Exact_solution_t0")

I get an error for the line "Expression("kn(0,q)",q=q)" - RuntimeError: In instant.recompile: The module did not compile with command 'make VERBOSE=1'

Please help me with this.

Answer 1

Here is a code that will solve this problem.

  from dolfin import *
  import numpy

  # System size definition 
  # system = (x1, x1+Lx) x (y1, y1+Ly)
  Lx = 2.0; Ly = 2.0
  x1 = 10.0; y1 = 10.0

  c0 = 1.0
  D0 = 1.0
  x0 = 1.5
  y0 = 1.5

  # Create mesh and define function space
  dx0 = dx1 = 0.1
  nx = int(Lx/dx0) 
  ny = int(Ly/dx1)
  mesh = RectangleMesh(Point(x1, y1), Point(x1+Lx, y1+Ly), nx, ny,"right/left")
  V = FunctionSpace(mesh, 'Lagrange', 1)

  code = '''
  #include <math.h>
  #include <boost/math/special_functions/bessel.hpp>
  using boost::math::cyl_bessel_i;
  using boost::math::cyl_bessel_j;
  using boost::math::cyl_bessel_k;
  using boost::math::cyl_neumann;

  namespace dolfin {
       class MyFun : public Expression
      {
              double c0,D0,x0,y0;
              public:
                        MyFun(): Expression() {};
              void eval(Array<double>& values, const Array<double>& x) const {
                      double f = D0*sqrt( logf(x[0]/x0)*logf(x[0]/x0) + logf(x[1]/y0)*logf(x[1]/y0) )/c0 ;
                      values[0] = cyl_bessel_k(0,f);
              }

              void update(double _c0, double _D0, double _x0, double _y0) {
                    c0 = _c0;
                    D0 = _D0;
                    x0 = _x0;
                    y0 = _y0;
            }
      };
}'''

  f=Expression(code)
  f.update(c0,D0,x0,y0)

  ss = interpolate(f,V)
  plot(ss, interactive=True)

Written with help from this forum post.