gen_f90.py

"""Generate F90 pre-computed SDC integration matrices."""

# Copyright (c) 2012.  Matthew Emmett.  All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions
# are met:
#
#   1. Redistributions of source code must retain the above copyright
#      notice, this list of conditions and the following disclaimer.
#
#   2. Redistributions in binary form must reproduce the above
#      copyright notice, this list of conditions and the following
#      disclaimer in the documentation and/or other materials provided
#      with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
# FOR A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE
# COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.

import collections
import sdcquad as sdc


# compute "refined" S matrices too?
compute_refinements = True

# set number of digits
digits = 40


##
## set global mpmath precision
##

from sympy.mpmath import mp
from sympy.mpmath import nstr
mp.dps    = digits
mp.pretty = True


##
## compute sdc nodes and smats
##

# build a dictionary (quads) for the various quadrature types:
#
# * G:  Gauss-Legendre,
# * U:  Uniform,
# * GL: Gauss-Lobatto,
# * GR: Gauss-Radau, and
# * CC: Clenshaw-Curtis

quads = {}
for q in [ 'G', 'GL', 'GR', 'CC' ]:
  quads[q] = {}

  # start with n0 nodes, and double the nodes maxsteps times
  for n0, maxsteps in [ (2, 5), (4, 2) ]:
    n = n0
    refine = [ 1 ]

    for i in range(maxsteps+1):
      quads[q][n] = {}

      for r in refine:
        print q, n, r
        try:
          (nodes, mask) = sdc.nodes(q, n)
          smat = sdc.smat(nodes[::r], mask[::r])
          qmat = sdc.qmat(nodes[::r], mask[::r])
          quads[q][n][r] = (nodes[::r], mask[::r], smat, qmat)
        except:
          print 'SKIPPED'
          pass

      if compute_refinements:
        refine.append(2**(i+1))

      n = 2*n-1


quads['U'] = {}
for n in range(2, 34):
  print 'U', n
  (nodes, mask) = sdc.nodes('U', n)
  smat = sdc.smat(nodes, mask)
  qmat = sdc.qmat(nodes, mask)
  quads['U'][n] = { 1: (nodes, mask, smat, qmat) }


##
## write fortran module 'quadrature.f90'
##

def fstr(x):
  s = nstr(x, n=digits)
  return s + '_pfdp'


with open('quadrature.f90', 'w') as f:
  f.write('! this file was generated by "mksmat.py"\n')
  f.write('module quadrature\ncontains\n')

  for q in quads:
    subroutine = 'sdcquad%s' % q
    f.writelines([ 'subroutine %s(nnodes, refine, smat, qmat, nodes, mask)\n' % subroutine,
                   'implicit none\n',
                   'integer, intent(in)  :: nnodes, refine\n',
                   'real(pfdp), intent(out) :: smat((nnodes-1)/refine,(nnodes-1)/refine+1)\n',
                   'real(pfdp), intent(out) :: qmat((nnodes-1)/refine,(nnodes-1)/refine+1)\n',
                   'real(pfdp), intent(out) :: nodes((nnodes-1)/refine+1)\n',
                   'logical, intent(out) :: mask((nnodes-1)/refine+1)\n',
                   '\n',
                   'select case(nnodes)\n' ])

    for n in sorted(quads[q]):
      f.writelines(['case(%d)\n' % n,
                    '  select case(refine)\n'])

      for r in sorted(quads[q][n]):
        f.writelines(['  case(%d)\n' % r])

        (nodes, mask, smat, qmat) = quads[q][n][r]

        for k, node in enumerate(nodes):
          f.write('    nodes(%d) = %s\n' % (k+1, fstr(node)))

        for k, node in enumerate(mask):
          if node:
            s = '.true.'
          else:
            s = '.false.'
          f.write('    mask(%d) = %s\n' % (k+1, s))

        for i, row in enumerate(smat):
          for j, entry in enumerate(row):
            f.write('    smat(%d,%d) = %s\n'
                    % (i+1, j+1, fstr(entry)))

        for i, row in enumerate(qmat):
          for j, entry in enumerate(row):
            f.write('    qmat(%d,%d) = %s\n'
                    % (i+1, j+1, fstr(entry)))

      f.writelines(['  end select\n'])

    f.writelines([ 'end select\n',
                   'end subroutine %s\n' % subroutine ])

  f.write('end module quadrature\n\n')