diff --git a/.flake8 b/.flake8 new file mode 100644 index 0000000..d9ad0b4 --- /dev/null +++ b/.flake8 @@ -0,0 +1,5 @@ +[flake8] +ignore = E203, E266, E501, W503, F403, F401 +max-line-length = 79 +max-complexity = 18 +select = B,C,E,F,W,T4,B9 diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 08902c8..b4159cb 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -1,24 +1,27 @@ name: CI on: - push: - branches: - - master - pull_request: + push: + branches: + - master + pull_request: jobs: - test: - name: QBMMTest - runs-on: ubuntu-latest - - steps: - - uses: actions/checkout@v2 - - name: Set up Python 3.8 - uses: actions/setup-python@v2 - with: - python-version: 3.8 - - name: Install dependencis - run: | - python -m pip install --upgrade pip - pip install -r requirements.txt - - name: "Test QBMM" - run: pytest test/test_qbmm.py + build: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v2 + - uses: actions/setup-python@v2 + with: + python-version: '3.8' + - name: Install pip + run: | + python -m pip install --upgrade pip + python -m pip install pipreqs + pipreqs --force + - name: Install requirements + run: python -m pip install -r requirements.txt + - name: Test QBMM + run: cd test && pytest test_qbmm.py + + + diff --git a/.gitignore b/.gitignore index 4f53ab9..0b3bbdd 100644 --- a/.gitignore +++ b/.gitignore @@ -4,6 +4,8 @@ tags *.ipynb_checkpoints D +src/output _build __pycache__ -mmatica_qbmmlib + +profile* diff --git a/README.md b/README.md index f0fef56..f287898 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,6 @@ # PyQBMMlib +![CI](https://github.com/sbryngelson/PyQBMMlib/workflows/CI/badge.svg) [![License: GPL v3](https://img.shields.io/badge/License-GPLv3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0) [![Documentation Status](https://readthedocs.org/projects/pyqbmmlib/badge/?version=latest)](https://pyqbmmlib.readthedocs.io/en/latest/?badge=latest) [![DOI](https://zenodo.org/badge/doi/10.1016/j.softx.2020.100615.svg)](http://dx.doi.org/10.1016/j.softx.2020.100615) diff --git a/docs/conf.py b/docs/conf.py index 139254e..6535e59 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -14,21 +14,22 @@ # import os import sys -sys.path.insert(0, os.path.abspath('../src/')) -sys.path.insert(0, os.path.abspath('../utils/')) -sys.path.insert(0, os.path.abspath('../test/')) + +sys.path.insert(0, os.path.abspath("../src/")) +sys.path.insert(0, os.path.abspath("../utils/")) +sys.path.insert(0, os.path.abspath("../test/")) # -- Project information ----------------------------------------------------- -project = u'PyQBMMlib' -copyright = u'2020, Spencer Bryngelson' -author = u'Spencer H. Bryngelson and Esteban Cisneros' +project = u"PyQBMMlib" +copyright = u"2020, Spencer Bryngelson" +author = u"Spencer H. Bryngelson and Esteban Cisneros" # The short X.Y version -version = u'' +version = u"" # The full version, including alpha/beta/rc tags -release = u'' +release = u"" # -- General configuration --------------------------------------------------- @@ -41,25 +42,25 @@ # extensions coming with Sphinx (named 'sphinx.ext.*') or your custom # ones. extensions = [ - 'sphinx.ext.autodoc', - 'sphinx.ext.doctest', - 'sphinx.ext.coverage', - 'sphinx.ext.mathjax', - 'sphinx.ext.inheritance_diagram', - 'sphinx.ext.napoleon' + "sphinx.ext.autodoc", + "sphinx.ext.doctest", + "sphinx.ext.coverage", + "sphinx.ext.mathjax", + "sphinx.ext.inheritance_diagram", + "sphinx.ext.napoleon", ] # Add any paths that contain templates here, relative to this directory. -templates_path = ['_templates'] +templates_path = ["_templates"] # The suffix(es) of source filenames. # You can specify multiple suffix as a list of string: # # source_suffix = ['.rst', '.md'] -source_suffix = '.rst' +source_suffix = ".rst" # The master toctree document. -master_doc = 'index' +master_doc = "index" # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. @@ -71,7 +72,7 @@ # List of patterns, relative to source directory, that match files and # directories to ignore when looking for source files. # This pattern also affects html_static_path and html_extra_path. -exclude_patterns = [u'_build', 'Thumbs.db', '.DS_Store'] +exclude_patterns = [u"_build", "Thumbs.db", ".DS_Store"] # The name of the Pygments (syntax highlighting) style to use. pygments_style = None @@ -82,7 +83,7 @@ # The theme to use for HTML and HTML Help pages. See the documentation for # a list of builtin themes. # -html_theme = 'classic' +html_theme = "classic" # html_theme = 'alabaster' # Theme options are theme-specific and customize the look and feel of a theme @@ -94,7 +95,7 @@ # Add any paths that contain custom static files (such as style sheets) here, # relative to this directory. They are copied after the builtin static files, # so a file named "default.css" will overwrite the builtin "default.css". -html_static_path = ['_static'] +html_static_path = ["_static"] # Custom sidebar templates, must be a dictionary that maps document names # to template names. @@ -110,7 +111,7 @@ # -- Options for HTMLHelp output --------------------------------------------- # Output file base name for HTML help builder. -htmlhelp_basename = 'PyQBMMlibdoc' +htmlhelp_basename = "PyQBMMlibdoc" # -- Options for LaTeX output ------------------------------------------------ @@ -119,15 +120,12 @@ # The paper size ('letterpaper' or 'a4paper'). # # 'papersize': 'letterpaper', - # The font size ('10pt', '11pt' or '12pt'). # # 'pointsize': '10pt', - # Additional stuff for the LaTeX preamble. # # 'preamble': '', - # Latex figure (float) alignment # # 'figure_align': 'htbp', @@ -137,8 +135,7 @@ # (source start file, target name, title, # author, documentclass [howto, manual, or own class]). latex_documents = [ - (master_doc, 'PyQBMMlib.tex', u'PyQBMMlib Documentation', - u'SHB', 'manual'), + (master_doc, "PyQBMMlib.tex", u"PyQBMMlib Documentation", u"SHB", "manual"), ] @@ -146,10 +143,7 @@ # One entry per manual page. List of tuples # (source start file, name, description, authors, manual section). -man_pages = [ - (master_doc, 'pyqbmmlib', u'PyQBMMlib Documentation', - [author], 1) -] +man_pages = [(master_doc, "pyqbmmlib", u"PyQBMMlib Documentation", [author], 1)] # -- Options for Texinfo output ---------------------------------------------- @@ -158,9 +152,15 @@ # (source start file, target name, title, author, # dir menu entry, description, category) texinfo_documents = [ - (master_doc, 'PyQBMMlib', u'PyQBMMlib Documentation', - author, 'PyQBMMlib', 'One line description of project.', - 'Miscellaneous'), + ( + master_doc, + "PyQBMMlib", + u"PyQBMMlib Documentation", + author, + "PyQBMMlib", + "One line description of project.", + "Miscellaneous", + ), ] @@ -179,7 +179,7 @@ # epub_uid = '' # A list of files that should not be packed into the epub file. -epub_exclude_files = ['search.html'] +epub_exclude_files = ["search.html"] # -- Extension configuration ------------------------------------------------- @@ -187,7 +187,7 @@ # -- Options for intersphinx extension --------------------------------------- # Example configuration for intersphinx: refer to the Python standard library. -intersphinx_mapping = {'https://docs.python.org/': None} +intersphinx_mapping = {"https://docs.python.org/": None} # -- Options for todo extension ---------------------------------------------- diff --git a/inputs/example_2d.yaml b/inputs/example_2d.yaml index 13c9788..294a8a7 100644 --- a/inputs/example_2d.yaml +++ b/inputs/example_2d.yaml @@ -6,17 +6,17 @@ qbmm: flow: no governing_dynamics: " - x - xdot" num_internal_coords: 2 - num_quadrature_nodes: 9 + num_quadrature_nodes: 4 method: "chyqmom" adaptive: no max_skewness: 30 checks: no advancer: - method: "RK23" + method: "RK3" time_step: 1.0e-5 final_time: 15.0 - error_tol: 1.0e-7 + error_tol: 1.0e-6 num_steps: 10000 num_steps_print: 100 num_steps_write: 100 diff --git a/notebooks/Untitled.ipynb b/notebooks/Untitled.ipynb new file mode 100644 index 0000000..1e3f0e2 --- /dev/null +++ b/notebooks/Untitled.ipynb @@ -0,0 +1,92 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "some_string\n", + "y is empty\n" + ] + } + ], + "source": [ + "x = \"some_string\"\n", + "y = \"\"\n", + "\n", + "if x:\n", + " print(x)\n", + "else:\n", + " print(\"x is empty\")\n", + " \n", + "if y:\n", + " print(y)\n", + "else:\n", + " print(\"y is empty\")\n", + " " + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "['CavityFlame.aux', 'CavityFlame.blg']" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "project_name = \"CavityFlame\"\n", + "[project_name+ext for ext in [\".aux\",\".blg\"]]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 2", + "language": "python", + "name": "python2" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.3" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/notebooks/Untitled1.ipynb b/notebooks/Untitled1.ipynb new file mode 100644 index 0000000..c1ca2cf --- /dev/null +++ b/notebooks/Untitled1.ipynb @@ -0,0 +1,200 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import itertools" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "X = np.zeros(402)\n", + "X[-1] = 1" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(48,)" + ] + }, + "execution_count": 15, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "X[:48].shape\n", + "X[-48:].shape" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", + " 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", + " 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.])" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "X[-48:]" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0 0\n", + "1 1\n", + "2 2\n", + "3 3\n", + "4 4\n", + "5 5\n", + "6 6\n", + "7 7\n", + "8 8\n", + "9 9\n", + "10 10\n", + "11 11\n", + "12 12\n", + "13 13\n", + "14 14\n", + "15 15\n" + ] + } + ], + "source": [ + "num_moments = 16\n", + "num_nodes = 27\n", + "\n", + "for i_moment, i_node in enumerate(range(num_moments), range(num_nodes)):\n", + " print(i_moment, i_node)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[ 1. 2. 3. 4. 5. 6. 7. 8. 9. 10.]\n" + ] + } + ], + "source": [ + "X = np.linspace(1, 10, 10)\n", + "print(X)" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0 , 0 , m+n+1 1\n", + "0 , 1 , m+n+1 2\n", + "0 , 2 , m+n+1 3\n", + "0 , 3 , m+n+1 4\n", + "1 , 0 , m+n+1 2\n", + "1 , 1 , m+n+1 3\n", + "1 , 2 , m+n+1 4\n", + "1 , 3 , m+n+1 5\n" + ] + } + ], + "source": [ + "num_moments = 2\n", + "num_nodes = 4\n", + "\n", + "flux = np.zeros(num_moments)\n", + "\n", + "for m, n in itertools.product(range(num_moments), range(num_nodes)):\n", + " print(m,\", \",n,\", \",\"m+n+1\",m+n+1)\n", + " flux[m] += m+n+1" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([10., 14.])" + ] + }, + "execution_count": 16, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "flux" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 2", + "language": "python", + "name": "python2" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.3" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..3f64c0d --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,10 @@ +[tool.black] +exclude = ''' +/( + \.git + | \.github + | docs + | inputs + | notebooks +)/ +''' diff --git a/requirements.txt b/requirements.txt index 19434e2..1099ca7 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,92 +1,6 @@ -appdirs==1.4.4 -appnope @ file:///opt/concourse/worker/volumes/live/0291c9e1-4b15-459f-623e-2770f55be269/volume/appnope_1594338395037/work -argon2-cffi @ file:///opt/concourse/worker/volumes/live/59af29ac-4890-416e-7ab7-794f8d6f7ecd/volume/argon2-cffi_1596828548321/work -async-generator==1.10 -attrs @ file:///tmp/build/80754af9/attrs_1604765588209/work -backcall==0.2.0 -bleach @ file:///tmp/build/80754af9/bleach_1600439572647/work -Cantera==2.5.0b1 -certifi==2020.6.20 -cffi==1.14.0 -chardet==3.0.4 -conda==4.9.2 -conda-package-handling==1.7.0+0.g7c4a471.dirty -cryptography==2.9.2 -cycler==0.10.0 -decorator==4.4.2 -defusedxml==0.6.0 -entrypoints==0.3 -h5py @ file:///opt/concourse/worker/volumes/live/ca975206-a8f5-4947-64c0-5af2e33e601d/volume/h5py_1593454180658/work -idna==2.9 -importlib-metadata @ file:///tmp/build/80754af9/importlib-metadata_1602276842396/work -iniconfig==1.1.1 -ipykernel @ file:///opt/concourse/worker/volumes/live/88f541d3-5a27-498f-7391-f2e50ca36560/volume/ipykernel_1596206680118/work/dist/ipykernel-5.3.4-py3-none-any.whl -ipython @ file:///opt/concourse/worker/volumes/live/26969e8f-c9f7-42dc-6ffb-b3effd424c49/volume/ipython_1604101242376/work -ipython-genutils==0.2.0 -ipywidgets @ file:///tmp/build/80754af9/ipywidgets_1601490159889/work -jedi @ file:///opt/concourse/worker/volumes/live/c63c70ef-654a-4d26-6b78-37aef096f225/volume/jedi_1596490811751/work -Jinja2==2.11.2 -jsonschema @ file:///tmp/build/80754af9/jsonschema_1602607155483/work -jupyter==1.0.0 -jupyter-client @ file:///tmp/build/80754af9/jupyter_client_1601311786391/work -jupyter-console @ file:///tmp/build/80754af9/jupyter_console_1598884538475/work -jupyter-core==4.6.3 -jupyterlab-pygments @ file:///tmp/build/80754af9/jupyterlab_pygments_1601490720602/work -kiwisolver==1.3.1 -Mako==1.1.3 -MarkupSafe @ file:///opt/concourse/worker/volumes/live/cb778296-98db-45ad-411e-6f726e102dc3/volume/markupsafe_1594371638608/work -matplotlib==3.3.3 -mistune @ file:///opt/concourse/worker/volumes/live/95802d64-d39c-491b-74ce-b9326880ca54/volume/mistune_1594373201816/work -mkl-fft==1.2.0 -mkl-random==1.1.1 -mkl-service==2.3.0 -mpmath==1.1.0 -nbclient @ file:///tmp/build/80754af9/nbclient_1602783176460/work -nbconvert @ file:///opt/concourse/worker/volumes/live/2b9c1d93-d0fd-432f-7d93-66c93d81b614/volume/nbconvert_1601914875037/work -nbformat @ file:///tmp/build/80754af9/nbformat_1602783287752/work -nest-asyncio @ file:///tmp/build/80754af9/nest-asyncio_1605115881283/work -notebook @ file:///opt/concourse/worker/volumes/live/be0f3504-189d-4bae-4e57-c5d6da73ffcd/volume/notebook_1601501605350/work -numpy @ file:///opt/concourse/worker/volumes/live/5572694e-967a-4c0c-52cf-b53d43e72de9/volume/numpy_and_numpy_base_1603491881791/work -packaging==20.4 -pandocfilters @ file:///opt/concourse/worker/volumes/live/c330e404-216d-466b-5327-8ce8fe854d3a/volume/pandocfilters_1605120442288/work -parso==0.7.0 -pexpect @ file:///opt/concourse/worker/volumes/live/8701bb20-ad87-46c7-5108-30c178cf97e5/volume/pexpect_1594383388344/work -pickleshare @ file:///opt/concourse/worker/volumes/live/93ec39d8-05bb-4f84-7efc-98735bc39b70/volume/pickleshare_1594384101884/work -Pillow==8.0.1 -pluggy==0.13.1 -prometheus-client==0.8.0 -prompt-toolkit @ file:///tmp/build/80754af9/prompt-toolkit_1602688806899/work -ptyprocess==0.6.0 -py==1.9.0 -pycosat==0.6.3 -pycparser==2.20 -Pygments @ file:///tmp/build/80754af9/pygments_1604103097372/work -pymbolic==2020.1 -pyOpenSSL==19.1.0 -pyparsing==2.4.7 -pyrometheus==2020.1 -pyrsistent @ file:///opt/concourse/worker/volumes/live/ff11f3f0-615b-4508-471d-4d9f19fa6657/volume/pyrsistent_1600141727281/work -PySocks==1.7.1 -pytest==6.1.2 -python-dateutil==2.8.1 -pytools==2020.4.3 -pyzmq==19.0.2 -qtconsole @ file:///tmp/build/80754af9/qtconsole_1600870028330/work -QtPy==1.9.0 -requests==2.23.0 -ruamel-yaml==0.15.87 -scipy==1.5.4 -Send2Trash==1.5.0 -six==1.15.0 +numpy==1.19.4 sympy==1.7.1 -terminado==0.9.1 -testpath==0.4.4 -toml==0.10.2 -tornado==6.0.4 -tqdm==4.46.0 -traitlets @ file:///tmp/build/80754af9/traitlets_1602787416690/work -urllib3==1.25.8 -wcwidth @ file:///tmp/build/80754af9/wcwidth_1593447189090/work -webencodings==0.5.1 -widgetsnbextension==3.5.1 -zipp @ file:///tmp/build/80754af9/zipp_1604001098328/work +scipy==1.5.4 +numba==0.52.0 +pytest==6.2.0 +PyYAML==5.3.1 diff --git a/src/advancer.py b/src/advancer.py index f142931..57bf4f6 100644 --- a/src/advancer.py +++ b/src/advancer.py @@ -10,23 +10,26 @@ """ import sys -sys.path.append('../utils/') + +sys.path.append("../utils/") from stats_util import * from pretty_print_util import * -from qbmm_manager import * +from simulation_domain import * import numpy as np import csv +import h5py + class time_advancer: """This class advances the moment transport equations in time. - - A ``config`` dictionary is required to usethis class. - Example files are provided in ``inputs/``. + + A ``config`` dictionary is required to usethis class. + Example files are provided in ``inputs/``. The dictionary is used to set the following variables: - :ivar method: Time integration scheme (``euler`` or ``RK23``) - :ivar time_step: Time integration time step + :ivar method: Time integration scheme (``euler``, ``RK2, or ``RK3``) + :ivar time_step: Time integration time step :ivar final_time: Final integration time :ivar num_steps: Number of integration steps :ivar num_steps_print: Report-status frequency @@ -34,8 +37,8 @@ class time_advancer: :ivar output_dir: Directory to which output is written :ivar output_id: An ID for each run's output file - Create an advancer object simply by typing - + Create an advancer object simply by typing + >>> advancer = time_advancer( config ) Before running, you'll need to initialize the state. One way you may do this is through: @@ -44,69 +47,101 @@ class time_advancer: for a 1D problem, where ``mu`` and ``sigma`` are user-specified. Now all you have to do is: - >>> advancer.run() + >>> advancer.run() """ def __init__(self, config): """Constructor - + :param config: Configuration :type config: dict """ - - + self.method = config['advancer']['method'] self.time_step = config['advancer']['time_step'] self.final_time = config['advancer']['final_time'] - self.error_tol = config['advancer']['error_tol'] self.num_steps = config['advancer']['num_steps'] self.num_steps_print = config['advancer']['num_steps_print'] self.num_steps_write = config['advancer']['num_steps_write'] self.output_dir = config['advancer']['output_dir'] self.output_id = config['advancer']['output_id'] self.write_to = config['advancer']['write_to'] - - self.qbmm_mgr = qbmm_manager( config ) - - self.num_dim = self.qbmm_mgr.num_moments - self.indices = self.qbmm_mgr.indices - - self.state = np.zeros( self.num_dim ) - self.rhs = np.zeros( self.num_dim ) - self.stage_state = np.zeros( [3, self.num_dim] ) - self.stage_k = np.zeros( [3, self.num_dim] ) + if 'project' in config['advancer']: + self.project = config['advancer']['project'] + else: + self.project = True - if self.method == 'Euler': + self.domain = simulation_domain(config) + + if self.method == "Euler": self.advance = self.advance_euler - elif self.method == 'RK23': + self.num_stages = 1 + elif self.method == "RK2": + self.advance = self.advance_RK2 + self.num_stages = 2 + elif self.method == "RK23": self.advance = self.advance_RK23 + self.num_stages = 3 + + num_points = self.domain.num_points + num_moments = self.domain.qbmm_mgr.num_moments + num_coords = self.domain.qbmm_mgr.num_coords + self.state = np.zeros([num_points, num_moments]) + self.rhs = np.zeros([num_points, num_moments]) + self.stage_state = np.zeros([self.num_stages, num_points, num_moments]) + self.stage_k = np.zeros([self.num_stages, num_points, num_moments]) + + self.max_time_step = 1.e5 + self.min_time_step = self.time_step + if "error_tol" in config["advancer"]: + if self.method == "Euler": + print("Euler method is non-adaptive. Ignoring error_tol.") + self.adaptive = False + elif self.method == "RK2": + print("To run RK2 adaptively, specify cfl, ignoring error_tol.") + self.adaptive = False + else: + self.adaptive = True + self.error_tol = config["advancer"]["error_tol"] + self.adapt_time_step = self.adapt_time_step_RK23 + else: + self.adaptive = False + + if "cfl" in config["advancer"]: + if self.method == "RK23": + print("RK23 adapts time step per error_tol, ignoring cfl") + else: + self.cfl = config["advancer"]["cfl"] + self.adaptive = True + self.adapt_time_step = self.adapt_time_step_cfl + else: + self.cfl = 1 + + # Report print('advancer: init: Configuration options ready') print('\t method = %s' % self.method) print('\t time_step = %.4E' % self.time_step) print('\t final_time = %.4E' % self.final_time) - print('\t error_tol = %.4E' % self.error_tol) + #print('\t error_tol = %.4E' % self.error_tol) print('\t num_steps_print = %i' % self.num_steps_print) print('\t num_steps_write = %i' % self.num_steps_write) print('\t output_dir = %s' % self.output_dir) print('\t output_id = %s' % self.output_id) print('\t write_to = %s' % self.write_to) - self.max_time_step = 1.e5 - self.min_time_step = self.time_step - + # Output options self.file_name = self.output_dir + 'qbmm_state_' + self.output_id if self.write_to == 'txt': self.file_name += '.dat' self.write_to_file = self.write_to_txt - elif self.write_to == 'h5': - self.file_name += '.h5' + elif self.write_to == "h5": + self.file_name += ".h5" self.write_to_file = self.write_to_h5 - return - + def initialize_state(self, init_state): """ This function initializes the advancer state @@ -114,13 +149,20 @@ def initialize_state(self, init_state): :param init_state: Initial condition :type init_state: array like """ - + self.state = init_state - print('state', self.state) - + print("state", self.state) + return + def initialize_state_jets(self): + """ + This function intializes the advancer state to R.O. Fox's jet conditions + """ + self.domain.initialize_state_jets(self.state) + return + def initialize_state_gaussian_trivar(self, mu1, mu2, mu3, sig1, sig2, sig3): """ This function initializes the state to the raw moments of a trivariate Gaussian distribution. @@ -133,19 +175,17 @@ def initialize_state_gaussian_trivar(self, mu1, mu2, mu3, sig1, sig2, sig3): :param sig3: Standard deviation of coordinate 3 :type mu1: float :type mu2: float - :type mu3: floag + :type mu3: float :type sig1: float :type sig2: float :type sig3: float """ - - self.state = raw_gaussian_moments_trivar( self.indices, mu1, mu2, mu3, sig1, sig2, sig3 ) + self.state[:] = raw_gaussian_moments_trivar(self.qbmm_mgr.indices, mu1, mu2, mu3, sig1, sig2, sig3) message = 'advancer: initialize_trigaussian: ' - f_array_pretty_print( message, 'state', self.state ) + f_array_pretty_print(message, 'state', self.state) return - def initialize_state_gaussian_bivar(self, mu1, mu2, sigma1, sigma2): """ This function initializes the state to the raw moments of a bivariate Gaussian distribution. @@ -159,11 +199,9 @@ def initialize_state_gaussian_bivar(self, mu1, mu2, sigma1, sigma2): :type sig1: float :type sig2: float """ - - - self.state = raw_gaussian_moments_bivar( self.indices, mu1, mu2, sigma1, sigma2 ) + self.state[:] = raw_gaussian_moments_bivar(self.qbmm_mgr.indices, mu1, mu2, sigma1, sigma2) message = 'advancer: initialize_bigaussian: ' - f_array_pretty_print( message, 'state', self.state ) + f_array_pretty_print(message, 'state', self.state) return def initialize_state_gaussian_univar(self, mu, sigma): @@ -174,81 +212,145 @@ def initialize_state_gaussian_univar(self, mu, sigma): :param sigma: Standard deviation """ - self.state = raw_gaussian_moments_univar( self.num_dim, mu, sigma) - message = 'advancer: initialize_gaussian: ' - f_array_pretty_print( message, 'state', self.state ) + self.state[:] = raw_gaussian_moments_univar(self.num_dim, mu, sigma) + message = "advancer: initialize_gaussian: " + f_array_pretty_print(message, "state", self.state) + return + + def initialize_flow(self): + """ + This function initializes flow + """ + self.simulation_domain return def advance_euler(self): """ This function advances the state with an explicit Euler scheme """ - # self.rhs += self.qbmm_mgr.evaluate_rhs( self.state ) + # Stage 1: { y_1, k_1 } = f( t_n, y_0 ) + self.stage_state[0] = self.state.copy() + self.stage_k[0] = self.domain.compute_rhs(self.stage_state[0]) + + # Project back to moment set using computed weights/abscissas + if self.project: + self.domain.project(self.stage_state[0]) + + # Updates + self.state = self.stage_state[0] + self.time_step * self.stage_k[0] + + self.time_step = self.cfl * self.domain.grid_spacing / self.domain.max_abscissa() + + return + + + def advance_RK2(self): + """ + This function advances the state with a Runge--Kutta 2 scheme + """ + + # Stage 1: + self.stage_state[0] = self.state.copy() + self.stage_k[0] = self.domain.compute_rhs(self.stage_state[0]) + + if self.project: + self.domain.project(self.stage_state[0]) - print('advancer: advance: Euler time advancer') + self.stage_state[1] = self.stage_state[0] + self.time_step * self.stage_k[0] - self.rhs = self.qbmm_mgr.compute_rhs( self.state ) + # Stage 2: + self.stage_k[1] = self.domain.compute_rhs(self.stage_state[1]) - ### state_{n+1} = proj_{n} + rhs_{n} - self.state += self.time_step * self.rhs + if self.project: + self.domain.project(self.stage_state[1]) - return + # Update + self.state = 0.5*(self.stage_state[0] + self.stage_state[1] + + self.time_step*self.stage_k[1]) - def advance_RK23(self): + return + + def advance_RK3(self): """ This function advances the state with a Runge--Kutta 2/3 scheme """ # Stage 1: { y_1, k_1 } = f( t_n, y_0 ) self.stage_state[0] = self.state.copy() - time = self.time - self.qbmm_mgr.compute_rhs( self.stage_state[0], self.stage_k[0] ) + self.qbmm_mgr.compute_rhs(self.stage_state[0], self.stage_k[0]) self.stage_state[1] = self.stage_state[0] + self.time_step * self.stage_k[0] # Stage 2: { y_2, k_2 } = f( t_n, y_1 + dt * k_1 ) - time = self.time + self.time_step - self.qbmm_mgr.compute_rhs( self.stage_state[1], self.stage_k[1] ) - test_state = 0.5 * ( self.stage_state[0] + ( self.stage_state[1] + self.time_step * self.stage_k[1] ) ) + self.qbmm_mgr.compute_rhs(self.stage_state[1], self.stage_k[1]) + test_state = 0.5 * ( + self.stage_state[0] + + (self.stage_state[1] + self.time_step * self.stage_k[1]) + ) # Stage 3: { y_3, k_3 } = f( t_n + 0.5 * dt, ... ) - self.stage_state[2] = 0.75 * self.stage_state[0] + 0.25 * ( self.stage_state[1] + self.time_step * self.stage_k[1] ) - self.qbmm_mgr.compute_rhs( self.stage_state[2], self.stage_k[2] ) + self.stage_state[2] = 0.75 * self.stage_state[0] + 0.25 * ( + self.stage_state[1] + self.time_step * self.stage_k[1] + ) + self.qbmm_mgr.compute_rhs(self.stage_state[2], self.stage_k[2]) # Updates - self.state = ( self.stage_state[0] + 2.0 * ( self.stage_state[2] + self.time_step * self.stage_k[2] ) ) / 3.0 - self.rk_error = np.linalg.norm( self.state - test_state ) / np.linalg.norm( self.state ) + self.state = ( + self.stage_state[0] + + 2.0 * (self.stage_state[2] + self.time_step * self.stage_k[2]) + ) / 3.0 + + if self.adaptive: + self.ts_error = np.linalg.norm(self.state - test_state) / np.linalg.norm( + self.state + ) + + return + - def adapt_time_step(self): + def adapt_time_step_cfl(self): + """ + This function adapts the time step for flow problems based on the CFL condition + """ + # print('before t_step: ', self.time_step) + self.time_step = self.cfl*self.domain.grid_spacing/self.domain.max_abscissa() + # print(' biggest abscissa ', self.domain.max_abscissa()) + # print('after t_step: ', self.time_step) + return + + + def adapt_time_step_RK23(self): """ This function adapts the time step according to user-specified tolerance """ - error_fraction = np.sqrt( 0.5 * self.error_tol / self.rk_error ) - time_step_factor = min( max( error_fraction, 0.3 ), 2.0 ) - new_time_step = time_step_factor * self.time_step - new_time_step = min( max( 0.9 * new_time_step, self.min_time_step ), self.max_time_step ) - self.time_step = new_time_step + error_fraction = np.sqrt(0.5 * self.error_tol / self.ts_error) + time_step_factor = min(max(error_fraction, 0.3), 2.0) + new_time_step = time_step_factor * self.time_step + new_time_step = min( + max(0.9 * new_time_step, self.min_time_step), self.max_time_step + ) + self.time_step = new_time_step return - + def run(self): """ Advancer driver """ - print('advancer: run: Preparing to step') - + print("advancer: run: Preparing to step") + self.time = 0.0 self.report_step(0) self.write_step(0) i_step = 0 - step = True + step = True while step == True: self.advance() - i_step += 1 + i_step += 1 self.time += self.time_step if i_step % self.num_steps_print == 0: @@ -256,15 +358,16 @@ def run(self): if i_step % self.num_steps_write == 0: self.write_step(i_step) - - self.adapt_time_step() + + if self.adaptive: + self.adapt_time_step() if self.time > self.final_time or i_step >= self.num_steps: step = False - - print('advancer: run: stepping ends') - print('advancer: Number of steps:',i_step) - + + print("advancer: run: stepping ends") + print("advancer: Number of steps:", i_step) + return def report_step(self, i_step): @@ -275,8 +378,17 @@ def report_step(self, i_step): :type i_step: int """ - message = 'advancer: step = ' + str(i_step) + ' ... time = ' + '{:.16E}'.format(self.time) + ' ... time_step = ' + '{:.16E}'.format(self.time_step) + ' ... ' - f_array_pretty_print( message, 'state', self.state ) + message = ( + "advancer: step = " + + str(i_step) + + " ... time = " + + "{:.16E}".format(self.time) + + " ... time_step = " + + "{:.16E}".format(self.time_step) + + " ... " + ) + print(message) + #f_array_pretty_print(message, "state", self.state) def write_step(self, i_step): """ @@ -286,10 +398,10 @@ def write_step(self, i_step): :type i_step: int """ - message = 'advancer: step = ' + str(i_step) + ' ... Writing to file' + message = "advancer: step = " + str(i_step) + " ... Writing to file" # print(message) - self.write_to_file( i_step ) - + self.write_to_file(i_step) + def write_to_txt(self, i_step): """ This function writes the current state to a txt file @@ -297,19 +409,29 @@ def write_to_txt(self, i_step): :param i_step: Current step :type i_step: int """ - write_flag = 'a' + write_flag = "a" if i_step == 0: - write_flag = 'w' - - with open( self.file_name, write_flag ) as file_id: - # csv.writer( file_id, delimiter=' ' ).writerow( self.state ) - csv.writer( file_id, delimiter=' ' ).writerow( [self.time] + self.state.tolist() ) - + write_flag = "w" + + with open(self.file_name, write_flag) as file_id: + # csv.writer( file_id, delimiter=' ' ).writerow( self.state ) + csv.writer(file_id, delimiter=" ").writerow( + [self.time] + self.state.tolist() + ) + return - - def write_to_h5(self): + + def write_to_h5(self, i_step): """ - This function writes the current state to a h5 file, but is not implemented yet. + This function writes the current state to a h5 file. It is recommended for flow problems. """ - print('advancer: write_to_h5: not implemented yet') + if i_step == 0: + write_flag = "w" + else: + write_flag = "a" + with h5py.File(self.file_name, write_flag) as f: + g = f.create_group("step_"+str(i_step)) + g["time"] = self.time + g["time_step"] = self.time_step + g["state"] = self.state return diff --git a/src/config_manager.py b/src/config_manager.py index c56d725..5cca422 100644 --- a/src/config_manager.py +++ b/src/config_manager.py @@ -1,17 +1,17 @@ import yaml -class config_manager(): +class config_manager: def __init__(self, config_file): - - iret = self.load_config( config_file ) + + iret = self.load_config(config_file) return def load_config(self, config_file): - print('config_mgr: Loading config from %s' % config_file ) - with open( config_file ) as cf: - self.config = yaml.load( cf, Loader = yaml.Loader ) + print("config_mgr: Loading config from %s" % config_file) + with open(config_file) as cf: + self.config = yaml.load(cf, Loader=yaml.Loader) return diff --git a/src/devel_driver.py b/src/devel_driver.py index 8cf1ef3..4b67d47 100644 --- a/src/devel_driver.py +++ b/src/devel_driver.py @@ -12,147 +12,244 @@ from advancer import * from config_manager import * +from simulation_domain import * import sys -sys.path.append('../utils/') + +sys.path.append("../utils/") from stats_util import * -from euler_util import * from jets_util import * from pretty_print_util import * import cProfile -def flow_example(): +import warnings +warnings.filterwarnings('error') + +def flow_example_3d(): """ - This driver solves a flow-coupled problem. + This driver solves a flow-coupled problem. Currently, it only computes moment fluxes, but work is underway to solve the compressible flow equations. """ # In development + + cfl = 0.4 + N = 50000 + dx = 1./N + U_max = 1 + config = {} config['qbmm'] = {} - config['advancer'] = {} + config["advancer"] = {} + config["domain"] = {} - config['qbmm']['flow'] = True - config['qbmm']['governing_dynamics'] = '' - config['qbmm']['num_internal_coords'] = 3 - config['qbmm']['num_quadrature_nodes'] = 27 - config['qbmm']['method'] = 'chyqmom' - config['qbmm']['adaptive'] = False - config['qbmm']['max_skewness'] = 30 + config["qbmm"]["internal_dynamics"] = "" + config["qbmm"]["num_coords"] = 3 + config["qbmm"]["num_nodes"] = 27 + config["qbmm"]["method"] = "chyqmom" + config["qbmm"]["adaptive"] = False + config["qbmm"]["max_skewness"] = 30 - qbmm_mgr = qbmm_manager( config ) - indices = qbmm_mgr.indices + config["domain"]["flow"] = True + config["domain"]["num_dim"] = 1 + config["domain"]["num_points"] = N + 2 + config["domain"]["grid_extents"] = [0, 1] + + # config["advancer"]["method"] = "Euler" + config["advancer"]["method"] = "RK2" + config["advancer"]["time_step"] = cfl * dx / U_max + config["advancer"]["cfl"] = cfl + # config["advancer"]["project"] = False + config["advancer"]["final_time"] = 0.1 + config["advancer"]["num_steps"] = 10000 + config["advancer"]["num_steps_print"] = 1 #1000 + config["advancer"]["num_steps_write"] = 1000 #1000 + config["advancer"]["output_dir"] = "output/" + config["advancer"]["output_id"] = "example_flow_compiled" + config["advancer"]["write_to"] = "h5" - # Initial condition - # mu1 = 1.0 - # mu2 = 1.0 - # mu3 = 1.0 - # sigma1 = 0.1 - # sigma2 = 0.1 - # sigma3 = 0.1 - # moments = raw_gaussian_moments_trivar( indices, mu1, mu2, mu3, - # sigma1, sigma2, sigma3 ) - - moments_left, moments_right = jet_initialize_moments( qbmm_mgr ) - - xi_left, wts_left = qbmm_mgr.moment_invert( moments_left, indices ) - xi_right, wts_right = qbmm_mgr.moment_invert( moments_right, indices ) - - print(wts_left) - print(wts_right) + advancer = time_advancer(config) + advancer.initialize_state_jets() + advancer.run() + + +def flow_example_2d(): + """ + This driver solves a flow-coupled problem. + Currently, it only computes moment fluxes, but work is underway to solve the compressible flow equations. + """ + # In development + + cfl = 0.4 + N = 400 + dx = 1./N + U_max = 1 - print(xi_left) - print(xi_right) + config = {} + config['qbmm'] = {} + config["advancer"] = {} + config["domain"] = {} - flux = moment_fluxes( indices, wts_left, wts_right, xi_left, xi_right ) + config["qbmm"]["internal_dynamics"] = "" + config["qbmm"]["num_coords"] = 2 + config["qbmm"]["num_nodes"] = 9 + config["qbmm"]["method"] = "chyqmom" + config["qbmm"]["adaptive"] = False + config["qbmm"]["max_skewness"] = 30 - print(flux) + config["domain"]["flow"] = True + config["domain"]["num_dim"] = 1 + config["domain"]["num_points"] = N + 2 + config["domain"]["grid_extents"] = [0, 1] + + # config["advancer"]["method"] = "Euler" + config["advancer"]["method"] = "RK2" + config["advancer"]["time_step"] = cfl * dx / U_max + config["advancer"]["cfl"] = cfl + # config["advancer"]["project"] = False + config["advancer"]["final_time"] = 0.6 + config["advancer"]["num_steps"] = 10000 + config["advancer"]["num_steps_print"] = 1 #1000 + config["advancer"]["num_steps_write"] = 10 #1000 + config["advancer"]["output_dir"] = "output/" + config["advancer"]["output_id"] = "example_flow_compiled" + config["advancer"]["write_to"] = "h5" + + advancer = time_advancer(config) + advancer.initialize_state_jets() + advancer.run() + + + +def flow_example_1d(): + """ + This driver solves a flow-coupled problem. + Currently, it only computes moment fluxes, but work is underway to solve the compressible flow equations. + """ + # In development + + cfl = 0.4 + N = 40 + dx = 1./N + U_max = 1 + + config = {} + config['qbmm'] = {} + config["advancer"] = {} + config["domain"] = {} + + config["qbmm"]["internal_dynamics"] = "" + config["qbmm"]["num_coords"] = 1 + config["qbmm"]["num_nodes"] = 3 + # config["qbmm"]["method"] = "chyqmom" + config["qbmm"]["method"] = "hyqmom" + config["qbmm"]["adaptive"] = False + config["qbmm"]["max_skewness"] = 30 + + config["domain"]["flow"] = True + config["domain"]["num_dim"] = 1 + config["domain"]["num_points"] = N + 2 + config["domain"]["grid_extents"] = [0, 1] + + config["advancer"]["method"] = "Euler" + # config["advancer"]["method"] = "RK2" + config["advancer"]["time_step"] = cfl * dx / U_max + config["advancer"]["cfl"] = cfl + # config["advancer"]["project"] = False + config["advancer"]["final_time"] = 0.1 + config["advancer"]["num_steps"] = 10000 + config["advancer"]["num_steps_print"] = 1 #1000 + config["advancer"]["num_steps_write"] = 10 #1000 + config["advancer"]["output_dir"] = "output/" + config["advancer"]["output_id"] = "example_flow_compiled" + config["advancer"]["write_to"] = "h5" + + advancer = time_advancer(config) + advancer.initialize_state_jets() + advancer.run() - return def advance_example(config): """ - This driver solves the moments transport equations. - It is independent of governing dynamics, which are specified in the ``config`` dictionary. - It is constrained to problems with 1 and 2 internal coordinates by the initial condition. - The initial condition is Gaussian. + This driver solves the moments transport equations. + It is independent of governing dynamics, which are specified in the ``config`` dictionary. + It is constrained to problems with 1 and 2 internal coordinates by the initial condition. + The initial condition is Gaussian. :param config: Configuration dictionary - :type config: dict + :type config: dict """ - advancer = time_advancer( config ) + advancer = time_advancer(config) # Initialize condition - num_dim = config['qbmm']['num_internal_coords'] - mu = config['init_condition']['mu'] - sigma = config['init_condition']['sigma'] + num_dim = config["qbmm"]["num_coords"] + mu = config["init_condition"]["mu"] + sigma = config["init_condition"]["sigma"] if num_dim == 1: - advancer.initialize_state_gaussian_univar( mu, sigma ) + advancer.initialize_state_gaussian_univar(mu, sigma) elif num_dim == 2: mu1 = mu[0] mu2 = mu[1] sigma1 = sigma[0] sigma2 = sigma[1] - advancer.initialize_state_gaussian_bivar( mu1, mu2, sigma1, sigma2 ) + advancer.initialize_state_gaussian_bivar(mu1, mu2, sigma1, sigma2) # Run advancer.run() - - return + def advance_example2dp1(): # In development! config = {} - config['qbmm'] = {} - config['advancer'] = {} - - config['qbmm']['flow'] = False - config['qbmm']['governing_dynamics'] = ' - x - xdot - r0' - config['qbmm']['num_internal_coords'] = 2 - config['qbmm']['num_quadrature_nodes'] = 4 - config['qbmm']['method'] = 'chyqmom' - config['qbmm']['adaptive'] = False - config['qbmm']['max_skewness'] = 30 - config['qbmm']['polydisperse'] = True - config['qbmm']['num_poly_nodes'] = 5 - config['qbmm']['poly_symbol'] = 'r0' - - config['advancer']['method'] = 'RK23' - config['advancer']['time_step'] = 1.e-5 - config['advancer']['final_time'] = 30. - config['advancer']['error_tol'] = 1.0e-5 - config['advancer']['num_steps'] = 20000 - config['advancer']['num_steps_print'] = 1000 - config['advancer']['num_steps_write'] = 1000 - config['advancer']['output_dir'] = 'D/' - config['advancer']['output_id'] = 'example_2D' - config['advancer']['write_to'] = 'txt' - - advancer = time_advancer( config ) + config["qbmm"] = {} + config["advancer"] = {} + + config["qbmm"]["flow"] = False + config["qbmm"]["internal_dynamics"] = " - x - xdot - r0" + config["qbmm"]["num_coords"] = 2 + config["qbmm"]["num_nodes"] = 4 + config["qbmm"]["method"] = "chyqmom" + config["qbmm"]["adaptive"] = False + config["qbmm"]["max_skewness"] = 30 + config["qbmm"]["polydisperse"] = True + config["qbmm"]["num_poly_nodes"] = 5 + config["qbmm"]["poly_symbol"] = "r0" + + config["advancer"]["method"] = "RK23" + config["advancer"]["time_step"] = 1.0e-5 + config["advancer"]["final_time"] = 30.0 + config["advancer"]["error_tol"] = 1.0e-5 + config["advancer"]["num_steps"] = 20000 + config["advancer"]["num_steps_print"] = 1000 + config["advancer"]["num_steps_write"] = 1000 + config["advancer"]["output_dir"] = "D/" + config["advancer"]["output_id"] = "example_2D" + config["advancer"]["write_to"] = "txt" + + advancer = time_advancer(config) # Initial condition - mu1 = 1.0 - mu2 = 0.0 - mu3 = 0.1 + mu1 = 1.0 + mu2 = 0.0 + mu3 = 0.1 sigma1 = 0.1 sigma2 = 0.2 - advancer.initialize_state_gaussian_bivar( mu1, mu2, sigma1, sigma2 ) + advancer.initialize_state_gaussian_bivar(mu1, mu2, sigma1, sigma2) advancer.run() - - return - -if __name__ == '__main__': - np.set_printoptions( formatter = { 'float': '{: 0.4E}'.format } ) - nargs = len( sys.argv ) +if __name__ == "__main__": + + np.set_printoptions(formatter={"float": "{: 0.4E}".format}) + + nargs = len(sys.argv) if nargs == 2: config_file = sys.argv[1] - config_mgr = config_manager( config_file ) - config = config_mgr.get_config() - advance_example( config ) + config_mgr = config_manager(config_file) + config = config_mgr.get_config() + advance_example(config) ### Complete workflow for devel driver: ### 1. Check whether input file exists ### 2. If yes, run advance_example, then stop @@ -160,7 +257,10 @@ def advance_example2dp1(): ### 4. If argv matches case, run, then stop ### 5. If argv does not match case, then exit else: - print('devel_driver: no config file supplied') + flow_example_3d() + # flow_example_2d() + # flow_example_1d() + #print('devel_driver: no config file supplied') + exit() - diff --git a/src/inversion.py b/src/inversion.py index d4f74d8..2d6c3b3 100644 --- a/src/inversion.py +++ b/src/inversion.py @@ -1,10 +1,9 @@ import numpy as np -from scipy import sqrt, zeros -from scipy.linalg import eig import math -# from numba import jit -# @jit(nopython=True) +from numba import njit + +@njit def sign(q): if q > 0: return 1 @@ -13,9 +12,10 @@ def sign(q): else: return -1 -def wheeler(moments, adaptive = False): + +def wheeler(moments, adaptive=False): """ - This function inverts moments into 1D quadrature weights and abscissas using adaptive Wheeler algorithm. + This function inverts moments into 1D quadrature weights and abscissas using adaptive Wheeler algorithm. :param moments: Statistical moments of the transported PDF :type moments: array like @@ -23,19 +23,19 @@ def wheeler(moments, adaptive = False): :rtype: array like """ - n = len( moments ) // 2 + n = len(moments) // 2 # Adaptivity parameters - rmax = 1e-8 - eabs = 1e-8 + rmax = 1e-8 + eabs = 1e-8 cutoff = 0 - + # Check if moments are unrealizable. if moments[0] <= 0: print("Wheeler: Moments are NOT realizable (moment[0] <= 0.0). Run failed.") exit() - if n == 1 or ( adaptive and moments[0] < rmax ): + if n == 1 or (adaptive and moments[0] < rmax): w = moments[0] x = moments[1] / moments[0] return x, w @@ -45,47 +45,57 @@ def wheeler(moments, adaptive = False): # Construct recurrence matrix ind = n - a = np.zeros( n ) - b = np.zeros( n ) - sigma = np.zeros( ( 2*ind+1, 2*ind+1 ) ) + a = np.zeros(n) + b = np.zeros(n) + sigma = np.zeros((2 * ind + 1, 2 * ind + 1)) - for i in range(1,2*ind+1): - sigma[1,i] = nu[i-1] + for i in range(1, 2 * ind + 1): + sigma[1, i] = nu[i - 1] - a[0] = nu[1]/nu[0] + a[0] = nu[1] / nu[0] b[0] = 0 - for k in range(2,ind+1): - for l in range(k,2*ind-k+2): - sigma[k, l] = sigma[k-1, l+1]-a[k-2]*sigma[k-1, l]-b[k-2]*sigma[k-2, l] - a[k-1] = sigma[k, k+1]/sigma[k, k]-sigma[k-1, k]/sigma[k-1, k-1] - b[k-1] = sigma[k, k]/sigma[k-1, k-1] + for k in range(2, ind + 1): + for l in range(k, 2 * ind - k + 2): + sigma[k, l] = ( + sigma[k - 1, l + 1] + - a[k - 2] * sigma[k - 1, l] + - b[k - 2] * sigma[k - 2, l] + ) + a[k - 1] = sigma[k, k + 1] / sigma[k, k] - sigma[k - 1, k] / sigma[k - 1, k - 1] + b[k - 1] = sigma[k, k] / sigma[k - 1, k - 1] # Find maximum n using diagonal element of sigma if adaptive: - for k in range(ind,1,-1): - if sigma[k,k] <= cutoff: - n = k-1 + for k in range(ind, 1, -1): + if sigma[k, k] <= cutoff: + n = k - 1 if n == 1: w = moments[0] - x = moments[1]/moments[0] + x = moments[1] / moments[0] return x, w # Use maximum n to re-calculate recurrence matrix - a = np.zeros( n ) - b = np.zeros( n ) - w = np.zeros( n ) - x = np.zeros( n ) - sigma = np.zeros( (2*n+1, 2*n+1) ) - sigma[1,1:] = nu - + a = np.zeros(n) + b = np.zeros(n) + w = np.zeros(n) + x = np.zeros(n) + sigma = np.zeros((2 * n + 1, 2 * n + 1)) + sigma[1, 1:] = nu + a[0] = nu[1] / nu[0] b[0] = 0 - for k in range(2, n+1): - for l in range(k, 2*n-k+2): - sigma[k, l] = sigma[k-1, l+1]-a[k-2]*sigma[k-1, l]-b[k-2]*sigma[k-2, l] - a[k-1] = sigma[k, k+1]/sigma[k, k]-sigma[k-1, k]/sigma[k-1, k-1] - b[k-1] = sigma[k, k]/sigma[k-1, k-1] + for k in range(2, n + 1): + for l in range(k, 2 * n - k + 2): + sigma[k, l] = ( + sigma[k - 1, l + 1] + - a[k - 2] * sigma[k - 1, l] + - b[k - 2] * sigma[k - 2, l] + ) + a[k - 1] = ( + sigma[k, k + 1] / sigma[k, k] - sigma[k - 1, k] / sigma[k - 1, k - 1] + ) + b[k - 1] = sigma[k, k] / sigma[k - 1, k - 1] # Check if moments are unrealizable if b.min() < 0: @@ -93,45 +103,46 @@ def wheeler(moments, adaptive = False): exit() # Setup Jacobi matrix for n-point quadrature, adapt n using rmin and eabs - for n1 in range(n,0,-1): + for n1 in range(n, 0, -1): if n1 == 1: w = moments[0] - x = moments[1]/moments[0] + x = moments[1] / moments[0] return x, w # Jacobi matrix - sqrt_b = np.sqrt( b[1:] ) - jacobi = np.diag( a ) + np.diag( sqrt_b, -1 ) + np.diag( sqrt_b, 1 ) - + sqrt_b = np.sqrt(b[1:]) + jacobi = np.diag(a) + np.diag(sqrt_b, -1) + np.diag(sqrt_b, 1) + # Compute weights and abscissas - eigenvalues, eigenvectors = np.linalg.eig( jacobi ) + eigenvalues, eigenvectors = np.linalg.eig(jacobi) idx = eigenvalues.argsort() - x = eigenvalues[idx].real - eigenvectors = eigenvectors[:,idx].real - w = moments[0]*eigenvectors[0,:]**2 + x = eigenvalues[idx].real + eigenvectors = eigenvectors[:, idx].real + w = moments[0] * eigenvectors[0, :] ** 2 # Adaptive conditions. When both satisfied, return the results. if adaptive: - dab = zeros(n1) - mab = zeros(n1) + dab = np.zeros(n1) + mab = np.zeros(n1) - for i in range(n1-1,0,-1): - dab[i] = min(abs(x[i]-x[0:i])) - mab[i] = max(abs(x[i]-x[0:i])) + for i in range(n1 - 1, 0, -1): + dab[i] = min(abs(x[i] - x[0:i])) + mab[i] = max(abs(x[i] - x[0:i])) mindab = min(dab[1:n1]) maxmab = max(mab[1:n1]) if n1 == 2: maxmab = 1 - if min(w)/max(w) > rmax and mindab/maxmab > eabs: + if min(w) / max(w) > rmax and mindab / maxmab > eabs: return x, w else: return x, w -def hyperbolic(moments, max_skewness = 30, checks = True): + +def hyperbolic(moments, max_skewness=30, checks=True): """ - This is a driver for hyperbolic qmom. + This is a driver for hyperbolic qmom. It calls :func:`hyqmom2` if ``len(moments) = 2``, or :func:`hyqmom3` if ``len(moments) = 3`` :param moments: Statistical moments of the transported PDF @@ -141,22 +152,23 @@ def hyperbolic(moments, max_skewness = 30, checks = True): :return: Abscissas, weights :rtype: array like """ - - num_moments = len( moments ) + + num_moments = len(moments) if num_moments == 3: - return hyqmom2( moments ) - elif num_moments == 5: - return hyqmom3( moments, max_skewness, checks ) + return hyqmom2(moments) + elif num_moments == 5: + return hyqmom3(moments, max_skewness, checks) else: - print('inversion: hyperbolic: incorrect number of moments(%i)' % num_moments) + print("inversion: hyperbolic: incorrect number of moments(%i)" % num_moments) return -# @jit(nopython=True) + +@njit def hyqmom2(moments): """ This function inverts moments into a two-node quadrature rule. - :param moments: Statistical moments of the transported PDF + :param moments: Statistical moments of the transported PDF :type moments: array like :return: Abscissas, weights :rtype: array like @@ -165,23 +177,24 @@ def hyqmom2(moments): n = 2 w = np.zeros(n) x = np.zeros(n) - - w[0] = moments[0]/2. + + w[0] = moments[0] / 2.0 w[1] = w[0] - bx = moments[1]/moments[0] - d2 = moments[2]/moments[0] - c2 = d2 - bx**2. + bx = moments[1] / moments[0] + d2 = moments[2] / moments[0] + c2 = d2 - bx ** 2.0 - if c2 < 10**(-12): - c2 = 10**(-12) - x[0] = bx - math.sqrt(c2) - x[1] = bx + math.sqrt(c2) + if c2 < 10 ** (-12): + c2 = 10 ** (-12) + x[0] = bx - np.sqrt(c2) + x[1] = bx + np.sqrt(c2) - return x,w + return x, w -# @jit(nopython=True) -def hyqmom3(moments, max_skewness = 30, checks = True): + +@njit +def hyqmom3(moments, max_skewness=30, checks=True): """ This function inverts moments into a three-node quadrature rule. @@ -192,111 +205,113 @@ def hyqmom3(moments, max_skewness = 30, checks = True): :return: Abscissas, weights :rtype: array like """ - + # print('moments in:', moments) n = 3 - etasmall = 10**(-10) - verysmall = 10**(-14) - realsmall = 10**(-14) + etasmall = 10 ** (-10) + verysmall = 10 ** (-14) + realsmall = 10 ** (-14) - w = np.zeros(n) - x = np.zeros(n) - xp = np.zeros(n) + w = np.zeros(n) + x = np.zeros(n) + xp = np.zeros(n) xps = np.zeros(n) rho = np.zeros(n) if moments[0] <= verysmall and checks: w[1] = moments[0] - return x,w + return x, w - bx = moments[1]/moments[0] - d2 = moments[2]/moments[0] - d3 = moments[3]/moments[0] - d4 = moments[4]/moments[0] - c2 = d2 - bx**2 - c3 = d3 - 3*bx*d2 + 2*bx**3 - c4 = d4 - 4*bx*d3 + 6*(bx**2)*d2 - 3*bx**4 + bx = moments[1] / moments[0] + d2 = moments[2] / moments[0] + d3 = moments[3] / moments[0] + d4 = moments[4] / moments[0] + c2 = d2 - bx ** 2 + c3 = d3 - 3 * bx * d2 + 2 * bx ** 3 + c4 = d4 - 4 * bx * d3 + 6 * (bx ** 2) * d2 - 3 * bx ** 4 if checks: - if c2 < 0: + if c2 < 0: if c2 < -verysmall: - print("Error: c2 negative in three node HYQMOM") - return + raise Exception("Error: c2 negative in three node HYQMOM") else: - realizable = c2*c4 - c2**3 - c3**2 + realizable = c2 * c4 - c2 ** 3 - c3 ** 2 if realizable < 0: if c2 >= etasmall: - q = c3/math.sqrt(c2)/c2 - eta = c4/c2/c2 + q = c3 / np.sqrt(c2) / c2 + eta = c4 / c2 / c2 if abs(q) > verysmall: - slope = (eta - 3)/q - det = 8 + slope**2 - qp = 0.5*(slope + math.sqrt(det)) - qm = 0.5*(slope - math.sqrt(det)) - if q > 0: + slope = (eta - 3) / q + det = 8 + slope ** 2 + qp = 0.5 * (slope + np.sqrt(det)) + qm = 0.5 * (slope - np.sqrt(det)) + if q > 0: q = qp else: q = qm else: q = 0 - eta = q**2 + 1 - c3 = q*math.sqrt(c2)*c2 - c4 = eta*c2**2 - if realizable < -10.**(-6): - print("Error: c4 small in HYQMOM3") - return - else: - c3 = 0. - c4 = c2**2. - - scale = math.sqrt(c2) + eta = q ** 2 + 1 + c3 = q * np.sqrt(c2) * c2 + c4 = eta * c2 ** 2 + if realizable < -(10.0 ** (-6)): + raise Exception("Error: c4 small in HYQMOM3") + else: + c3 = 0.0 + c4 = c2 ** 2.0 + + scale = np.sqrt(c2) if checks and c2 < etasmall: - q = 0 + q = 0 eta = 1 else: - q = c3/math.sqrt(c2)/c2 - eta = c4/c2/c2 + q = c3 / np.sqrt(c2) / c2 + eta = c4 / c2 / c2 - if q**2 > max_skewness**2: - slope = (eta - 3)/q + # print('q=',q) + # print('maxskew=',max_skewness) + if q ** 2 > max_skewness ** 2: + slope = (eta - 3) / q if q > 0: q = max_skewness else: q = -max_skewness - eta = 3 + slope*q + eta = 3 + slope * q if checks: - realizable = eta - 1 - q**2 + realizable = eta - 1 - q ** 2 if realizable < 0: - eta = 1 + q**2 + eta = 1 + q ** 2 - xps[0] = (q - math.sqrt(4*eta-3*q**2))/2. - xps[1] = 0. - xps[2] = (q + math.sqrt(4*eta-3*q**2))/2. + xps[0] = (q - np.sqrt(4 * eta - 3 * q ** 2)) / 2.0 + xps[1] = 0.0 + xps[2] = (q + np.sqrt(4 * eta - 3 * q ** 2)) / 2.0 - dem = 1./math.sqrt(4*eta - 3*q**2) - prod = -xps[0]*xps[2] - prod = max(prod,1+realsmall) + dem = 1.0 / np.sqrt(4 * eta - 3 * q ** 2) + prod = -xps[0] * xps[2] + prod = max(prod, 1 + realsmall) - rho[0] = -dem/xps[0] - rho[1] = 1 - 1/prod - rho[2] = dem/xps[2] + rho[0] = -dem / xps[0] + rho[1] = 1 - 1 / prod + rho[2] = dem / xps[2] srho = np.sum(rho) - rho = rho/srho - if min(rho) < 0 : - print("Error: Negative weight in HYQMOM") - return - scales = np.sum(rho*xps**2)/np.sum(rho) - for i in range(n): - xp[i] = xps[i]*scale/math.sqrt(scales) + rho = rho / srho + if min(rho) < 0: + raise Exception("Error: Negative weight in HYQMOM") - w = moments[0]*rho + scales = np.sum(rho * xps ** 2) / np.sum(rho) + xp = xps * scale / np.sqrt(scales) + + w = moments[0] * rho x = xp - x = bx + x + x = bx + x + + return x, w - return x,w -def conditional_hyperbolic(moments, indices, max_skewness = 30, checks = True): + +@njit +def conditional_hyperbolic(moments, indices, max_skewness=30, checks=True): """ This function inverts moments into a two-node quadrature rule. @@ -305,18 +320,18 @@ def conditional_hyperbolic(moments, indices, max_skewness = 30, checks = True): :return: Abscissas, weights :rtype: array like """ - - num_dim = len(indices) - if num_dim == 6: - return chyqmom4( moments, indices ) - elif num_dim == 10: - return chyqmom9( moments, indices, max_skewness, checks) - elif num_dim == 16: - return chyqmom27( moments, indices, max_skewness, checks ) + # num_dim = len(indices) + + if num_dim == 16: + return chyqmom27(moments, indices, max_skewness, checks) + if num_dim == 10: + return chyqmom9(moments, indices, max_skewness, checks) + # if num_dim == 6: + # return chyqmom4(moments, indices) -# @jit(nopython=True) -def chyqmom4(moments, indices, max_skewness = 30): +@njit +def chyqmom4(moments, indices, max_skewness=30): # normalidx = indices.tolist() # mom00 = moments[normalidx.index([0,0])] @@ -338,64 +353,53 @@ def chyqmom4(moments, indices, max_skewness = 30): x = np.zeros(n) y = np.zeros(n) - bx = mom10/mom00 - by = mom01/mom00 - d20 = mom20/mom00 - d11 = mom11/mom00 - d02 = mom02/mom00 + bx = mom10 / mom00 + by = mom01 / mom00 + d20 = mom20 / mom00 + d11 = mom11 / mom00 + d02 = mom02 / mom00 - c20 = d20 - bx**2. - c11 = d11 - bx*by - c02 = d02 - by**2. + c20 = d20 - bx ** 2.0 + c11 = d11 - bx * by + c02 = d02 - by ** 2.0 M1 = np.array([1, 0, c20]) - xp, rho = hyqmom2(M1) - yf = c11*xp/c20 - mu2avg = c02 - np.sum(rho*yf**2) - mu2avg = max(mu2avg,0) - mu2 = mu2avg - M3 = np.array([ 1, 0, mu2 ]) + xp, rho = hyqmom2(M1) + yf = c11 * xp / c20 + mu2avg = c02 - np.sum(rho * yf ** 2) + mu2avg = max(mu2avg, 0) + mu2 = mu2avg + M3 = np.array([1, 0, mu2]) xp3, rh3 = hyqmom2(M3) - yp21 = xp3[0] - yp22 = xp3[1] - rho21 = rh3[0] - rho22 = rh3[1] - - w[0] = rho[0]*rho21 - w[1] = rho[0]*rho22 - w[2] = rho[1]*rho21 - w[3] = rho[1]*rho22 - w = mom00*w - - x[0] = xp[0] - x[1] = xp[0] - x[2] = xp[1] - x[3] = xp[1] - x = bx + x - - y[0] = yf[0] + yp21 - y[1] = yf[0] + yp22 - y[2] = yf[1] + yp21 - y[3] = yf[1] + yp22 - y = by + y - - x = [ x, y ] - return x,w - -# @jit(nopython=True) -def chyqmom9(moments, indices, max_skewness = 30, checks = True): + yp21 = xp3[0] + yp22 = xp3[1] + rho21 = rh3[0] + rho22 = rh3[1] + + w[0] = rho[0] * rho21 + w[1] = rho[0] * rho22 + w[2] = rho[1] * rho21 + w[3] = rho[1] * rho22 + w = mom00 * w + + x[0] = xp[0] + x[1] = xp[0] + x[2] = xp[1] + x[3] = xp[1] + x = bx + x - # normalidx = indices.tolist() - # mom00 = moments[normalidx.index([0,0])] - # mom10 = moments[normalidx.index([1,0])] - # mom01 = moments[normalidx.index([0,1])] - # mom20 = moments[normalidx.index([2,0])] - # mom11 = moments[normalidx.index([1,1])] - # mom02 = moments[normalidx.index([0,2])] - # mom30 = moments[normalidx.index([3,0])] - # mom03 = moments[normalidx.index([0,3])] - # mom40 = moments[normalidx.index([4,0])] - # mom04 = moments[normalidx.index([0,4])] + y[0] = yf[0] + yp21 + y[1] = yf[0] + yp22 + y[2] = yf[1] + yp21 + y[3] = yf[1] + yp22 + y = by + y + + x = [x, y] + return x, w + + +@njit +def chyqmom9(moments, indices, max_skewness=30, checks=True): mom00 = moments[0] mom10 = moments[1] @@ -413,141 +417,136 @@ def chyqmom9(moments, indices, max_skewness = 30, checks = True): x = np.zeros(n) y = np.zeros(n) - csmall = 10.**(-10) - verysmall = 10.**(-14) - - bx = mom10/mom00 - by = mom01/mom00 - d20 = mom20/mom00 - d11 = mom11/mom00 - d02 = mom02/mom00 - d30 = mom30/mom00 - d03 = mom03/mom00 - d40 = mom40/mom00 - d04 = mom04/mom00 - - c20 = d20 - bx**2. - c11 = d11 - bx*by - c02 = d02 - by**2. - c30 = d30 - 3.*bx*d20 + 2.*bx**3. - c03 = d03 - 3.*by*d02 + 2.*by**3. - c40 = d40 - 4.*bx*d30 + 6*(bx**2.)*d20 - 3.*bx**(4) - c04 = d04 - 4.*by*d03 + 6*(by**2.)*d02 - 3.*by**(4) + abscissas = np.zeros((n, 2)) + + csmall = 10.0 ** (-10) + verysmall = 10.0 ** (-14) + + if mom00 < verysmall and checks: + abscissas[:, 0] = x[:] + abscissas[:, 1] = y[:] + w[3] = mom00 + return abscissas, w + + bx = mom10 / mom00 + by = mom01 / mom00 + d20 = mom20 / mom00 + d11 = mom11 / mom00 + d02 = mom02 / mom00 + d30 = mom30 / mom00 + d03 = mom03 / mom00 + d40 = mom40 / mom00 + d04 = mom04 / mom00 + + c20 = d20 - bx ** 2.0 + c11 = d11 - bx * by + c02 = d02 - by ** 2.0 + c30 = d30 - 3.0 * bx * d20 + 2.0 * bx ** 3.0 + c03 = d03 - 3.0 * by * d02 + 2.0 * by ** 3.0 + c40 = d40 - 4.0 * bx * d30 + 6 * (bx ** 2.0) * d20 - 3.0 * bx ** (4) + c04 = d04 - 4.0 * by * d03 + 6 * (by ** 2.0) * d02 - 3.0 * by ** (4) M1 = np.array([1, 0, c20, c30, c40]) - xp, rho = hyqmom3(M1,max_skewness,checks) + xp, rho = hyqmom3(M1, max_skewness, checks) if checks and c20 < csmall: - rho[0] = 0. - rho[1] = 1. - rho[2] = 0. - yf = 0*xp + rho[0] = 0.0 + rho[1] = 1.0 + rho[2] = 0.0 + yf = 0 * xp M2 = np.array([1, 0, c02, c03, c04]) - xp2, rho2 = hyqmom3(M2,max_skewness,checks) - yp21 = xp2[0] - yp22 = xp2[1] - yp23 = xp2[2] - rho21 = rho2[0] - rho22 = rho2[1] + xp2, rho2 = hyqmom3(M2, max_skewness, checks) + yp21 = xp2[0] + yp22 = xp2[1] + yp23 = xp2[2] + rho21 = rho2[0] + rho22 = rho2[1] rho23 = rho2[2] else: - yf = c11*xp/c20 - mu2avg = c02 - np.sum(rho*(yf**2.)) - mu2avg = max(mu2avg, 0.) - mu2 = mu2avg - mu3 = 0*mu2 - mu4 = mu2**2. + yf = c11 * xp / c20 + mu2avg = c02 - np.sum(rho * (yf ** 2.0)) + mu2avg = max(mu2avg, 0.0) + mu2 = mu2avg + mu3 = 0 * mu2 + mu4 = mu2 ** 2.0 if mu2 > csmall: - q = (c03 - np.sum(rho*(yf**3.)))/mu2**(3./2.) - eta = (c04 - np.sum(rho*(yf**4.)) - 6*np.sum(rho*(yf**2.))*mu2)/mu2**2. - if eta < (q**2 + 1): + q = (c03 - np.sum(rho * (yf ** 3.0))) / mu2 ** (3.0 / 2.0) + eta = ( + c04 - np.sum(rho * (yf ** 4.0)) - 6 * np.sum(rho * (yf ** 2.0)) * mu2 + ) / mu2 ** 2.0 + if eta < (q ** 2 + 1): if abs(q) > verysmall: - slope = (eta - 3.)/q - det = 8. + slope**2. - qp = 0.5*(slope + math.sqrt(det)) - qm = 0.5*(slope - math.sqrt(det)) + slope = (eta - 3.0) / q + det = 8.0 + slope ** 2.0 + qp = 0.5 * (slope + np.sqrt(det)) + qm = 0.5 * (slope - np.sqrt(det)) if q > 0: q = qp else: q = qm - else: + else: q = 0 - eta = q**2 + 1 + eta = q ** 2 + 1 - mu3 = q*mu2**(3./2.) - mu4 = eta*mu2**2. + mu3 = q * mu2 ** (3.0 / 2.0) + mu4 = eta * mu2 ** 2.0 M3 = np.array([1, 0, mu2, mu3, mu4]) - xp3, rh3 = hyqmom3(M3,max_skewness,checks) - yp21 = xp3[0] - yp22 = xp3[1] - yp23 = xp3[2] - rho21 = rh3[0] - rho22 = rh3[1] + xp3, rh3 = hyqmom3(M3, max_skewness, checks) + yp21 = xp3[0] + yp22 = xp3[1] + yp23 = xp3[2] + rho21 = rh3[0] + rho22 = rh3[1] rho23 = rh3[2] - w[0] = rho[0]*rho21 - w[1] = rho[0]*rho22 - w[2] = rho[0]*rho23 - w[3] = rho[1]*rho21 - w[4] = rho[1]*rho22 - w[5] = rho[1]*rho23 - w[6] = rho[2]*rho21 - w[7] = rho[2]*rho22 - w[8] = rho[2]*rho23 - w = mom00*w - - x[0] = xp[0] + w[0] = rho[0] * rho21 + w[1] = rho[0] * rho22 + w[2] = rho[0] * rho23 + w[3] = rho[1] * rho21 + w[4] = rho[1] * rho22 + w[5] = rho[1] * rho23 + w[6] = rho[2] * rho21 + w[7] = rho[2] * rho22 + w[8] = rho[2] * rho23 + w = mom00 * w + + x[0] = xp[0] x[1] = xp[0] - x[2] = xp[0] - x[3] = xp[1] + x[2] = xp[0] + x[3] = xp[1] x[4] = xp[1] x[5] = xp[1] - x[6] = xp[2] + x[6] = xp[2] x[7] = xp[2] - x[8] = xp[2] - x = bx + x + x[8] = xp[2] + x = bx + x - y[0] = yf[0] + yp21 - y[1] = yf[0] + yp22 - y[2] = yf[0] + yp23 - y[3] = yf[1] + yp21 + y[0] = yf[0] + yp21 + y[1] = yf[0] + yp22 + y[2] = yf[0] + yp23 + y[3] = yf[1] + yp21 y[4] = yf[1] + yp22 - y[5] = yf[1] + yp23 + y[5] = yf[1] + yp23 y[6] = yf[2] + yp21 - y[7] = yf[2] + yp22 - y[8] = yf[2] + yp23 - y = by + y + y[7] = yf[2] + yp22 + y[8] = yf[2] + yp23 + y = by + y - x = [ x, y ] - return x,w + abscissas[:, 0] = x[:] + abscissas[:, 1] = y[:] + # SHB Note: If running 0D case then the abscissas output might not work -# @jit(nopython=True) -def chyqmom27(moments, indices, max_skewness = 30, checks = True): + return abscissas, w - # Indices used for calling chyqmom9 - RF_idx = np.array( [ - [0,0], [1,0], [0,1], [2,0], [1,1], - [0,2], [3,0], [0,3], [4,0], [0,4] ] - ) - # normalidx = indices.tolist() - # m000 = moments[normalidx.index([0,0,0])] - # m100 = moments[normalidx.index([1,0,0])] - # m010 = moments[normalidx.index([0,1,0])] - # m001 = moments[normalidx.index([0,0,1])] - # m200 = moments[normalidx.index([2,0,0])] - # m110 = moments[normalidx.index([1,1,0])] - # m101 = moments[normalidx.index([1,0,1])] - # m020 = moments[normalidx.index([0,2,0])] - # m011 = moments[normalidx.index([0,1,1])] - # m002 = moments[normalidx.index([0,0,2])] - # m300 = moments[normalidx.index([3,0,0])] - # m030 = moments[normalidx.index([0,3,0])] - # m003 = moments[normalidx.index([0,0,3])] - # m400 = moments[normalidx.index([4,0,0])] - # m040 = moments[normalidx.index([0,4,0])] - # m004 = moments[normalidx.index([0,0,4])] +@njit +def chyqmom27(moments, indices, max_skewness=30, checks=True): + + # Indices used for calling chyqmom9 + RF_idx = np.array( + [[0, 0], [1, 0], [0, 1], [2, 0], [1, 1], [0, 2], [3, 0], [0, 3], [4, 0], [0, 4]] + ) m000 = moments[0] m100 = moments[1] @@ -566,90 +565,91 @@ def chyqmom27(moments, indices, max_skewness = 30, checks = True): m040 = moments[14] m004 = moments[15] - small = 1.e-10 - isosmall = 1.e-14 - csmall = 1.e-10 - wsmall = 1.e-4 - verysmall = 1.e-14 + small = 1.0e-10 + isosmall = 1.0e-14 + csmall = 1.0e-10 + wsmall = 1.0e-4 + verysmall = 1.0e-14 n = 27 w = np.zeros(n) - x = np.zeros(n) - y = np.zeros(n) - z = np.zeros(n) + abscissas = np.zeros((n, 3)) + Yf = np.zeros(3) + Zf = np.zeros((3, 3)) + W = np.zeros(n) if m000 <= verysmall and checks: - w[12] = m000 - return - - bx = m100/m000 - by = m010/m000 - bz = m001/m000 - - if checks and m000 <= isosmall: - d200 = m200/m000 - d020 = m020/m000 - d002 = m002/m000 - d300 = m300/m000 - d030 = m030/m000 - d003 = m003/m000 - d400 = m400/m000 - d040 = m040/m000 - d004 = m004/m000 - - c200 = d200 - bx**2 - c020 = d020 - by**2 - c002 = d002 - bz**2 - c300 = d300 - 3*bx*d200 + 2*bx**3 - c030 = d030 - 3*by*d020 + 2*by**3 - c003 = d003 - 3*bz*d002 + 2*bz**3 - c400 = d400 - 4*bx*d300 + 6*(bx**2)*d200 - 3*bx**4 - c040 = d040 - 4*by*d030 + 6*(by**2)*d020 - 3*by**4 - c004 = d004 - 4*bz*d003 + 6*(bz**2)*d002 - 3*bz**4 + W[12] = m000 + return abscissas, W + + bx = m100 / m000 + by = m010 / m000 + bz = m001 / m000 + + if checks and m000 <= isosmall: + d200 = m200 / m000 + d020 = m020 / m000 + d002 = m002 / m000 + d300 = m300 / m000 + d030 = m030 / m000 + d003 = m003 / m000 + d400 = m400 / m000 + d040 = m040 / m000 + d004 = m004 / m000 + + c200 = d200 - bx ** 2 + c020 = d020 - by ** 2 + c002 = d002 - bz ** 2 + c300 = d300 - 3 * bx * d200 + 2 * bx ** 3 + c030 = d030 - 3 * by * d020 + 2 * by ** 3 + c003 = d003 - 3 * bz * d002 + 2 * bz ** 3 + c400 = d400 - 4 * bx * d300 + 6 * (bx ** 2) * d200 - 3 * bx ** 4 + c040 = d040 - 4 * by * d030 + 6 * (by ** 2) * d020 - 3 * by ** 4 + c004 = d004 - 4 * bz * d003 + 6 * (bz ** 2) * d002 - 3 * bz ** 4 c110 = 0 - c101 = 0 + c101 = 0 c011 = 0 else: - d200 = m200/m000 - d110 = m110/m000 - d101 = m101/m000 - d020 = m020/m000 - d011 = m011/m000 - d002 = m002/m000 - d300 = m300/m000 - d030 = m030/m000 - d003 = m003/m000 - d400 = m400/m000 - d040 = m040/m000 - d004 = m004/m000 - - c200 = d200 - bx**2 - c110 = d110 - bx*by - c101 = d101 - bx*bz - c020 = d020 - by**2 - c011 = d011 - by*bz - c002 = d002 - bz**2 - c300 = d300 - 3*bx*d200 + 2*bx**3 - c030 = d030 - 3*by*d020 + 2*by**3 - c003 = d003 - 3*bz*d002 + 2*bz**3 - c400 = d400 - 4*bx*d300 + 6*bx**2*d200 - 3*bx**4 - c040 = d040 - 4*by*d030 + 6*by**2*d020 - 3*by**4 - c004 = d004 - 4*bz*d003 + 6*bz**2*d002 - 3*bz**4 + d200 = m200 / m000 + d110 = m110 / m000 + d101 = m101 / m000 + d020 = m020 / m000 + d011 = m011 / m000 + d002 = m002 / m000 + d300 = m300 / m000 + d030 = m030 / m000 + d003 = m003 / m000 + d400 = m400 / m000 + d040 = m040 / m000 + d004 = m004 / m000 + + c200 = d200 - bx ** 2 + c110 = d110 - bx * by + c101 = d101 - bx * bz + c020 = d020 - by ** 2 + c011 = d011 - by * bz + c002 = d002 - bz ** 2 + c300 = d300 - 3 * bx * d200 + 2 * bx ** 3 + c030 = d030 - 3 * by * d020 + 2 * by ** 3 + c003 = d003 - 3 * bz * d002 + 2 * bz ** 3 + c400 = d400 - 4 * bx * d300 + 6 * bx ** 2 * d200 - 3 * bx ** 4 + c040 = d040 - 4 * by * d030 + 6 * by ** 2 * d020 - 3 * by ** 4 + c004 = d004 - 4 * bz * d003 + 6 * bz ** 2 * d002 - 3 * bz ** 4 if c200 <= 0 and checks: c200 = 0 c300 = 0 c400 = 0 - if c200*c400 < (c200**3 + c300**2) and checks: - q = c300/c200**(3./2.) - eta = c400/c200**2 + if c200 * c400 < (c200 ** 3 + c300 ** 2) and checks: + q = c300 / c200 ** (3.0 / 2.0) + eta = c400 / c200 ** 2 if abs(q) > verysmall: - slope = (eta - 3.)/q - det = 8 + slope**2 - qp = 0.5*( slope + math.sqrt(det) ) - qm = 0.5*( slope - math.sqrt(det) ) + slope = (eta - 3.0) / q + det = 8 + slope ** 2 + qp = 0.5 * (slope + np.sqrt(det)) + qm = 0.5 * (slope - np.sqrt(det)) if q > 0: q = qp else: @@ -657,604 +657,572 @@ def chyqmom27(moments, indices, max_skewness = 30, checks = True): else: q = 0 - eta = q**2 + 1 - c300 = q*c200**(3./2.) - c400 = eta*c200**2. + eta = q ** 2 + 1 + c300 = q * c200 ** (3.0 / 2.0) + c400 = eta * c200 ** 2.0 if c020 <= 0 and checks: c020 = 0 c030 = 0 c040 = 0 - if c200*c400 < (c200**3 + c300**2) and checks: - q = c300/c200**(3/2) - eta = c400/c200**2 + if c200 * c400 < (c200 ** 3 + c300 ** 2) and checks: + q = c300 / c200 ** (3 / 2) + eta = c400 / c200 ** 2 if abs(q) > verysmall: - slope = (eta - 3)/q - det = 8 + slope**2 - qp = 0.5*( slope + math.sqrt(det) ) - qm = 0.5*( slope - math.sqrt(det) ) + slope = (eta - 3) / q + det = 8 + slope ** 2 + qp = 0.5 * (slope + np.sqrt(det)) + qm = 0.5 * (slope - np.sqrt(det)) if sign(q) == 1: - q = qp + q = qp else: - q = qm + q = qm else: - q = 0 - eta = q**2 + 1 - c300 = q*c200**(3/2) - c400 = eta*c200**2 + q = 0 + eta = q ** 2 + 1 + c300 = q * c200 ** (3 / 2) + c400 = eta * c200 ** 2 if c020 <= 0 and checks: - c020 = 0 - c030 = 0 - c040 = 0 + c020 = 0 + c030 = 0 + c040 = 0 - if c020*c040 < (c020**3 + c030**2) and checks: - q = c030/c020**(3/2) - eta = c040/c020**2 + if c020 * c040 < (c020 ** 3 + c030 ** 2) and checks: + q = c030 / c020 ** (3 / 2) + eta = c040 / c020 ** 2 if abs(q) > verysmall: - slope = (eta - 3)/q - det = 8 + slope**2 - qp = 0.5*( slope + math.sqrt(det) ) - qm = 0.5*( slope - math.sqrt(det) ) + slope = (eta - 3) / q + det = 8 + slope ** 2 + qp = 0.5 * (slope + np.sqrt(det)) + qm = 0.5 * (slope - np.sqrt(det)) if sign(q) == 1: - q = qp + q = qp else: - q = qm + q = qm else: - q = 0 - eta = q**2 + 1 - c030 = q*c020**(3/2) - c040 = eta*c020**2 + q = 0 + eta = q ** 2 + 1 + c030 = q * c020 ** (3 / 2) + c040 = eta * c020 ** 2 if c002 <= 0 and checks: - c002 = 0 - c003 = 0 - c004 = 0 + c002 = 0 + c003 = 0 + c004 = 0 - if c002*c004 < (c002**3 + c003**2) and checks: - q = c003/c002**(3/2) - eta = c004/c002**2 + if c002 * c004 < (c002 ** 3 + c003 ** 2) and checks: + q = c003 / c002 ** (3 / 2) + eta = c004 / c002 ** 2 if abs(q) > verysmall: - slope = (eta - 3)/q - det = 8 + slope**2 - qp = 0.5*( slope + math.sqrt(det) ) - qm = 0.5*( slope - math.sqrt(det) ) + slope = (eta - 3) / q + det = 8 + slope ** 2 + qp = 0.5 * (slope + np.sqrt(det)) + qm = 0.5 * (slope - np.sqrt(det)) if sign(q) == 1: - q = qp + q = qp else: - q = qm + q = qm else: - q = 0 - eta = q**2 + 1 - c003 = q*c002**(3/2) - c004 = eta*c002**2 - - M1 = np.array([ 1, 0, c200, c300, c400 ]) - xp, rho = hyqmom3(M1,max_skewness,checks) - - rho11 = 0 - rho12 = 1 - rho13 = 0 - rho21 = 0 - rho23 = 0 - rho31 = 0 - rho32 = 1 - rho33 = 0 - - yp11 = 0 - yp12 = 0 - yp13 = 0 - yp21 = 0 - yp22 = 0 - yp23 = 0 - yp31 = 0 - yp32 = 0 - yp33 = 0 - - Yf = np.zeros(3) - - rho111 = 0 - rho112 = 1 - rho113 = 0 - rho121 = 0 - rho122 = 1 - rho123 = 0 - rho131 = 0 - rho132 = 1 - rho133 = 0 - rho211 = 0 - rho212 = 1 - rho213 = 0 - rho221 = 0 - rho222 = 1 - rho223 = 0 - rho231 = 0 - rho232 = 1 - rho233 = 0 - rho311 = 0 - rho312 = 1 - rho313 = 0 - rho321 = 0 - rho322 = 1 - rho323 = 0 - rho331 = 0 - rho332 = 1 - rho333 = 0 - - zp111 = 0 - zp112 = 0 - zp113 = 0 - zp121 = 0 - zp122 = 0 - zp123 = 0 - zp131 = 0 - zp132 = 0 - zp133 = 0 - zp211 = 0 - zp212 = 0 - zp213 = 0 - zp221 = 0 - zp222 = 0 - zp223 = 0 - zp231 = 0 - zp232 = 0 - zp233 = 0 - zp311 = 0 - zp312 = 0 - zp313 = 0 - zp321 = 0 - zp322 = 0 - zp323 = 0 - zp331 = 0 - zp332 = 0 - zp333 = 0 - - Zf = np.zeros((3,3)) + q = 0 + eta = q ** 2 + 1 + c003 = q * c002 ** (3 / 2) + c004 = eta * c002 ** 2 + + M1 = np.array([1, 0, c200, c300, c400]) + xp, rho = hyqmom3(M1, max_skewness, checks) + + rho11 = 0 + rho12 = 1 + rho13 = 0 + rho21 = 0 + rho23 = 0 + rho31 = 0 + rho32 = 1 + rho33 = 0 + + yp11 = 0 + yp12 = 0 + yp13 = 0 + yp21 = 0 + yp22 = 0 + yp23 = 0 + yp31 = 0 + yp32 = 0 + yp33 = 0 + + rho111 = 0 + rho112 = 1 + rho113 = 0 + rho121 = 0 + rho122 = 1 + rho123 = 0 + rho131 = 0 + rho132 = 1 + rho133 = 0 + rho211 = 0 + rho212 = 1 + rho213 = 0 + rho221 = 0 + rho222 = 1 + rho223 = 0 + rho231 = 0 + rho232 = 1 + rho233 = 0 + rho311 = 0 + rho312 = 1 + rho313 = 0 + rho321 = 0 + rho322 = 1 + rho323 = 0 + rho331 = 0 + rho332 = 1 + rho333 = 0 + + zp111 = 0 + zp112 = 0 + zp113 = 0 + zp121 = 0 + zp122 = 0 + zp123 = 0 + zp131 = 0 + zp132 = 0 + zp133 = 0 + zp211 = 0 + zp212 = 0 + zp213 = 0 + zp221 = 0 + zp222 = 0 + zp223 = 0 + zp231 = 0 + zp232 = 0 + zp233 = 0 + zp311 = 0 + zp312 = 0 + zp313 = 0 + zp321 = 0 + zp322 = 0 + zp323 = 0 + zp331 = 0 + zp332 = 0 + zp333 = 0 if c200 <= csmall and checks: - if c020 <= csmall: - M0 = np.array([ 1, 0, c002, c003, c004]) - Z0, W0 = hyqmom3(M0,max_skewness,checks) - - rho[0] = 0 - rho[1] = 1 - rho[2] = 0 - rho22 = 1 - rho221 = W0[0] - rho222 = W0[1] - rho223 = W0[2] - xp = 0*xp - zp221 = Z0[0] - zp222 = Z0[1] - zp223 = Z0[2] + if c020 <= csmall: + M0 = np.array([1, 0, c002, c003, c004]) + Z0, W0 = hyqmom3(M0, max_skewness, checks) + + rho[0] = 0 + rho[1] = 1 + rho[2] = 0 + rho22 = 1 + rho221 = W0[0] + rho222 = W0[1] + rho223 = W0[2] + xp = 0 * xp + zp221 = Z0[0] + zp222 = Z0[1] + zp223 = Z0[2] else: - M1 = np.array([ - 1, 0, 0, - c020, c011, c002, - c030, c003, - c040, c004 - ]) - Q1, W1 = chyqmom9(M1,RF_idx,max_skewness,checks) - Y1 = Q1[0] - Z1 = Q1[1] - - rho[0] = 0 - rho[1] = 1 - rho[2] = 0 - rho12 = 0 - rho21 = 1 - rho22 = 1 - rho23 = 1 - rho31 = 0 - rho211 = W1[0] - rho212 = W1[1] - rho213 = W1[2] - rho221 = W1[3] - rho222 = W1[4] - rho223 = W1[5] - rho231 = W1[6] - rho232 = W1[7] - rho233 = W1[8] - - xp = 0*xp - yp21 = Y1[0] - yp22 = Y1[4] - yp23 = Y1[8] - zp211 = Z1[0] - zp212 = Z1[1] - zp213 = Z1[2] - zp221 = Z1[3] - zp222 = Z1[4] - zp223 = Z1[5] - zp231 = Z1[6] - zp232 = Z1[7] - zp233 = Z1[8] + M1 = np.array([1, 0, 0, c020, c011, c002, c030, c003, c040, c004]) + Q1, W1 = chyqmom9(M1, RF_idx, max_skewness, checks) + Y1 = Q1[:,0] + Z1 = Q1[:,1] + + rho[0] = 0 + rho[1] = 1 + rho[2] = 0 + rho12 = 0 + rho21 = 1 + rho22 = 1 + rho23 = 1 + rho31 = 0 + rho211 = W1[0] + rho212 = W1[1] + rho213 = W1[2] + rho221 = W1[3] + rho222 = W1[4] + rho223 = W1[5] + rho231 = W1[6] + rho232 = W1[7] + rho233 = W1[8] + + xp = 0 * xp + yp21 = Y1[0] + yp22 = Y1[4] + yp23 = Y1[8] + zp211 = Z1[0] + zp212 = Z1[1] + zp213 = Z1[2] + zp221 = Z1[3] + zp222 = Z1[4] + zp223 = Z1[5] + zp231 = Z1[6] + zp232 = Z1[7] + zp233 = Z1[8] elif c020 <= csmall and checks: - M2 = np.array([ - 1, 0, 0, - c200, c101, c002, - c300, c003, - c400, c004 - ]) - Q2, W2 = chyqmom9(M2,RF_idx,max_skewness,checks) - X2 = Q2[0] - Z2 = Q2[1] - - rho[0] = 1 - rho[1] = 1 - rho[2] = 1 - rho12 = 1 - rho22 = 1 - rho32 = 1 - rho121 = W2[0] - rho122 = W2[1] - rho123 = W2[2] - rho221 = W2[3] - rho222 = W2[4] - rho223 = W2[5] - rho321 = W2[6] - rho322 = W2[7] - rho323 = W2[8] - xp[0] = X2[0] - xp[1] = X2[4] - xp[2] = X2[8] - zp121 = Z2[0] - zp122 = Z2[1] - zp123 = Z2[2] - zp221 = Z2[3] - zp222 = Z2[4] - zp223 = Z2[5] - zp321 = Z2[6] - zp322 = Z2[7] - zp323 = Z2[8] + M2 = np.array([1, 0, 0, c200, c101, c002, c300, c003, c400, c004]) + Q2, W2 = chyqmom9(M2, RF_idx, max_skewness, checks) + X2 = Q2[:,0] + Z2 = Q2[:,1] + + rho[0] = 1 + rho[1] = 1 + rho[2] = 1 + rho12 = 1 + rho22 = 1 + rho32 = 1 + rho121 = W2[0] + rho122 = W2[1] + rho123 = W2[2] + rho221 = W2[3] + rho222 = W2[4] + rho223 = W2[5] + rho321 = W2[6] + rho322 = W2[7] + rho323 = W2[8] + xp[0] = X2[0] + xp[1] = X2[4] + xp[2] = X2[8] + zp121 = Z2[0] + zp122 = Z2[1] + zp123 = Z2[2] + zp221 = Z2[3] + zp222 = Z2[4] + zp223 = Z2[5] + zp321 = Z2[6] + zp322 = Z2[7] + zp323 = Z2[8] elif c002 <= csmall and checks: - M3 = np.array([ - 1, 0, 0, - c200, c110, c020, - c300, c030, - c400, c040 - ]) - Q3, W3 = chyqmom9(M3,RF_idx,max_skewness,checks) - X3 = Q3[0] - Y3 = Q3[1] - - rho[0] = 1 - rho[1] = 1 - rho[2] = 1 - rho11 = W3[0] - rho12 = W3[1] - rho13 = W3[2] - rho21 = W3[3] - rho22 = W3[4] - rho23 = W3[5] - rho31 = W3[6] - rho32 = W3[7] + M3 = np.array([1, 0, 0, c200, c110, c020, c300, c030, c400, c040]) + Q3, W3 = chyqmom9(M3, RF_idx, max_skewness, checks) + X3 = Q3[:,0] + Y3 = Q3[:,1] + + rho[0] = 1 + rho[1] = 1 + rho[2] = 1 + rho11 = W3[0] + rho12 = W3[1] + rho13 = W3[2] + rho21 = W3[3] + rho22 = W3[4] + rho23 = W3[5] + rho31 = W3[6] + rho32 = W3[7] rho33 = W3[8] - xp[0] = X3[0] - xp[1] = X3[4] - xp[2] = X3[8] - yp11 = Y3[0] - yp12 = Y3[1] - yp13 = Y3[2] - yp21 = Y3[3] - yp22 = Y3[4] - yp23 = Y3[5] - yp31 = Y3[6] - yp32 = Y3[7] - yp33 = Y3[8] + xp[0] = X3[0] + xp[1] = X3[4] + xp[2] = X3[8] + yp11 = Y3[0] + yp12 = Y3[1] + yp13 = Y3[2] + yp21 = Y3[3] + yp22 = Y3[4] + yp23 = Y3[5] + yp31 = Y3[6] + yp32 = Y3[7] + yp33 = Y3[8] else: - M4 = np.array([ - 1, 0, 0, - c200, c110, c020, - c300, c030, - c400, c040]) - Q4, W4 = chyqmom9(M4,RF_idx,max_skewness,checks) - X4 = Q4[0] - Y4 = Q4[1] - - rho11 = W4[0]/(W4[0] + W4[1] + W4[2]) - rho12 = W4[1]/(W4[0] + W4[1] + W4[2]) - rho13 = W4[2]/(W4[0] + W4[1] + W4[2]) - rho21 = W4[3]/(W4[3] + W4[4] + W4[5]) - rho22 = W4[4]/(W4[3] + W4[4] + W4[5]) - rho23 = W4[5]/(W4[3] + W4[4] + W4[5]) - rho31 = W4[6]/(W4[6] + W4[7] + W4[8]) - rho32 = W4[7]/(W4[6] + W4[7] + W4[8]) - rho33 = W4[8]/(W4[6] + W4[7] + W4[8]) - - Yf[0] = rho11*Y4[0] + rho12*Y4[1] + rho13*Y4[2] - Yf[1] = rho21*Y4[3] + rho22*Y4[4] + rho23*Y4[5] - Yf[2] = rho31*Y4[6] + rho32*Y4[7] + rho33*Y4[8] - - yp11 = Y4[0] - Yf[0] - yp12 = Y4[1] - Yf[0] - yp13 = Y4[2] - Yf[0] - yp21 = Y4[3] - Yf[1] - yp22 = Y4[4] - Yf[1] - yp23 = Y4[5] - Yf[1] - yp31 = Y4[6] - Yf[2] - yp32 = Y4[7] - Yf[2] - yp33 = Y4[8] - Yf[2] - scale1 = math.sqrt(c200) - scale2 = math.sqrt(c020) - Rho1 = np.diag(rho) - Rho2 = np.array([ - [rho11, rho12, rho13], - [rho21, rho22, rho23], - [rho31, rho32, rho33] - ]) - Yp2 = np.array([ - [yp11, yp12, yp13], - [yp21, yp22, yp23], - [yp31, yp32, yp33] - ]) - Yp2s = Yp2/scale2 - RAB = Rho1*Rho2 - XAB = np.array([ - [xp[0], xp[1], xp[2]], - [xp[0], xp[1], xp[2]], - [xp[0], xp[1], xp[2]] - ]) - XABs = XAB/scale1 - YAB = Yp2 + np.diag(Yf)*np.ones(3) - YABs = YAB/scale2 - C01 = np.multiply(RAB,YABs) - Yc0 = np.ones(3) - Yc1 = XABs - Yc2 = Yp2s - A1 = np.sum(np.multiply(C01,Yc1)) - A2 = np.sum(np.multiply(C01,Yc2)) - - c101s = c101/scale1 - c011s = c011/scale2 - if c101s**2 >= c002*(1 - small) : - c101s = sign(c101s)*math.sqrt(c002) - elif c011s**2 >= c002*(1 - small) : - c110s = c110/scale1/scale2 - c011s = sign(c011s)*math.sqrt(c002) - c101s = c110s*c011s + M4 = np.array([1, 0, 0, c200, c110, c020, c300, c030, c400, c040]) + Q4, W4 = chyqmom9(M4, RF_idx, max_skewness, checks) + X4 = Q4[:,0] + Y4 = Q4[:,1] + + rho11 = W4[0] / (W4[0] + W4[1] + W4[2]) + rho12 = W4[1] / (W4[0] + W4[1] + W4[2]) + rho13 = W4[2] / (W4[0] + W4[1] + W4[2]) + rho21 = W4[3] / (W4[3] + W4[4] + W4[5]) + rho22 = W4[4] / (W4[3] + W4[4] + W4[5]) + rho23 = W4[5] / (W4[3] + W4[4] + W4[5]) + rho31 = W4[6] / (W4[6] + W4[7] + W4[8]) + rho32 = W4[7] / (W4[6] + W4[7] + W4[8]) + rho33 = W4[8] / (W4[6] + W4[7] + W4[8]) + + Yf[0] = rho11 * Y4[0] + rho12 * Y4[1] + rho13 * Y4[2] + Yf[1] = rho21 * Y4[3] + rho22 * Y4[4] + rho23 * Y4[5] + Yf[2] = rho31 * Y4[6] + rho32 * Y4[7] + rho33 * Y4[8] + + yp11 = Y4[0] - Yf[0] + yp12 = Y4[1] - Yf[0] + yp13 = Y4[2] - Yf[0] + yp21 = Y4[3] - Yf[1] + yp22 = Y4[4] - Yf[1] + yp23 = Y4[5] - Yf[1] + yp31 = Y4[6] - Yf[2] + yp32 = Y4[7] - Yf[2] + yp33 = Y4[8] - Yf[2] + scale1 = np.sqrt(c200) + scale2 = np.sqrt(c020) + Rho1 = np.diag(rho) + Rho2 = np.array( + [[rho11, rho12, rho13], [rho21, rho22, rho23], [rho31, rho32, rho33]] + ) + Yp2 = np.array([[yp11, yp12, yp13], [yp21, yp22, yp23], [yp31, yp32, yp33]]) + Yp2s = Yp2 / scale2 + RAB = Rho1 * Rho2 + XAB = np.array( + [[xp[0], xp[1], xp[2]], [xp[0], xp[1], xp[2]], [xp[0], xp[1], xp[2]]] + ) + XABs = XAB / scale1 + YAB = Yp2 + np.diag(Yf) * np.ones(3) + YABs = YAB / scale2 + C01 = np.multiply(RAB, YABs) + Yc0 = np.ones(3) + Yc1 = XABs + Yc2 = Yp2s + A1 = np.sum(np.multiply(C01, Yc1)) + A2 = np.sum(np.multiply(C01, Yc2)) + + c101s = c101 / scale1 + c011s = c011 / scale2 + if c101s ** 2 >= c002 * (1 - small): + c101s = sign(c101s) * np.sqrt(c002) + elif c011s ** 2 >= c002 * (1 - small): + c110s = c110 / scale1 / scale2 + c011s = sign(c011s) * np.sqrt(c002) + c101s = c110s * c011s b0 = 0 b1 = c101s b2 = 0 if A2 < wsmall: - b2 = ( c011s - A1*b1 )/A2 + b2 = (c011s - A1 * b1) / A2 - Zf = b0*Yc0 + b1*Yc1 + b2*Yc2 - SUM002 = np.sum(np.multiply(RAB,Zf**2)) - mu2 = c002 - SUM002 - mu2 = max(0,mu2) - q = 0 ; eta = 1 + Zf = b0 * Yc0 + b1 * Yc1 + b2 * Yc2 + SUM002 = np.sum(np.multiply(RAB, Zf ** 2)) + mu2 = c002 - SUM002 + mu2 = max(0, mu2) + q = 0 + eta = 1 if mu2 > csmall: - SUM1 = mu2**(3/2) - SUM3 = np.sum(np.multiply(RAB,Zf**3)) - q = ( c003 - SUM3 )/SUM1 - SUM2 = mu2**2 - SUM4 = np.sum(np.multiply(RAB,Zf**4)) + 6*SUM002*mu2 - eta = ( c004 - SUM4 )/SUM2 - if eta < (q**2 + 1): + SUM1 = mu2 ** (3 / 2) + SUM3 = np.sum(np.multiply(RAB, Zf ** 3)) + q = (c003 - SUM3) / SUM1 + SUM2 = mu2 ** 2 + SUM4 = np.sum(np.multiply(RAB, Zf ** 4)) + 6 * SUM002 * mu2 + eta = (c004 - SUM4) / SUM2 + if eta < (q ** 2 + 1): if abs(q) > verysmall: - slope = (eta - 3)/q - det = 8 + slope**2 - qp = 0.5*( slope + math.sqrt(det) ) - qm = 0.5*( slope - math.sqrt(det) ) + slope = (eta - 3) / q + det = 8 + slope ** 2 + qp = 0.5 * (slope + np.sqrt(det)) + qm = 0.5 * (slope - np.sqrt(det)) if sign(q) == 1: - q = qp + q = qp else: - q = qm + q = qm else: - q = 0 - eta = q**2 + 1 - mu3 = q*mu2**(3/2) - mu4 = eta*mu2**2 + q = 0 + eta = q ** 2 + 1 + mu3 = q * mu2 ** (3 / 2) + mu4 = eta * mu2 ** 2 M5 = np.array([1, 0, mu2, mu3, mu4]) - xp11, rh11 = hyqmom3(M5,max_skewness,checks) - - rho111 = rh11[0] - rho112 = rh11[1] - rho113 = rh11[2] - - zp111 = xp11[0] - zp112 = xp11[1] - zp113 = xp11[2] - - rh12 = rh11 - xp12 = xp11 - rho121 = rh12[0] - rho122 = rh12[1] - rho123 = rh12[2] - - zp121 = xp12[0] - zp122 = xp12[1] - zp123 = xp12[2] - - rh13 = rh11 - xp13 = xp11 - rho131 = rh13[0] - rho132 = rh13[1] - rho133 = rh13[2] - - zp131 = xp13[0] - zp132 = xp13[1] - zp133 = xp13[2] - - rh21 = rh11 - xp21 = xp11 - zp211 = xp21[0] - zp212 = xp21[1] - zp213 = xp21[2] - - rho211 = rh21[0] - rho212 = rh21[1] - rho213 = rh21[2] - - rh22 = rh11 - xp22 = xp11 - zp221 = xp22[0] - zp222 = xp22[1] - zp223 = xp22[2] - - rho221 = rh22[0] - rho222 = rh22[1] - rho223 = rh22[2] - - rh23 = rh11 - xp23 = xp11 - zp231 = xp23[0] - zp232 = xp23[1] - zp233 = xp23[2] - - rho231 = rh23[0] - rho232 = rh23[1] - rho233 = rh23[2] - - rh31 = rh11 - xp31 = xp11 - rho311 = rh31[0] - rho312 = rh31[1] - rho313 = rh31[2] - - zp311 = xp31[0] - zp312 = xp31[1] - zp313 = xp31[2] - - rh32 = rh11 - xp32 = xp11 - rho321 = rh32[0] - rho322 = rh32[1] - rho323 = rh32[2] - - zp321 = xp32[0] - zp322 = xp32[1] - zp323 = xp32[2] - - rh33 = rh11 - xp33 = xp11 - rho331 = rh33[0] - rho332 = rh33[1] - rho333 = rh33[2] - - zp331 = xp33[0] - zp332 = xp33[1] - zp333 = xp33[2] - - W = np.zeros( 27 ) - W[0] = rho[0]*rho11*rho111 - W[1] = rho[0]*rho11*rho112 - W[2] = rho[0]*rho11*rho113 - W[3] = rho[0]*rho12*rho121 - W[4] = rho[0]*rho12*rho122 - W[5] = rho[0]*rho12*rho123 - W[6] = rho[0]*rho13*rho131 - W[7] = rho[0]*rho13*rho132 - W[8] = rho[0]*rho13*rho133 - W[9] = rho[1]*rho21*rho211 - W[10]= rho[1]*rho21*rho212 - W[11]= rho[1]*rho21*rho213 - W[12]= rho[1]*rho22*rho221 - W[13]= rho[1]*rho22*rho222 - W[14]= rho[1]*rho22*rho223 - W[15]= rho[1]*rho23*rho231 - W[16]= rho[1]*rho23*rho232 - W[17]= rho[1]*rho23*rho233 - W[18]= rho[2]*rho31*rho311 - W[19]= rho[2]*rho31*rho312 - W[20]= rho[2]*rho31*rho313 - W[21]= rho[2]*rho32*rho321 - W[22]= rho[2]*rho32*rho322 - W[23]= rho[2]*rho32*rho323 - W[24]= rho[2]*rho33*rho331 - W[25]= rho[2]*rho33*rho332 - W[26]= rho[2]*rho33*rho333 - W = m000*W - - abscissas = np.zeros( (27,3) ) - abscissas[0,0] = xp[0] - abscissas[1,0] = xp[0] - abscissas[2,0] = xp[0] - abscissas[3,0] = xp[0] - abscissas[4,0] = xp[0] - abscissas[5,0] = xp[0] - abscissas[6,0] = xp[0] - abscissas[7,0] = xp[0] - abscissas[8,0] = xp[0] - abscissas[9,0] = xp[1] - abscissas[10,0]= xp[1] - abscissas[11,0]= xp[1] - abscissas[12,0]= xp[1] - abscissas[13,0]= xp[1] - abscissas[14,0]= xp[1] - abscissas[15,0]= xp[1] - abscissas[16,0]= xp[1] - abscissas[17,0]= xp[1] - abscissas[18,0]= xp[2] - abscissas[19,0]= xp[2] - abscissas[20,0]= xp[2] - abscissas[21,0]= xp[2] - abscissas[22,0]= xp[2] - abscissas[23,0]= xp[2] - abscissas[24,0]= xp[2] - abscissas[25,0]= xp[2] - abscissas[26,0]= xp[2] - abscissas[:,0] += bx - - abscissas[0,1] = Yf[0]+yp11 - abscissas[1,1] = Yf[0]+yp11 - abscissas[2,1] = Yf[0]+yp11 - abscissas[3,1] = Yf[0]+yp12 - abscissas[4,1] = Yf[0]+yp12 - abscissas[5,1] = Yf[0]+yp12 - abscissas[6,1] = Yf[0]+yp13 - abscissas[7,1] = Yf[0]+yp13 - abscissas[8,1] = Yf[0]+yp13 - abscissas[9,1] = Yf[1]+yp21 - abscissas[10,1]= Yf[1]+yp21 - abscissas[11,1]= Yf[1]+yp21 - abscissas[12,1]= Yf[1]+yp22 - abscissas[13,1]= Yf[1]+yp22 - abscissas[14,1]= Yf[1]+yp22 - abscissas[15,1]= Yf[1]+yp23 - abscissas[16,1]= Yf[1]+yp23 - abscissas[17,1]= Yf[1]+yp23 - abscissas[18,1]= Yf[2]+yp31 - abscissas[19,1]= Yf[2]+yp31 - abscissas[20,1]= Yf[2]+yp31 - abscissas[21,1]= Yf[2]+yp32 - abscissas[22,1]= Yf[2]+yp32 - abscissas[23,1]= Yf[2]+yp32 - abscissas[24,1]= Yf[2]+yp33 - abscissas[25,1]= Yf[2]+yp33 - abscissas[26,1]= Yf[2]+yp33 - abscissas[:,1] += by - - abscissas[0,2] = Zf[0,0]+zp111 - abscissas[1,2] = Zf[0,0]+zp112 - abscissas[2,2] = Zf[0,0]+zp113 - abscissas[3,2] = Zf[0,1]+zp121 - abscissas[4,2] = Zf[0,1]+zp122 - abscissas[5,2] = Zf[0,1]+zp123 - abscissas[6,2] = Zf[0,2]+zp131 - abscissas[7,2] = Zf[0,2]+zp132 - abscissas[8,2] = Zf[0,2]+zp133 - abscissas[9,2] = Zf[1,0]+zp211 - abscissas[10,2]= Zf[1,0]+zp212 - abscissas[11,2]= Zf[1,0]+zp213 - abscissas[12,2]= Zf[1,1]+zp221 - abscissas[13,2]= Zf[1,1]+zp222 - abscissas[14,2]= Zf[1,1]+zp223 - abscissas[15,2]= Zf[1,2]+zp231 - abscissas[16,2]= Zf[1,2]+zp232 - abscissas[17,2]= Zf[1,2]+zp233 - abscissas[18,2]= Zf[2,0]+zp311 - abscissas[19,2]= Zf[2,0]+zp312 - abscissas[20,2]= Zf[2,0]+zp313 - abscissas[21,2]= Zf[2,1]+zp321 - abscissas[22,2]= Zf[2,1]+zp322 - abscissas[23,2]= Zf[2,1]+zp323 - abscissas[24,2]= Zf[2,2]+zp331 - abscissas[25,2]= Zf[2,2]+zp332 - abscissas[26,2]= Zf[2,2]+zp333 - abscissas[:,2] += bz + xp11, rh11 = hyqmom3(M5, max_skewness, checks) + + rho111 = rh11[0] + rho112 = rh11[1] + rho113 = rh11[2] + + zp111 = xp11[0] + zp112 = xp11[1] + zp113 = xp11[2] + + rh12 = rh11 + xp12 = xp11 + rho121 = rh12[0] + rho122 = rh12[1] + rho123 = rh12[2] + + zp121 = xp12[0] + zp122 = xp12[1] + zp123 = xp12[2] + + rh13 = rh11 + xp13 = xp11 + rho131 = rh13[0] + rho132 = rh13[1] + rho133 = rh13[2] + + zp131 = xp13[0] + zp132 = xp13[1] + zp133 = xp13[2] + + rh21 = rh11 + xp21 = xp11 + zp211 = xp21[0] + zp212 = xp21[1] + zp213 = xp21[2] + + rho211 = rh21[0] + rho212 = rh21[1] + rho213 = rh21[2] + + rh22 = rh11 + xp22 = xp11 + zp221 = xp22[0] + zp222 = xp22[1] + zp223 = xp22[2] + + rho221 = rh22[0] + rho222 = rh22[1] + rho223 = rh22[2] + + rh23 = rh11 + xp23 = xp11 + zp231 = xp23[0] + zp232 = xp23[1] + zp233 = xp23[2] + + rho231 = rh23[0] + rho232 = rh23[1] + rho233 = rh23[2] + + rh31 = rh11 + xp31 = xp11 + rho311 = rh31[0] + rho312 = rh31[1] + rho313 = rh31[2] + + zp311 = xp31[0] + zp312 = xp31[1] + zp313 = xp31[2] + + rh32 = rh11 + xp32 = xp11 + rho321 = rh32[0] + rho322 = rh32[1] + rho323 = rh32[2] + + zp321 = xp32[0] + zp322 = xp32[1] + zp323 = xp32[2] + + rh33 = rh11 + xp33 = xp11 + rho331 = rh33[0] + rho332 = rh33[1] + rho333 = rh33[2] + + zp331 = xp33[0] + zp332 = xp33[1] + zp333 = xp33[2] + + W[0] = rho[0] * rho11 * rho111 + W[1] = rho[0] * rho11 * rho112 + W[2] = rho[0] * rho11 * rho113 + W[3] = rho[0] * rho12 * rho121 + W[4] = rho[0] * rho12 * rho122 + W[5] = rho[0] * rho12 * rho123 + W[6] = rho[0] * rho13 * rho131 + W[7] = rho[0] * rho13 * rho132 + W[8] = rho[0] * rho13 * rho133 + W[9] = rho[1] * rho21 * rho211 + W[10] = rho[1] * rho21 * rho212 + W[11] = rho[1] * rho21 * rho213 + W[12] = rho[1] * rho22 * rho221 + W[13] = rho[1] * rho22 * rho222 + W[14] = rho[1] * rho22 * rho223 + W[15] = rho[1] * rho23 * rho231 + W[16] = rho[1] * rho23 * rho232 + W[17] = rho[1] * rho23 * rho233 + W[18] = rho[2] * rho31 * rho311 + W[19] = rho[2] * rho31 * rho312 + W[20] = rho[2] * rho31 * rho313 + W[21] = rho[2] * rho32 * rho321 + W[22] = rho[2] * rho32 * rho322 + W[23] = rho[2] * rho32 * rho323 + W[24] = rho[2] * rho33 * rho331 + W[25] = rho[2] * rho33 * rho332 + W[26] = rho[2] * rho33 * rho333 + W = m000 * W + + abscissas[0, 0] = xp[0] + abscissas[1, 0] = xp[0] + abscissas[2, 0] = xp[0] + abscissas[3, 0] = xp[0] + abscissas[4, 0] = xp[0] + abscissas[5, 0] = xp[0] + abscissas[6, 0] = xp[0] + abscissas[7, 0] = xp[0] + abscissas[8, 0] = xp[0] + abscissas[9, 0] = xp[1] + abscissas[10, 0] = xp[1] + abscissas[11, 0] = xp[1] + abscissas[12, 0] = xp[1] + abscissas[13, 0] = xp[1] + abscissas[14, 0] = xp[1] + abscissas[15, 0] = xp[1] + abscissas[16, 0] = xp[1] + abscissas[17, 0] = xp[1] + abscissas[18, 0] = xp[2] + abscissas[19, 0] = xp[2] + abscissas[20, 0] = xp[2] + abscissas[21, 0] = xp[2] + abscissas[22, 0] = xp[2] + abscissas[23, 0] = xp[2] + abscissas[24, 0] = xp[2] + abscissas[25, 0] = xp[2] + abscissas[26, 0] = xp[2] + abscissas[:, 0] += bx + + abscissas[0, 1] = Yf[0] + yp11 + abscissas[1, 1] = Yf[0] + yp11 + abscissas[2, 1] = Yf[0] + yp11 + abscissas[3, 1] = Yf[0] + yp12 + abscissas[4, 1] = Yf[0] + yp12 + abscissas[5, 1] = Yf[0] + yp12 + abscissas[6, 1] = Yf[0] + yp13 + abscissas[7, 1] = Yf[0] + yp13 + abscissas[8, 1] = Yf[0] + yp13 + abscissas[9, 1] = Yf[1] + yp21 + abscissas[10, 1] = Yf[1] + yp21 + abscissas[11, 1] = Yf[1] + yp21 + abscissas[12, 1] = Yf[1] + yp22 + abscissas[13, 1] = Yf[1] + yp22 + abscissas[14, 1] = Yf[1] + yp22 + abscissas[15, 1] = Yf[1] + yp23 + abscissas[16, 1] = Yf[1] + yp23 + abscissas[17, 1] = Yf[1] + yp23 + abscissas[18, 1] = Yf[2] + yp31 + abscissas[19, 1] = Yf[2] + yp31 + abscissas[20, 1] = Yf[2] + yp31 + abscissas[21, 1] = Yf[2] + yp32 + abscissas[22, 1] = Yf[2] + yp32 + abscissas[23, 1] = Yf[2] + yp32 + abscissas[24, 1] = Yf[2] + yp33 + abscissas[25, 1] = Yf[2] + yp33 + abscissas[26, 1] = Yf[2] + yp33 + abscissas[:, 1] += by + + abscissas[0, 2] = Zf[0, 0] + zp111 + abscissas[1, 2] = Zf[0, 0] + zp112 + abscissas[2, 2] = Zf[0, 0] + zp113 + abscissas[3, 2] = Zf[0, 1] + zp121 + abscissas[4, 2] = Zf[0, 1] + zp122 + abscissas[5, 2] = Zf[0, 1] + zp123 + abscissas[6, 2] = Zf[0, 2] + zp131 + abscissas[7, 2] = Zf[0, 2] + zp132 + abscissas[8, 2] = Zf[0, 2] + zp133 + abscissas[9, 2] = Zf[1, 0] + zp211 + abscissas[10, 2] = Zf[1, 0] + zp212 + abscissas[11, 2] = Zf[1, 0] + zp213 + abscissas[12, 2] = Zf[1, 1] + zp221 + abscissas[13, 2] = Zf[1, 1] + zp222 + abscissas[14, 2] = Zf[1, 1] + zp223 + abscissas[15, 2] = Zf[1, 2] + zp231 + abscissas[16, 2] = Zf[1, 2] + zp232 + abscissas[17, 2] = Zf[1, 2] + zp233 + abscissas[18, 2] = Zf[2, 0] + zp311 + abscissas[19, 2] = Zf[2, 0] + zp312 + abscissas[20, 2] = Zf[2, 0] + zp313 + abscissas[21, 2] = Zf[2, 1] + zp321 + abscissas[22, 2] = Zf[2, 1] + zp322 + abscissas[23, 2] = Zf[2, 1] + zp323 + abscissas[24, 2] = Zf[2, 2] + zp331 + abscissas[25, 2] = Zf[2, 2] + zp332 + abscissas[26, 2] = Zf[2, 2] + zp333 + abscissas[:, 2] += bz return abscissas, W diff --git a/src/qbmm_manager.py b/src/qbmm_manager.py index 9e43632..738e84c 100644 --- a/src/qbmm_manager.py +++ b/src/qbmm_manager.py @@ -9,30 +9,33 @@ """ import sys -sys.path.append('../utils/') - +sys.path.append("../utils/") from inversion import * from pretty_print_util import * import numpy as np import sympy as smp from sympy.parsing.sympy_parser import parse_expr -try: - import numba - from nquad import * -except: - print('Did not find numba! Install it for significant speedups.') - from quad import * +import numba +from nquad import * + +# try: +# import numba +# from nquad import * +# except: +# print("qbmm_mgr: Did not find numba! Install it for significant speedups.") +# from quad import * + class qbmm_manager: """ - This class manages the computation of moment-transport RHS. - It is meant to be called from within :class:`time_advancer`, with which it interfaces through :func:`compute_rhs`. - The ``config`` dictionary carries values for the following variables: + This class manages the computation of moment-transport RHS. + It is meant to be called from within :class:`time_advancer`, with which it interfaces through :func:`compute_rhs`. + The ``config`` dictionary carries values for the following variables: - :ivar governing dynamics: Governing internal dynamics - :ivar num_internal_coords: Number of internal coordinates - :ivar num_quadrature_nodes: Number of quadrature nodes + :ivar internal dynamics: Optional governing internal dynamics + :ivar num_internal_coords: Number of internal coordinates + :ivar num_nodes: Number of quadrature nodes :ivar method: Inversion method (``qmom``, ``hyqmom``, ``chyqmom``) :ivar adaptive: Adaptivity flag for ``method = qmom`` (Wheeler) :ivar max_skewness: Maximum skewness for ``method = hyqmom or chyqmom`` (hyperbolic or conditional hyperbolic) @@ -46,65 +49,59 @@ def __init__(self, config): :type config: dict """ - qbmm_config = config['qbmm'] - self.governing_dynamics = qbmm_config['governing_dynamics'] - self.num_internal_coords = qbmm_config['num_internal_coords'] - self.num_quadrature_nodes = qbmm_config['num_quadrature_nodes'] + qbmm_config = config["qbmm"] + self.num_coords = qbmm_config["num_coords"] + self.num_nodes = qbmm_config["num_nodes"] # self.poly = config['qbmm']['polydisperse'] # if self.poly: # self.num_poly_nodes = config['qbmm']['num_poly_nodes'] # self.poly_symbol = config['qbmm']['poly_symbol'] - if 'flow' in qbmm_config: - self.flow = qbmm_config['flow'] + if "internal_dynamics" in qbmm_config: + self.internal_dynamics = qbmm_config["internal_dynamics"] else: - self.flow = False + self.internal_dynamics = "" - if 'adaptive' in qbmm_config: - self.adaptive = qbmm_config['adaptive'] + if "adaptive" in qbmm_config: + self.adaptive = qbmm_config["adaptive"] else: self.adaptive = False - if 'method' in qbmm_config: - self.method = qbmm_config['method'] - else: - if self.num_internal_coords == 1: - self.method = 'hyqmom' - elif (self.num_internal_coords == 2 - or - self.num_internal_coords == 3): - self.method = 'chyqmom' - - iret = self.set_inversion( config ) + if "method" in qbmm_config: + self.method = qbmm_config["method"] + else: + if self.num_coords == 1: + self.method = "hyqmom" + elif self.num_coords == 2 or self.num_coords == 3: + self.method = "chyqmom" + + iret = self.set_inversion(config) if iret == 1: - print('qbmm_mgr: init: Configuration failed') + print("qbmm_mgr: init: Configuration failed") return # Report config - print('qbmm_mgr: init: Configuration options ready:') - print('\t flow = %s' % self.flow ) - print('\t governing_dynamics = %s' % self.governing_dynamics) - print('\t num_internal_coords = %i' % self.num_internal_coords) - print('\t method = %s' % self.method) + print("qbmm_mgr: init: Configuration options ready:") + print("\t internal_dynamics = %s" % self.internal_dynamics) + print("\t num_coords = %i" % self.num_coords) + print("\t method = %s" % self.method) # Report method-specific config - if self.method == 'qmom': - print( '\t adaptive = %s' % str( self.adaptive ) ) - if self.method == 'hyqmom' or self.method == 'chyqmom': - print( '\t max_skewness = %i' % self.max_skewness ) + if self.method == "qmom": + print("\t adaptive = %s" % str(self.adaptive)) + if self.method == "hyqmom" or self.method == "chyqmom": + print("\t max_skewness = %i" % self.max_skewness) # Determine moment indices self.moment_indices() - print( '\t num_moments = %i' % self.num_moments ) + print("\t num_moments = %i" % self.num_moments) # Determine coefficients & exponents from governing dynamics - if self.num_internal_coords < 3: - self.transport_terms() - - + if self.num_coords < 3 and self.internal_dynamics: + self.transport_terms() # RHS buffer - self.rhs = np.zeros( self.num_moments ) + self.rhs = np.zeros(self.num_moments) return @@ -114,214 +111,274 @@ def set_inversion(self, config): :param config: Configuration :type config: dict - """ - qbmm_config = config['qbmm'] + """ + qbmm_config = config["qbmm"] self.checks = True - if 'checks' in qbmm_config: - self.checks = qbmm_config['checks'] - - if self.num_internal_coords == 1: + if "checks" in qbmm_config: + self.checks = qbmm_config["checks"] + + if self.num_coords == 1: # self.moment_invert = self.moment_invert_1D # - if self.method == 'qmom': + if self.method == "qmom": # self.inversion_algorithm = wheeler - self.adaptive = False - if 'adaptive' in qbmm_config: - self.adaptive = qbmm_config['adaptive'] + self.adaptive = False + if "adaptive" in qbmm_config: + self.adaptive = qbmm_config["adaptive"] self.inversion_option = self.adaptive # - elif self.method == 'hyqmom': + elif self.method == "hyqmom": # self.inversion_algorithm = hyperbolic - self.max_skewness = 30 - if 'max_skewness' in qbmm_config: - self.max_skewness = qbmm_config['max_skewness'] + self.max_skewness = 30 + if "max_skewness" in qbmm_config: + self.max_skewness = qbmm_config["max_skewness"] self.inversion_option = self.max_skewness # else: - message = 'qbmm_mgr: set_inversion: Error: No method %s for num_internal_coords = 1' - print( message % self.method ) - return(1) + message = "qbmm_mgr: set_inversion: Error: No method %s for num_coords = 1" + print(message % self.method) + return 1 # - elif self.num_internal_coords == 2: + elif self.num_coords == 2: # self.moment_invert = self.moment_invert_2PD # - if self.method == 'chyqmom': + if self.method == "chyqmom": # self.moment_invert = conditional_hyperbolic self.inversion_algorithm = conditional_hyperbolic self.max_skewness = 30 - self.permutation = 12 - if 'max_skewness' in qbmm_config: - self.max_skewness = qbmm_config['max_skewness'] - if 'permutation' in qbmm_config: - self.permutation = qbmm_config['permutation'] + self.permutation = 12 + if "max_skewness" in qbmm_config: + self.max_skewness = qbmm_config["max_skewness"] + if "permutation" in qbmm_config: + self.permutation = qbmm_config["permutation"] self.inversion_option = self.max_skewness self.inversion_option = self.permutation self.inversion_option = self.checks # - elif self.num_internal_coords == 3: + elif self.num_coords == 3: # self.moment_invert = self.moment_invert_2PD # - if self.method == 'chyqmom': + if self.method == "chyqmom": # - self.moment_invert = conditional_hyperbolic - self.inversion_algorithm = conditional_hyperbolic + # self.moment_invert = conditional_hyperbolic + # self.inversion_algorithm = conditional_hyperbolic + self.moment_invert = chyqmom27 + self.inversion_algorithm = chyqmom27 self.max_skewness = 30 - if 'max_skewness' in qbmm_config: - self.max_skewness = qbmm_config['max_skewness'] + if "max_skewness" in qbmm_config: + self.max_skewness = qbmm_config["max_skewness"] self.inversion_option = self.max_skewness self.inversion_option = self.checks # else: - message = 'qbmm_mgr: set_inversion: Error: dimensionality %i unsupported' - print( message % self.num_internal_coords ) + message = "qbmm_mgr: set_inversion: Error: dimensionality %i unsupported" + print(message % self.num_coords) return(1) return(0) - + def moment_indices(self): """ This function sets moment indices according to dimensionality (num_coords and num_nodes) and method. """ - + ### self.num_moments = 0 # - if self.num_internal_coords == 1: + if self.num_coords == 1: # - if self.method == 'qmom': - self.indices = np.arange( 2 * self.num_quadrature_nodes ) - elif self.method == 'hyqmom': - self.indices = np.arange( 2 * ( self.num_quadrature_nodes - 1 ) + 1 ) + if self.method == "qmom": + self.indices = np.arange(2 * self.num_nodes) + elif self.method == "hyqmom": + # temp = np.arange(2 * (self.num_nodes - 1) + 1) + if self.num_nodes == 2: + self.indices = np.array([[0],[1],[2]]) + elif self.num_nodes == 3: + self.indices = np.array([[0],[1],[2],[3],[4]]) + else: + raise NotImplementedError + # - self.num_moments = len( self.indices ) + self.num_moments = len(self.indices) # - message = 'qbmm_mgr: moment_indices: ' - f_array_pretty_print( message, 'indices', self.indices ) - elif self.num_internal_coords == 2: + message = "qbmm_mgr: moment_indices: " + # f_array_pretty_print(message, "indices", self.indices) + print('self.indices',self.indices) + elif self.num_coords == 2: # - if self.method == 'chyqmom': - if self.num_quadrature_nodes == 4: - self.indices = np.array( [ [0,0], [1,0], [0,1], [2,0], [1,1], [0,2] ] ) - elif self.num_quadrature_nodes == 9: - self.indices = np.array( [ [0,0], [1,0], [0,1], [2,0], [1,1], [0,2], [3,0], [0,3], [4,0], [0,4] ] ) - else : - print( 'qbmm_mgr: moment_indices: Error: incorrect number of quadrature nodes (not 4 or 9), aborting... %i' % self.num_quadrature_nodes ) + if self.method == "chyqmom": + if self.num_nodes == 4: + self.indices = np.array( + [[0, 0], [1, 0], [0, 1], [2, 0], [1, 1], [0, 2]] + ) + elif self.num_nodes == 9: + self.indices = np.array( + [ + [0, 0], + [1, 0], + [0, 1], + [2, 0], + [1, 1], + [0, 2], + [3, 0], + [0, 3], + [4, 0], + [0, 4], + ] + ) + else: + print( + "qbmm_mgr: moment_indices: Error: incorrect number of quadrature nodes (not 4 or 9), aborting... %i" + % self.num_nodes + ) quit() else: - print( 'qbmm_mgr: moment_indices: Error: method is not chyqmom for 2 internal coordinates, aborting... %i' % self.method ) - quit() + print( + "qbmm_mgr: moment_indices: Error: method is not chyqmom for 2 internal coordinates, aborting... %i" + % self.method + ) + quit() # self.num_moments = self.indices.shape[0] # - elif self.num_internal_coords == 3: + elif self.num_coords == 3: # - if self.method == 'chyqmom': - if self.num_quadrature_nodes == 27: - self.indices = np.array( [ [0,0,0], [1,0,0], [0,1,0], [0,0,1], [2,0,0], [1,1,0], [1,0,1], [0,2,0], [0,1,1], [0,0,2], [3,0,0], [0,3,0], [0,0,3], [4,0,0], [0,4,0], [0,0,4] ] ) - else : - print( 'qbmm_mgr: moment_indices: Error: incorrect number of quadrature nodes (not 27), aborting... %i' % self.num_quadrature_nodes ) + if self.method == "chyqmom": + if self.num_nodes == 27: + self.indices = np.array( + [ + [0, 0, 0], + [1, 0, 0], + [0, 1, 0], + [0, 0, 1], + [2, 0, 0], + [1, 1, 0], + [1, 0, 1], + [0, 2, 0], + [0, 1, 1], + [0, 0, 2], + [3, 0, 0], + [0, 3, 0], + [0, 0, 3], + [4, 0, 0], + [0, 4, 0], + [0, 0, 4], + ] + ) + else: + print( + "qbmm_mgr: moment_indices: Error: incorrect number of quadrature nodes (not 27), aborting... %i" + % self.num_nodes + ) quit() - else : - print( 'qbmm_mgr: moment_indices: Error: Unsupported method, aborting...' ) + else: + print( + "qbmm_mgr: moment_indices: Error: Unsupported method, aborting..." + ) quit() self.num_moments = self.indices.shape[0] else: # - print('qbmm_mgr: moment_indices: Error: dimensionality %i unsupported' % self.num_internal_coords ) + print( + "qbmm_mgr: moment_indices: Error: dimensionality %i unsupported" + % self.num_coords + ) quit() # # # Todo: append indices for polydisperse direction r0 - # if self.num_internal_coords == 2 and self.poly: + # if self.num_coords == 2 and self.poly: # orig_idx = self.indices # self.indices = np.zeros( num_poly_nodes * len(orig_idx) ) # for j in range( num_poly_nodes ): # for i in range( len(orig_idx) ): # self.indices[i] = np.append( orig_idx[i], j ) - return + return def transport_terms(self): """ - This function determines the RHS in the moments equation for a given governing dynamics + This function determines the RHS in the moments equation for given internal dynamics """ - if self.num_internal_coords == 1: - x = smp.symbols( 'x' ) - l = smp.symbols( 'l', real = True ) - xdot = parse_expr( self.governing_dynamics ) - integrand = xdot * ( x ** ( l - 1 ) ) + if self.num_coords == 1: + x = smp.symbols("x") + l = smp.symbols("l", real=True) + xdot = parse_expr(self.internal_dynamics) + integrand = xdot * (x ** (l - 1)) self.symbolic_indices = l - elif self.num_internal_coords == 2: + elif self.num_coords == 2: # if self.poly: # r0 = smp.symbols( self.poly_symbol ) - x,xdot = smp.symbols( 'x xdot' ) - l,m = smp.symbols( 'l m', real = True ) - xddot = parse_expr( self.governing_dynamics ) - integrand = xddot * ( x ** l ) * ( xdot ** ( m - 1 ) ) - self.symbolic_indices = [l,m] + x, xdot = smp.symbols("x xdot") + l, m = smp.symbols("l m", real=True) + xddot = parse_expr(self.internal_dynamics) + integrand = xddot * (x ** l) * (xdot ** (m - 1)) + self.symbolic_indices = [l, m] - terms = smp.powsimp( smp.expand( integrand ) ).args - num_terms = len( terms ) + terms = smp.powsimp(smp.expand(integrand)).args + num_terms = len(terms) # Add constant term for 2+D problems total_num_terms = num_terms - if self.num_internal_coords == 2: + if self.num_coords == 2: total_num_terms += 1 # Initialize exponents and coefficients (weird, but works) - self.exponents = [[smp.symbols('a') for i in range(total_num_terms)] - for j in range(self.num_internal_coords)] - self.coefficients = [smp.symbols('a') for i in range(total_num_terms)] + self.exponents = [ + [smp.symbols("a") for i in range(total_num_terms)] + for j in range(self.num_coords) + ] + self.coefficients = [smp.symbols("a") for i in range(total_num_terms)] # Everything is simpler if now transferred into numpy arrays - self.exponents = np.array(self.exponents).T + self.exponents = np.array(self.exponents).T self.coefficients = np.array(self.coefficients).T # Loop over terms - for i in range( num_terms ): - self.exponents[i,0] = terms[i].as_coeff_exponent(x)[1] - if self.num_internal_coords == 1: - self.coefficients[i] = l * smp.poly( terms[i]).coeffs()[0] + for i in range(num_terms): + self.exponents[i, 0] = terms[i].as_coeff_exponent(x)[1] + if self.num_coords == 1: + self.coefficients[i] = l * smp.poly(terms[i]).coeffs()[0] else: - self.exponents[i,1] = terms[i].as_coeff_exponent(xdot)[1] - self.coefficients[i] = m * smp.poly( terms[i] ).coeffs()[0] + self.exponents[i, 1] = terms[i].as_coeff_exponent(xdot)[1] + self.coefficients[i] = m * smp.poly(terms[i]).coeffs()[0] # Add extra constant term if in 2D - if self.num_internal_coords == 2: - self.exponents[ num_terms, 0 ] = l - 1 - self.exponents[ num_terms, 1 ] = m + 1 - self.coefficients[ num_terms ] = l - - self.num_coefficients = len( self.coefficients ) - self.num_exponents = len( self.exponents ) - - message = 'qbmm_mgr: transport_terms: ' - for i in range( total_num_terms ): - sym_array_pretty_print( message, 'exponents', self.exponents[i,:] ) - - message = 'qbmm_mgr: transport_terms: ' - sym_array_pretty_print( message, 'coefficients', self.coefficients ) - - for i in range( self.num_coefficients ): - if self.num_internal_coords == 1: - self.coefficients[i] = smp.lambdify([l],self.coefficients[i]) - for j in range(self.num_internal_coords): - self.exponents[i,j] = smp.lambdify([l],self.exponents[i,j]) - elif self.num_internal_coords == 2: - self.coefficients[i] = smp.lambdify([l,m],self.coefficients[i]) - for j in range(self.num_internal_coords): - self.exponents[i,j] = smp.lambdify([l,m],self.exponents[i,j]) + if self.num_coords == 2: + self.exponents[num_terms, 0] = l - 1 + self.exponents[num_terms, 1] = m + 1 + self.coefficients[num_terms] = l + + self.num_coefficients = len(self.coefficients) + self.num_exponents = len(self.exponents) + + # message = 'qbmm_mgr: transport_terms: ' + # for i in range( total_num_terms ): + # sym_array_pretty_print( message, 'exponents', self.exponents[i,:] ) + + # message = 'qbmm_mgr: transport_terms: ' + # sym_array_pretty_print( message, 'coefficients', self.coefficients ) + + for i in range(self.num_coefficients): + if self.num_coords == 1: + self.coefficients[i] = smp.lambdify([l], self.coefficients[i]) + for j in range(self.num_coords): + self.exponents[i, j] = smp.lambdify([l], self.exponents[i, j]) + elif self.num_coords == 2: + self.coefficients[i] = smp.lambdify([l, m], self.coefficients[i]) + for j in range(self.num_coords): + self.exponents[i, j] = smp.lambdify([l, m], self.exponents[i, j]) return - + def moment_invert_1D(self, moments): """ This function inverts tracked moments into a quadrature rule in 1D @@ -335,9 +392,9 @@ def moment_invert_1D(self, moments): >>> xi, wts = qbmm_mgr.moment_inver(moments) - and qbmm_manager automatically selects moment_invert_1D based if ``num_internal_coords = 1`` + and qbmm_manager automatically selects moment_invert_1D based if ``num_coords = 1`` """ - return self.inversion_algorithm( moments, self.inversion_option ) + return self.inversion_algorithm(moments, self.inversion_option) def moment_invert_2PD(self, moments): """ @@ -352,33 +409,51 @@ def moment_invert_2PD(self, moments): >>> xi, wts = qbmm_mgr.moment_inver(moments) - and qbmm_manager automatically selects moment_invert_2PD based if ``num_internal_coords > 1`` + and qbmm_manager automatically selects moment_invert_2PD based if ``num_coords > 1`` """ - return self.inversion_algorithm( moments, self.indices, self.inversion_option ) - - def projection(self, weights, abscissas, indices): - """ - This function reconstructs moments (indices) from quadrature weights and abscissas - - :param weights: Quadrature weights - :param abscissas: Quadrature abscissas - :param indices: Full moment set indices - :type weights: array like - :type abscissas: array like - :type weights: array like - :return: projected moments - :rtype: array like - """ - abscissas = np.array(abscissas) - moments = np.zeros( len(indices) ) - for i_index in range( len(indices) ): - if self.num_internal_coords == 3: - moments[i_index] = quadrature_3d(weights, abscissas, indices[i_index], self.num_quadrature_nodes) - if self.num_internal_coords == 2: - moments[i_index] = quadrature_2d(weights, abscissas, indices[i_index], self.num_quadrature_nodes) - elif self.num_internal_coords == 1: - moments[i_index] = quadrature_1d(weights, abscissas, indices[i_index]) - return moments + return self.inversion_algorithm(moments, self.indices, self.inversion_option) + + # def projection(self, weights, abscrissas): + # """ + # This function reconstructs moments (indices) from quadrature weights and abscissas + + # :param weights: Quadrature weights + # :param abscissas: Quadrature abscissas + # :type weights: array like + # :type abscissas: array like + # :type weights: array like + # :return: projected moments + # :rtype: array like + # """ + # return self.projection(weights, abscissas, self.indices) + + # def projection(self, weights, abscissas, indices): + # """ + # This function reconstructs moments (indices) from quadrature weights and abscissas + + # :param weights: Quadrature weights + # :param abscissas: Quadrature abscissas + # :param indices: Full moment set indices + # :type weights: array like + # :type abscissas: array like + # :type weights: array like + # :return: projected moments + # :rtype: array like + # """ + # abscissas = np.array(abscissas) + # moments = np.zeros(len(indices)) + # for i in range(len(indices)): + # if self.num_coords == 3: + # moments[i] = quadrature_3d( + # weights, abscissas, indices[i], self.num_nodes + # ) + # if self.num_coords == 2: + # moments[i] = quadrature_2d( + # weights, abscissas, indices[i], self.num_nodes + # ) + # elif self.num_coords == 1: + # moments[i] = quadrature_1d(weights, abscissas, indices[i]) + # return moments def compute_rhs(self, moments, rhs): """ @@ -389,42 +464,76 @@ def compute_rhs(self, moments, rhs): """ # Compute abscissas and weights from moments - if self.num_internal_coords == 1: - abscissas, weights = self.moment_invert( moments ) + if self.num_coords == 1: + abscissas, weights = self.moment_invert(moments) else: - abscissas, weights = self.moment_invert( moments,self.indices ) + abscissas_temp, weights = self.moment_invert(moments, self.indices) + # New: Transpose + if self.num_coords == 2: + abscissas = [ + abscissas_temp[:,0], + abscissas_temp[:,1]] + elif self.num_coords == 3: + abscissas = [ + abscissas_temp[:,0], + abscissas_temp[:,1], + abscissas_temp[:,2]] # Loop over moments - for i_moment in range( self.num_moments ): + for i_moment in range(self.num_moments): # Evalue RHS terms - if self.num_internal_coords == 1: - exponents = [np.double(self.exponents[j,0](self.indices[i_moment])) - for j in range(self.num_exponents)] - coefficients = [np.double(self.coefficients[j](self.indices[i_moment]) ) - for j in range(self.num_coefficients)] - elif self.num_internal_coords == 2: - exponents = [ [ \ - np.double(self.exponents[j,0](self.indices[i_moment][0],self.indices[i_moment][1])), \ - np.double(self.exponents[j,1](self.indices[i_moment][0],self.indices[i_moment][1])) \ - ] for j in range(self.num_exponents)] - coefficients = [ \ - np.double(self.coefficients[j](self.indices[i_moment][0],self.indices[i_moment][1])) - for j in range(self.num_coefficients)] - else : - print('num_internal_coords',self.num_internal_coords,'not supported yet') + if self.num_coords == 1: + exponents = [ + np.double(self.exponents[j, 0](self.indices[i_moment])) + for j in range(self.num_exponents) + ] + coefficients = [ + np.double(self.coefficients[j](self.indices[i_moment])) + for j in range(self.num_coefficients) + ] + elif self.num_coords == 2: + exponents = [ + [ + np.double( + self.exponents[j, 0]( + self.indices[i_moment][0], self.indices[i_moment][1] + ) + ), + np.double( + self.exponents[j, 1]( + self.indices[i_moment][0], self.indices[i_moment][1] + ) + ), + ] + for j in range(self.num_exponents) + ] + coefficients = [ + np.double( + self.coefficients[j]( + self.indices[i_moment][0], self.indices[i_moment][1] + ) + ) + for j in range(self.num_coefficients) + ] + else: + print( + "num_coords", self.num_coords, "not supported yet" + ) quit() # Put them in numpy arrays - np_exponents = np.array( exponents ) - np_coefficients = np.array( coefficients ) + np_exponents = np.array(exponents) + np_coefficients = np.array(coefficients) # Project back to moments - rhs_moments = self.projection( weights, abscissas, np_exponents ) + rhs_moments = projection(weights, abscissas, np_exponents, + self.num_coords, self.num_nodes) # Compute RHS - rhs[i_moment] = np.dot( np_coefficients, rhs_moments ) + rhs[i_moment] = np.dot(np_coefficients, rhs_moments) # - projected_moments = self.projection( weights, abscissas, self.indices ) - for i_moment in range( self.num_moments ): - moments[i_moment] = projected_moments[i_moment] + projected_moments = projection(weights, abscissas, self.indices, + self.num_coords, self.num_nodes) + + for i_moment in range(self.num_moments): + moments[i_moment] = projected_moments[i_moment] # return - diff --git a/src/simulation_domain.py b/src/simulation_domain.py new file mode 100644 index 0000000..918c640 --- /dev/null +++ b/src/simulation_domain.py @@ -0,0 +1,181 @@ +import sys +sys.path.append('../utils/') + +from qbmm_manager import * +from stats_util import * +from jets_util import * +from itertools import product + +import numba +from nquad import * + +# try: +# import numba +# from nquad import * +# except: +# print("simulation_domain: Did not find numba! Install it for significant speedups.") +# from quad import * + +class simulation_domain(): + """ + This class handles grid and state for flow problems + """ + def __init__(self, config): + """ + Constructor + + :param config: Configuration + :type config: dict + """ + # Create a qbmm manager + self.qbmm_mgr = qbmm_manager(config) + + # Get config + self.num_dim = config["domain"]["num_dim"] + self.num_points = config["domain"]["num_points"] + self.extents = config["domain"]["grid_extents"] + + self.grid_spacing = (self.extents[1] - self.extents[0]) / (self.num_points - 2) + + if "flow" in config["domain"]: + self.flow = config["domain"]["flow"] + self.flux = np.zeros([self.num_points, self.qbmm_mgr.num_moments]) + else: + self.flow = False + + print("domain: init: Configuration options ready:") + print("\t num_dim = %i" % self.num_dim) + print("\t num_points = %i" % self.num_points) + print("\t extents = [%.4E, %.4E]" % (self.extents[0], self.extents[1])) + print("\t grid spacing = %.4E" % self.grid_spacing) + + # Initialize grid state & rhs + self.state = np.zeros([self.num_points, self.qbmm_mgr.num_moments]) + self.rhs = np.zeros([self.num_points, self.qbmm_mgr.num_moments]) + + # Initialize weights and abscissas + self.weights = np.zeros([self.num_points, self.qbmm_mgr.num_nodes]) + self.abscissas = np.zeros([self.num_points, self.qbmm_mgr.num_coords, self.qbmm_mgr.num_nodes]) + + return + + + def initialize_state_uniform(self, mu, sigma): + """ + Initialize grid state + """ + + raw_moments = raw_gaussian_moments_univar(self.qbmm_mgr.num_moments, mu, sigma) + for i_point in range(self.num_points): + self.state[i_point] = raw_moments + + print(self.state) + return + + + def initialize_state_jets(self, state): + """ + Initialize grid state to 1D crossing-jets + + :param state: domain moments + :type state: array like + """ + # print('shape num coords :',self.qbmm_mgr.num_coords) + wts_left, wts_right, xi_left, xi_right = jet_initialize_moments( + self.qbmm_mgr.num_coords, + self.qbmm_mgr.num_nodes) + + disc_loc = 0.125 + n_pt = len(self.weights) - 2 + disc_idx = int(n_pt * disc_loc) - 2 + print('Dislocation index is ', disc_idx, ' out of ', n_pt, ' points') + print("abscissas left: ", xi_left[0,:]) + print("abscissas right: ", xi_right[0,:]) + + # Populate weights + self.weights[:disc_idx] = wts_left + self.weights[-disc_idx:] = wts_right + # Populate abscissas + self.abscissas[:disc_idx] = xi_left + self.abscissas[-disc_idx:] = xi_right + + # Populate state + moments_left = projection(wts_left, xi_left, self.qbmm_mgr.indices, + self.qbmm_mgr.num_coords, self.qbmm_mgr.num_nodes) + moments_right = projection(wts_right, xi_right, self.qbmm_mgr.indices, + self.qbmm_mgr.num_coords, self.qbmm_mgr.num_nodes) + + state[:disc_idx] = moments_left + state[-disc_idx:] = moments_right + state[0] = moments_right + state[-1] = moments_left + + # print("Domain: state: ", state[0,:]) + # print("Domain: weights: ", self.weights[0,:]) + # print("Domain: abscissas: ", self.abscissas[:,0,:]) + # raise Exception + + return + + + def max_abscissa(self): + """ + Return the maximum value of the x-abscissa in the simulation domain + """ + return max(0.01,abs(self.abscissas[:,0,:]).max()) + + def compute_rhs(self, state): + """ + Compute moment transport fluxes, source terms. + + :param state: domain moments + :type state: array like + """ + if self.flow: + self.grid_inversion(state) + domain_get_fluxes(self.weights, self.abscissas, self.qbmm_mgr.indices, + self.num_points, self.qbmm_mgr.num_moments, + self.qbmm_mgr.num_nodes, self.flux) + self.rhs = self.flux / self.grid_spacing + if np.isnan(np.sum(self.rhs)): + raise Exception('NaN in RHS') + elif self.qbmm_mgr.internal_dynamics: + internal_rhs = np.zeros([self.num_points, self.qbmm_mgr.num_moments]) + self.qbmm_mgr.compute_rhs(state, internal_rhs) + self.rhs = internal_rhs + else: + return self.zeros(self.rhs.shape) + + return self.rhs + + + def grid_inversion(self,state): + """ + Invert moments to weights/abscissas over whole grid + + :param state: domain moments + :type state: array like + """ + if self.qbmm_mgr.num_coords == 3: + domain_invert_3d(state, self.qbmm_mgr.indices, + self.weights, self.abscissas, + self.num_points, self.qbmm_mgr.num_coords, + self.qbmm_mgr.num_nodes) + elif self.qbmm_mgr.num_coords == 2: + domain_invert_2d(state, self.qbmm_mgr.indices, + self.weights, self.abscissas, + self.num_points, self.qbmm_mgr.num_coords, + self.qbmm_mgr.num_nodes) + elif self.qbmm_mgr.num_coords == 1: + domain_invert_1d(state, self.qbmm_mgr.indices, + self.weights, self.abscissas, + self.num_points, self.qbmm_mgr.num_coords, + self.qbmm_mgr.num_nodes) + else: + raise NotImplementedError + + def project(self, state): + domain_project(state, self.qbmm_mgr.indices, + self.weights, self.abscissas, + self.num_points, self.qbmm_mgr.num_coords, + self.qbmm_mgr.num_nodes) diff --git a/src/timing_function_calls.ipynb b/src/timing_function_calls.ipynb new file mode 100644 index 0000000..dfce742 --- /dev/null +++ b/src/timing_function_calls.ipynb @@ -0,0 +1,69 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 18, + "id": "moderate-shock", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "f 0.07301426999998739\n", + "g 0.12769997799998123\n" + ] + } + ], + "source": [ + "import timeit\n", + "def f(): \n", + " return\n", + "print('f',timeit.timeit(f))\n", + "\n", + "def g(): \n", + " x = 5 + 4\n", + " y = x^4+2+6+x\n", + " return\n", + "print('g',timeit.timeit(g))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "educational-funeral", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "later-accident", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.2" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/test/test_cqmom.py b/test/test_cqmom.py deleted file mode 100644 index 7d33fba..0000000 --- a/test/test_cqmom.py +++ /dev/null @@ -1,33 +0,0 @@ -import sys -sys.path.append('../src/') -from qbmm_manager import * - -if __name__ == '__main__': - - np.set_printoptions( formatter = { 'float': '{: 0.4E}'.format } ) - - ### - ### Say hello - print('test_wheeler: Testing CQMOM for moment inversion') - - ### - ### QBMM Configuration - print('test_wheeler: Configuring and initializing qbmm') - - config = {} - config['governing_dynamics'] = ' dx + x = 1' - config['num_internal_coords'] = 2 - config['num_quadrature_nodes'] = 4 - config['method'] = 'cqmom' - config['adaptive'] = False - config['max_skewness'] = 30 - - ### - ### QBMM - qbmm_mgr = qbmm_manager( config ) - - indices = np.ones( [ 4, 2 ] ) - moments = np.ones( [ 4, 2 ] ) - qbmm_mgr.moment_invert( moments, indices ) - - exit diff --git a/test/test_project.py b/test/test_project.py deleted file mode 100644 index 24732be..0000000 --- a/test/test_project.py +++ /dev/null @@ -1,75 +0,0 @@ -import sys -sys.path.append('../src/') -from qbmm_manager import * -sys.path.append('../utils/') -from stats_util import * - -def test_project_1D(qbmm_mgr,mu,sig,tol): - init_moments = raw_gaussian_moments_univar( qbmm_mgr.num_moments, mu, sig ) - - abscissas, weights = qbmm_mgr.moment_invert( init_moments ) - projected_moments = qbmm_mgr.projection( weights, abscissas, qbmm_mgr.indices ) - - err = np.linalg.norm(projected_moments - init_moments) - success = ( err < tol ) - - return success - -def test_project_2D(qbmm_mgr,mu,sig,tol): - init_moments = raw_gaussian_moments_bivar( - qbmm_mgr.indices, - mu[0], mu[1], - sig[0], sig[1] - ) - - abscissas, weights = qbmm_mgr.moment_invert( init_moments, qbmm_mgr.indices ) - projected_moments = qbmm_mgr.projection( weights, abscissas, qbmm_mgr.indices ) - - err = np.linalg.norm(projected_moments - init_moments) - success = ( err < tol ) - - return success - -if __name__ == '__main__': - - np.set_printoptions( formatter = { 'float': '{: 0.4E}'.format } ) - - print('Unit test: test_projection') - - tol = 1.e-13 # Why? - success = True - - # 1D - print('Unit test: test 1d projection') - config = {} - config['qbmm'] = {} - config['qbmm']['governing_dynamics'] = '4*x - 2*x**2' - config['qbmm']['num_internal_coords'] = 1 - config['qbmm']['num_quadrature_nodes'] = 3 - config['qbmm']['method'] = 'qmom' - qbmm_mgr = qbmm_manager( config ) - - mu = 0.0 - sigma = 1.1 - success *= test_project_1D(qbmm_mgr,mu,sigma,tol) - - # 2D - print('Unit test: test 2d projection') - config = {} - config['qbmm'] = {} - config['qbmm']['governing_dynamics'] = ' - xdot - x ' - config['qbmm']['num_internal_coords'] = 2 - config['qbmm']['num_quadrature_nodes'] = 4 - config['qbmm']['method'] = 'chyqmom' - qbmm_mgr = qbmm_manager( config ) - - mu = [1.1,0.1] - sigma = [1.0,1.0] - success *= test_project_2D(qbmm_mgr,mu,sigma,tol) - - if success == True: - print('test_projection: ' + '\033[92m' + 'test passed' + '\033[0m') - else: - print('test_projection: ' + '\033[91m' + 'test failed ' + '\033[0m') - - exit() diff --git a/test/test_qbmm.py b/test/test_qbmm.py index b02765d..e671b1e 100644 --- a/test/test_qbmm.py +++ b/test/test_qbmm.py @@ -1,9 +1,54 @@ import sys -sys.path.append('../src/') -sys.path.append('../utils/') + +sys.path.append("../src/") +sys.path.append("../utils/") from qbmm_manager import * from stats_util import * import numpy.polynomial.hermite as hermite_poly +import pytest + + +def test_project(): + tol = 1.0e-10 + success = True + + # 1D + config = {} + config["qbmm"] = {} + config["qbmm"]["governing_dynamics"] = "4*x - 2*x**2" + config["qbmm"]["num_internal_coords"] = 1 + config["qbmm"]["num_quadrature_nodes"] = 3 + config["qbmm"]["method"] = "qmom" + qbmm_mgr = qbmm_manager(config) + + mu = 0.0 + sig = 1.1 + init_moments = raw_gaussian_moments_univar(qbmm_mgr.num_moments, mu, sig) + abscissas, weights = qbmm_mgr.moment_invert(init_moments) + projected_moments = qbmm_mgr.projection(weights, abscissas, qbmm_mgr.indices) + err_1D = np.linalg.norm(projected_moments - init_moments) + + # 2D + config = {} + config["qbmm"] = {} + config["qbmm"]["governing_dynamics"] = " - xdot - x " + config["qbmm"]["num_internal_coords"] = 2 + config["qbmm"]["num_quadrature_nodes"] = 4 + config["qbmm"]["method"] = "chyqmom" + qbmm_mgr = qbmm_manager(config) + + mu = [1.1, 0.1] + sig = [1.0, 1.0] + init_moments = raw_gaussian_moments_bivar( + qbmm_mgr.indices, mu[0], mu[1], sig[0], sig[1] + ) + abscissas, weights = qbmm_mgr.moment_invert(init_moments, qbmm_mgr.indices) + projected_moments = qbmm_mgr.projection(weights, abscissas, qbmm_mgr.indices) + err_2D = np.linalg.norm(projected_moments - init_moments) + + assert err_1D < tol + assert err_2D < tol + def test_wheeler(): """ @@ -12,49 +57,50 @@ def test_wheeler(): """ num_nodes = 4 config = {} - config['qbmm'] = {} - config['qbmm']['governing_dynamics'] = '4*x - 2*x**2' - config['qbmm']['num_internal_coords'] = 1 - config['qbmm']['num_quadrature_nodes'] = num_nodes - config['qbmm']['method'] = 'qmom' + config["qbmm"] = {} + config["qbmm"]["governing_dynamics"] = "4*x - 2*x**2" + config["qbmm"]["num_internal_coords"] = 1 + config["qbmm"]["num_quadrature_nodes"] = num_nodes + config["qbmm"]["method"] = "qmom" ### ### QBMM - qbmm_mgr = qbmm_manager( config ) - + qbmm_mgr = qbmm_manager(config) + ### ### Tests - + # Anticipate success tol = 1.0e-10 # Round-off error in moments computation success = True - + # Test 1 - mu = 5.0 + mu = 5.0 sigma = 1.0 ### ### Reference solution - sqrt_pi = np.sqrt( np.pi ) - sqrt_two = np.sqrt( 2.0 ) - h_abs, h_wts = hermite_poly.hermgauss( num_nodes ) + sqrt_pi = np.sqrt(np.pi) + sqrt_two = np.sqrt(2.0) + h_abs, h_wts = hermite_poly.hermgauss(num_nodes) g_abs = sqrt_two * sigma * h_abs + mu g_wts = h_wts / sqrt_pi - + ### ### QBMM - moments = raw_gaussian_moments_univar( qbmm_mgr.num_moments, mu, sigma ) - my_abs, my_wts = qbmm_mgr.moment_invert( moments ) + moments = raw_gaussian_moments_univar(qbmm_mgr.num_moments, mu, sigma) + my_abs, my_wts = qbmm_mgr.moment_invert(moments) ### ### Errors & Report diff_abs = my_abs - g_abs diff_wts = my_wts - g_wts - error_abs = np.linalg.norm( my_abs - g_abs ) - error_wts = np.linalg.norm( my_wts - g_wts ) + error_abs = np.linalg.norm(my_abs - g_abs) + error_wts = np.linalg.norm(my_wts - g_wts) assert error_abs < tol assert error_wts < tol + if __name__ == "__main__": if len(sys.argv) > 1: diff --git a/test/test_wheeler.py b/test/test_wheeler.py deleted file mode 100644 index fb985b7..0000000 --- a/test/test_wheeler.py +++ /dev/null @@ -1,120 +0,0 @@ -import sys -sys.path.append('../src/') -from qbmm_manager import * -import numpy.polynomial.hermite as hermite_poly - -def gauss_hermite(num_nodes): - """ - This function returns Gauss-Hermite abscissas & weights - from numpy's Hermite-polynomial engine - """ - h_abs, h_wts = hermite_poly.hermgauss( num_nodes ) - return h_abs, h_wts - -def generate_gaussian_moments(mu, sigma): - """ - This function returns the first 8 Gaussian moments, - per Marchisio & Fox, pg. 52 - """ - moments = np.zeros( 8 ) - moments[0] = 1.0 - moments[1] = mu - moments[2] = mu*mu + sigma*sigma - moments[3] = mu**3.0 + 3.0 * mu * sigma**2.0 - moments[4] = mu**4.0 + 6.0 * (mu**2.0) * (sigma**2.0) + 3.0 * sigma**4.0 - moments[5] = mu**5.0 + 10.0 * (mu**3.0) * (sigma**2.0) + 15.0 * mu * sigma**4.0 - moments[6] = mu**6.0 + 15.0 * (mu**4.0) * (sigma**2.0) + 45.0 * (mu**2.0) * (sigma**4.0) + 15.0 * sigma**6.0 - moments[7] = mu**7.0 + 21.0 * (mu**5.0) * (sigma**2.0) + 105.0 * (mu**3.0) * (sigma**4.0) + 105.0 * mu * sigma**6.0 - - return moments - -def test_wheeler(test, mu, sigma, qbmm_mgr, tol): - """ - This function tests QBMM Wheeler inversion by comparing - against numpy's Gauss-Hermite for given mu and sigma - """ - ### - ### Reference solution - num_nodes = 4 - sqrt_pi = np.sqrt( np.pi ) - sqrt_two = np.sqrt( 2.0 ) - h_abs, h_wts = gauss_hermite( num_nodes ) - g_abs = sqrt_two * sigma * h_abs + mu - g_wts = h_wts / sqrt_pi - - ### - ### QBMM - moments = generate_gaussian_moments( mu, sigma ) - my_abs, my_wts = qbmm_mgr.moment_invert( moments ) - - ### - ### Errors & Report - - error_abs = np.linalg.norm( my_abs - g_abs ) - error_wts = np.linalg.norm( my_wts - g_wts ) - - success = ( error_abs < tol and error_wts < tol ) - - if not success: - print('test_wheeler: Test(%i): Moments = [{:s}]'.format( ', '.join( [ '{:.4e}'.format(p) for p in moments ] ) ) % test ) - print('test_wheeler: Test(%i): Reference solution:' % test) - print('\t abscissas = [{:s}]'.format( ', '.join( [ '{:.4e}'.format(p) for p in g_abs ] ) ) ) - print('\t weights = [{:s}]'.format( ', '.join( [ '{:.4e}'.format(p) for p in g_wts ] ) ) ) - - print('test_wheeler: Test(%i): QBMM solution:' % test) - print('\t abscissas = [{:s}]'.format( ', '.join( [ '{:.4e}'.format(p) for p in my_abs ] ) ) ) - print('\t weights = [{:s}]'.format( ', '.join( [ '{:.4e}'.format(p) for p in my_wts ] ) ) ) - - print('test_wheeler: Test(%i): Errors:' % test) - print('\t abscissas: error = %.4E' % error_abs ) - print('\t weights: error = %.4E' % error_wts ) - - return success - -if __name__ == '__main__': - - np.set_printoptions( formatter = { 'float': '{: 0.4E}'.format } ) - - ### - ### Say hello - print('test_wheeler: Testing Wheeler algorithm for moment inversion') - - ### - ### QBMM Configuration - print('test_wheeler: Configuring and initializing qbmm') - - nnodes = [ 1, 2, 3, 4, 5 ] - for n in nnodes: - config = {} - config['qbmm'] = {} - config['qbmm']['governing_dynamics'] = '4*x - 2*x**2' - config['qbmm']['num_internal_coords'] = 1 - config['qbmm']['num_quadrature_nodes'] = n - config['qbmm']['method'] = 'qmom' - - ### - ### QBMM - qbmm_mgr = qbmm_manager( config ) - - ### - ### Tests - - # Anticipate success - tol = 1.0e-13 # Why? - success = True - - # Test 1 - mu = 0.0 - sigma = 1.0 - success *= test_wheeler( 1, mu, sigma, qbmm_mgr, tol ) - - ### Test 2 - mu = 5.0 - success *= test_wheeler( 2, mu, sigma, qbmm_mgr, tol ) - - if success == True: - print('test_wheeler: ' + '\033[92m' + 'test passed' + '\033[0m') - else: - print('test_wheeler: ' + '\033[91m' + 'test failed ' + '\033[0m') - - exit() diff --git a/utils/euler_util.py b/utils/euler_util.py deleted file mode 100644 index c053ac0..0000000 --- a/utils/euler_util.py +++ /dev/null @@ -1,46 +0,0 @@ -import numpy as np - -def local_flux(weight, abscissa, moment): - - flux = 1.0 - flux *= abscissa[0] ** moment[0] - flux *= abscissa[1] ** moment[1] - flux *= abscissa[2] ** moment[2] - # for i_coord in range( 3, num_coords ): flux *= abscissa[i_cord] ** moment[i_coord] - flux *= weight - return flux - -def moment_fluxes(indices, wts_left, wts_right, xi_left, xi_right): - """ - Computes moment fluxes - - inputs: - ------- - num_nodes: number of quadrature nodes, depends on inversion algorithm - indices: moment indices, size [ num_moments, num_internal_coords ] - wts_left: weights on the left side, size [ num_nodes ] - wts_right: weights on the right side, size [ num_nodes ] - xi_left: abscissas on the left side, size [ num_internal_coords, num_nodes ] - xi_right: abscissas on the right side, size [ num_internal_corods, num_nodes ] - """ - - num_moments = len( indices ) - num_coords, num_nodes = xi_left.shape - - flux = np.zeros( num_moments ) - - for i_moment in range( num_moments ): - for i_node in range( num_nodes ): - - # compute local fluxes - flux_left = local_flux( wts_left[i_node], xi_left[:,i_node], indices[i_moment,:] ) - flux_right = local_flux( wts_right[i_node],xi_right[:,i_node],indices[i_moment,:] ) - - # limiter (?) - flux_left = flux_left * max( xi_left[0,i_node], 0.0 ) - flux_right = flux_right * min( xi_right[0,i_node], 0.0 ) - - # quadrature - flux[i_moment] += flux_left + flux_right - - return flux diff --git a/utils/jets_util.py b/utils/jets_util.py index c2229f5..d536606 100644 --- a/utils/jets_util.py +++ b/utils/jets_util.py @@ -1,96 +1,188 @@ import sys -sys.path.append('../src/') + +sys.path.append("../src/") from qbmm_manager import * import numpy as np -def jet_initialize_moments(qbmm_mgr): - - num_nodes = qbmm_mgr.num_quadrature_nodes - - wts_left = np.zeros( num_nodes ) - wts_right = np.zeros( num_nodes ) - xi_left = np.zeros( [3,num_nodes] ) - xi_right = np.zeros( [3,num_nodes] ) - - num_left = 1.0 - num_right = 0.5 - - u_left = 0.01 - u_right = 1.0 - - wts_left[0] = 0.125 * num_left - wts_left[2] = 0.125 * num_left - wts_left[6] = 0.125 * num_left - wts_left[8] = 0.125 * num_left - wts_left[18] = 0.125 * num_left - wts_left[20] = 0.125 * num_left - wts_left[24] = 0.125 * num_left - wts_left[26] = 0.125 * num_left - - wts_right[0] = 0.125 * num_right - wts_right[2] = 0.125 * num_right - wts_right[6] = 0.125 * num_right - wts_right[8] = 0.125 * num_right - wts_right[18] = 0.125 * num_right - wts_right[20] = 0.125 * num_right - wts_right[24] = 0.125 * num_right - wts_right[26] = 0.125 * num_right - - xi_left[0,0] = u_left - xi_left[0,2] = u_left - xi_left[0,6] = u_left - xi_left[0,8] = u_left - xi_left[0,18] = u_left - xi_left[0,20] = u_left - xi_left[0,24] = u_left - xi_left[0,26] = u_left - - xi_right[0,0] = u_right - xi_right[0,2] = u_right - xi_right[0,6] = u_right - xi_right[0,8] = u_right - xi_right[0,18] = u_right - xi_right[0,20] = u_right - xi_right[0,24] = u_right - xi_right[0,26] = u_right - - xi_left[1,0] = u_left - xi_left[1,2] = u_left - xi_left[1,6] = -1.0 - xi_left[1,8] = -1.0 - xi_left[1,18] = 1.0 - xi_left[1,20] = 1.0 - xi_left[1,24] = -1.0 - xi_left[1,26] = -1.0 - - xi_right[1,0] = 1.0 - xi_right[1,2] = 1.0 - xi_right[1,6] = -1.0 - xi_right[1,8] = -1.0 - xi_right[1,18] = 1.0 - xi_right[1,20] = 1.0 - xi_right[1,24] = -1.0 - xi_right[1,26] = -1.0 - - xi_left[2,0] = -1.0 - xi_left[2,2] = 1.0 - xi_left[2,6] = -1.0 - xi_left[2,8] = 1.0 - xi_left[2,18] = -1.0 - xi_left[2,20] = 1.0 - xi_left[2,24] = -1.0 - xi_left[2,26] = 1.0 - xi_right[2,:] = xi_left[2,:] - - xi_left[1,:] += 0.5 * xi_left[0,:] - xi_right[1,:] += 0.5 * xi_right[0,:] - - xi_left[2,:] = 0.0 - xi_right[2,:] = 0.0 - - print(qbmm_mgr.indices) - - moments_left = qbmm_mgr.projection( wts_left, xi_left, qbmm_mgr.indices ) - moments_right = qbmm_mgr.projection( wts_right, xi_right, qbmm_mgr.indices ) - - return moments_left, moments_right + +def jet_initialize_moments(num_coords,num_nodes): + print('Initializing jet with', num_coords, 'coordinates') + if num_coords == 3: + return init_3d_jet(num_nodes) + elif num_coords == 2: + return init_2d_jet(num_nodes) + elif num_coords == 1: + return init_1d_jet(num_nodes) + +def init_3d_jet(num_nodes): + wts_left = np.zeros(num_nodes) + wts_right = np.zeros(num_nodes) + xi_left = np.zeros([3, num_nodes]) + xi_right = np.zeros([3, num_nodes]) + + w_left = 1 + w_right = 0.5 + + u_left = 0.01 + u_right = 1 + + wts_left[0] = 0.125 * w_left + wts_left[2] = 0.125 * w_left + wts_left[6] = 0.125 * w_left + wts_left[8] = 0.125 * w_left + wts_left[18] = 0.125 * w_left + wts_left[20] = 0.125 * w_left + wts_left[24] = 0.125 * w_left + wts_left[26] = 0.125 * w_left + + wts_right[0] = 0.125 * w_right + wts_right[2] = 0.125 * w_right + wts_right[6] = 0.125 * w_right + wts_right[8] = 0.125 * w_right + wts_right[18] = 0.125 * w_right + wts_right[20] = 0.125 * w_right + wts_right[24] = 0.125 * w_right + wts_right[26] = 0.125 * w_right + + xi_left[0, 0] = u_left + xi_left[0, 2] = u_left + xi_left[0, 6] = u_left + xi_left[0, 8] = u_left + xi_left[0, 18] = u_left + xi_left[0, 20] = u_left + xi_left[0, 24] = u_left + xi_left[0, 26] = u_left + + xi_right[0, 0] = u_right + xi_right[0, 2] = u_right + xi_right[0, 6] = u_right + xi_right[0, 8] = u_right + xi_right[0, 18] = u_right + xi_right[0, 20] = u_right + xi_right[0, 24] = u_right + xi_right[0, 26] = u_right + + xi_left[1, 0] = u_left + xi_left[1, 2] = u_left + xi_left[1, 6] = -1 + xi_left[1, 8] = -1 + xi_left[1, 18] = 1 + xi_left[1, 20] = 1 + xi_left[1, 24] = -1 + xi_left[1, 26] = -1 + + xi_right[1, 0] = 1 + xi_right[1, 2] = 1 + xi_right[1, 6] = -1 + xi_right[1, 8] = -1 + xi_right[1, 18] = 1 + xi_right[1, 20] = 1 + xi_right[1, 24] = -1 + xi_right[1, 26] = -1 + + xi_left[2, 0] = -1 + xi_left[2, 2] = 1 + xi_left[2, 6] = -1 + xi_left[2, 8] = 1 + xi_left[2, 18] = -1 + xi_left[2, 20] = 1 + xi_left[2, 24] = -1 + xi_left[2, 26] = 1 + xi_right[2, :] = xi_left[2, :] + + xi_left[1, :] += 0.5 * xi_left[0, :] + xi_right[1, :] += 0.5 * xi_right[0, :] + + xi_left[2, :] = 0.0 + xi_right[2, :] = 0.0 + + return wts_left, wts_right, xi_left, xi_right + + +def init_2d_jet(num_nodes): + wts_left = np.zeros(num_nodes) + wts_right = np.zeros(num_nodes) + xi_left = np.zeros([2, num_nodes]) + xi_right = np.zeros([2, num_nodes]) + + w_left = 1 + w_right = 0.5 + + u_left = 0.01 + u_right = 1 + + wts_left[0] = w_left / 4. + wts_left[2] = w_left / 4. + wts_left[6] = w_left / 4. + wts_left[8] = w_left / 4. + + wts_right[0] = w_right / 4. + wts_right[2] = w_right / 4. + wts_right[6] = w_right / 4. + wts_right[8] = w_right / 4. + + xi_left[0, 0] = u_left + xi_left[0, 2] = u_left + xi_left[0, 6] = u_left + xi_left[0, 8] = u_left + + xi_right[0, 0] = u_right + xi_right[0, 2] = u_right + xi_right[0, 6] = u_right + xi_right[0, 8] = u_right + + # Taking from y-direction in 3D jet + xi_left[1, 0] = u_left + xi_left[1, 2] = u_left + # xi_left[1, 0] = 1 + # xi_left[1, 2] = 1 + xi_left[1, 6] = -1 + xi_left[1, 8] = -1 + + # Try taking from z-direction instead + # xi_left[1, 0] = -1 + # xi_left[1, 2] = 1 + # xi_left[1, 6] = -1 + # xi_left[1, 8] = 1 + + xi_right[1, :] = xi_left[1, :] + + xi_left[1, :] += 0.5 * xi_left[0, :] + xi_right[1, :] += 0.5 * xi_right[0, :] + + # xi_left[1, :] = 0. + # xi_right[1, :] = 0. + + # xi_right[:, :] += np.random.rand()*1e-7 + # xi_left[:, :] += np.random.rand()*1e-7 + # wts_left[:] += np.random.rand()*1e-7 + # wts_right[ :] += np.random.rand()*1e-7 + + return wts_left, wts_right, xi_left, xi_right + + +def init_1d_jet(num_nodes): + wts_left = np.zeros(num_nodes) + wts_right = np.zeros(num_nodes) + xi_left = np.zeros([1, num_nodes]) + xi_right = np.zeros([1, num_nodes]) + + w_left = 1 + w_right = 0.5 + + u_left = 0.01 + u_right = 1 + + wts_left[0] = w_left/2. + wts_left[2] = w_left/2. + + wts_right[0] = w_right/2. + wts_right[2] = w_right/2. + + xi_left[0, 0] = u_left + xi_left[0, 2] = u_left + + xi_right[0, 0] = u_right + xi_right[0, 2] = u_right + + return wts_left, wts_right, xi_left, xi_right diff --git a/utils/nquad.py b/utils/nquad.py index 0170b40..2e6fe8e 100644 --- a/utils/nquad.py +++ b/utils/nquad.py @@ -1,32 +1,238 @@ -import math import numpy as np -from numba import jit - -@jit(nopython=True) -def quadrature_1d(weights, abscissas, moment_index): - """ - This function computes quadrature in 1D - Inputs: - - weights: quadrature weights - Return: - """ - xi_to_idx = abscissas ** moment_index - q = np.dot( weights, xi_to_idx ) +from inversion import hyqmom3, chyqmom9, chyqmom27 +from numba import njit, prange + +@njit +def projection( + weights, abscissas, indices, + num_coords, num_nodes): + + moments = np.zeros(indices.shape[0]) + ni = len(indices) + for i in range(ni): + if num_coords == 3: + moments[i] = quadrature_3d( + weights, abscissas, indices[i], num_nodes + ) + if num_coords == 2: + moments[i] = quadrature_2d( + weights, abscissas, indices[i], num_nodes + ) + if num_coords == 1: + moments[i] = quadrature_1d( + weights, abscissas, indices[i], num_nodes + ) + return moments + + +@njit +def quadrature_1d(weights, abscissas, moment_index, num_quadrature_nodes): + # xi_to_idx = abscissas ** moment_index + # q = np.dot(weights, xi_to_idx) + # return q + + q = 0.0 + for i in range(num_quadrature_nodes): + q += ( + weights[i] + * (abscissas[0][i] ** moment_index[0]) + ) return q -@jit(nopython=True) + +@njit def quadrature_2d(weights, abscissas, moment_index, num_quadrature_nodes): - q = 0. - for i in range( num_quadrature_nodes ): - q += weights[i]*(abscissas[0][i]**moment_index[0]) * \ - (abscissas[1][i]**moment_index[1]) + q = 0.0 + for i in range(num_quadrature_nodes): + q += ( + weights[i] + * (abscissas[0][i] ** moment_index[0]) + * (abscissas[1][i] ** moment_index[1]) + ) return q + +@njit def quadrature_3d(weights, abscissas, moment_index, num_quadrature_nodes): - q = 0. - for i in range( num_quadrature_nodes ): - q += weights[i] * \ - abscissas[0,i]**moment_index[0] * \ - abscissas[1,i]**moment_index[1] * \ - abscissas[2,i]**moment_index[2] + q = 0.0 + for i in range(num_quadrature_nodes): + q += ( + weights[i] + * abscissas[0, i] ** moment_index[0] + * abscissas[1, i] ** moment_index[1] + * abscissas[2, i] ** moment_index[2] + ) return q + + +@njit(parallel=True) +def flux_quadrature(wts, xi, indices, num_moments, num_nodes, num_points): + flux_min = np.zeros((num_points,num_moments,num_nodes)) + flux_max = np.zeros((num_points,num_moments,num_nodes)) + for i in prange(num_points): + for m in range(num_moments): + for n in range(num_nodes): + flux = ( + wts[i,n] + * xi[i, 0, n]**indices[m, 0] + ) + + if len(indices[0]) == 2: + flux *= ( + xi[i, 1, n]**indices[m, 1] + ) + + if len(indices[0]) == 3: + flux *= ( + xi[i, 1, n]**indices[m, 1] * + xi[i, 2, n]**indices[m, 2] + ) + + flux_min[i,m,n] = flux * min(xi[i, 0, n], 0.) + flux_max[i,m,n] = flux * max(xi[i, 0, n], 0.) + + return flux_min, flux_max + + +@njit +def domain_get_fluxes(weights, abscissas, indices, num_points, num_moments, num_nodes, flux): + + f_min, f_max = flux_quadrature( + weights, abscissas, + indices, num_moments, + num_nodes, num_points) + + f_sum = np.zeros_like(flux) + f_sum[1:-1,:] = np.sum(f_max[0:-2,:] + f_min[1:-1,:],axis=2) + flux[1:-2] = f_sum[1:-2] - f_sum[2:-1] + + +@njit(parallel=True) +def domain_project(state, indices, weights, abscissas, num_points, num_coords, num_nodes): + # Input: Weights, abscissas, indices, num_points, num_coords, num_nodes + # Output: State + + for i_point in prange(1, num_points-1): + state[i_point] = projection( + weights[i_point], + abscissas[i_point], + indices, + num_coords, + num_nodes) + + # Boundary conditions + state[0] = projection(weights[-2], abscissas[-2], indices, num_coords, num_nodes) + state[-1] = projection(weights[1], abscissas[1], indices, num_coords, num_nodes) + + +@njit(parallel=True) +def domain_invert_3d(state, indices, weights, abscissas, num_points, num_coords, num_nodes): + # Input: State, indices, num_points, num_coords, num_nodes + # Output: Weights, abscissas + + for i_point in prange(1, num_points-1): + # Invert + xi, wts = chyqmom27(state[i_point], indices) + abscissas[i_point] = xi.T + weights[i_point] = wts + + # Boundary conditions + xi, wts = chyqmom27(state[0], indices) + abscissas[0] = xi.T + weights[0] = wts + + xi, wts = chyqmom27(state[-1], indices) + abscissas[-1] = xi.T + weights[-1] = wts + + return weights, abscissas + +# @njit +def domain_invert_2d(state, indices, weights, abscissas, num_points, num_coords, num_nodes): + # Input: State, indices, num_points, num_coords, num_nodes + # Output: Weights, abscissas + + for i_point in range(1, num_points-1): + # Invert + xi, wts = chyqmom9(state[i_point], indices) + abscissas[i_point] = xi.T + weights[i_point] = wts + + # Boundary conditions + xi, wts = chyqmom9(state[0], indices) + abscissas[0] = xi.T + weights[0] = wts + + xi, wts = chyqmom9(state[-1], indices) + abscissas[-1] = xi.T + weights[-1] = wts + + return weights, abscissas + +# @njit +def domain_invert_1d(state, indices, weights, abscissas, num_points, num_coords, num_nodes): + # Input: State, indices, num_points, num_coords, num_nodes + # Output: Weights, abscissas + + for i_point in range(1, num_points-1): + # Invert + xi, wts = hyqmom3(state[i_point]) + abscissas[i_point] = xi.T + weights[i_point] = wts + + # Boundary conditions + xi, wts = hyqmom3(state[0]) + abscissas[0] = xi.T + weights[0] = wts + + xi, wts = hyqmom3(state[-1]) + abscissas[-1] = xi.T + weights[-1] = wts + + return weights, abscissas + +# @njit +# def domain_invert_2d(state, indices, weights, abscissas, num_points, num_coords, num_nodes): + +# for i_point in range(1, num_points-1): +# # print(' ipt', i_point) +# # print('State before', state[i_point]) +# # Invert +# xi, wts = chyqmom9(state[i_point], indices) +# abscissas[i_point] = xi.T +# weights[i_point] = wts +# # Project +# state[i_point] = projection(wts, xi.T, indices, num_coords, num_nodes) +# # print('State after', state[i_point]) +# if np.isnan(abscissas[i_point]).any() or \ +# np.isnan(weights[i_point]).any(): +# ii = np.where(np.isnan(abscissas[i_point]))[0] +# # print('state', iii, 'is nan:', state[iii]) +# print('abs:', abscissas[i_point]) +# print('w:',weights[i_point]) +# print('found nan in wts or absc') +# raise Exception() + +# # Boundary conditions +# state[0] = projection(weights[-2], abscissas[-2], indices, num_coords, num_nodes) +# state[-1] = projection(weights[1], abscissas[1], indices, num_coords, num_nodes) + +# # print('get pt0') +# # print(state[0]) +# xi, wts = chyqmom9(state[0], indices) +# abscissas[0] = xi.T +# weights[0] = wts +# # print('get pt-1') +# # print(state[-1]) +# xi, wts = chyqmom9(state[-1], indices) +# abscissas[-1] = xi.T +# weights[-1] = wts +# # print(' finish update quad 2d') + +# print('max = ', np.max(state[:,:])) +# if np.isnan(state).any(): +# print('found nan') +# ii = np.where(np.isnan(state))[0] +# iii = ii[0] +# print('state', iii, 'is nan:', state[iii]) +# raise Exception() diff --git a/utils/nquad_bak.py b/utils/nquad_bak.py new file mode 100644 index 0000000..1141cee --- /dev/null +++ b/utils/nquad_bak.py @@ -0,0 +1,158 @@ +import numpy as np +from inversion import chyqmom9, chyqmom27 +from numba import njit + +@njit +def projection( + weights, abscissas, indices, + num_coords, num_nodes): + + moments = np.zeros(indices.shape[0]) + ni = len(indices) + for i in range(ni): + moments[i] = quadrature( + weights, abscissas, indices[i], num_coords + ) + return moments + + +@njit +def quadrature(weights, abscissas, moment_index, num_coords): + if num_coords == 3: + return np.sum( + weights[:] + * abscissas[0, :] ** moment_index[0] + * abscissas[1, :] ** moment_index[1] + * abscissas[2, :] ** moment_index[2] + ) + elif num_coords == 2: + return np.sum( + weights[:] + * abscissas[0, :] ** moment_index[0] + * abscissas[1, :] ** moment_index[1] + ) + elif num_coords == 1: + return np.sum( + weights[:] + * abscissas[0, :] ** moment_index[0] + ) + +@njit +def quadrature_limited(weights, abscissas, moment_index, num_coords, num_nodes): + flux_min = np.zeros(num_nodes) + flux_max = np.zeros(num_nodes) + for n in range(num_nodes): + q = weights[n] * np.prod(abscissas[:,n]**moment_index[:]) + flux_min[n] = q * min(abscissas[0, n], 0.) + flux_max[n] = q * max(abscissas[0, n], 0.) + return np.sum(flux_min), np.sum(flux_max) + +@njit +def flux_quadrature(weights, abscissas, indices, num_moments, num_nodes, num_points): + flux_min = np.zeros((num_points,num_moments)) + flux_max = np.zeros((num_points,num_moments)) + num_coords = len(indices[0]) + for i in range(num_points): + for m in range(num_moments): + flux_min[i,m], flux_max[i,m] = quadrature_limited( + weights[i],abscissas[i],indices[m],num_coords,num_nodes) + return flux_min, flux_max + +@njit +def domain_get_fluxes(weights, abscissas, indices, num_points, num_moments, num_nodes, flux): + + f_min, f_max = flux_quadrature( + weights, abscissas, + indices, num_moments, + num_nodes, num_points) + + f_sum = np.zeros_like(flux) + f_sum[1:-1] = f_max[0:-2] + f_min[1:-1] + flux[1:-2] = f_sum[1:-2] - f_sum[2:-1] + + +@njit +def domain_project(state, indices, weights, abscissas, num_points, num_coords, num_nodes): + # Input: Weights, abscissas, indices, num_points, num_coords, num_nodes + # Output: State + + for i_point in range(1, num_points-1): + state[i_point] = projection( + weights[i_point], + abscissas[i_point], + indices, + num_coords, + num_nodes) + + # Boundary conditions + state[0] = projection(weights[-2], abscissas[-2], indices, num_coords, num_nodes) + state[-1] = projection(weights[1], abscissas[1], indices, num_coords, num_nodes) + + +@njit +def domain_invert_3d(state, indices, weights, abscissas, num_points, num_coords, num_nodes): + # Input: State, indices, num_points, num_coords, num_nodes + # Output: Weights, abscissas + + for i_point in range(1, num_points-1): + # Invert + xi, wts = chyqmom27(state[i_point], indices) + abscissas[i_point] = xi.T + weights[i_point] = wts + + # Boundary conditions + xi, wts = chyqmom27(state[0], indices) + abscissas[0] = xi.T + weights[0] = wts + + xi, wts = chyqmom27(state[-1], indices) + abscissas[-1] = xi.T + weights[-1] = wts + + return weights, abscissas + +@njit +def domain_invert_2d(state, indices, weights, abscissas, num_points, num_coords, num_nodes): + + for i_point in range(1, num_points-1): + # print(' ipt', i_point) + # print('State before', state[i_point]) + # Invert + xi, wts = chyqmom9(state[i_point], indices) + abscissas[i_point] = xi.T + weights[i_point] = wts + # Project + state[i_point] = projection(wts, xi.T, indices, num_coords, num_nodes) + # print('State after', state[i_point]) + if np.isnan(abscissas[i_point]).any() or \ + np.isnan(weights[i_point]).any(): + ii = np.where(np.isnan(abscissas[i_point]))[0] + # print('state', iii, 'is nan:', state[iii]) + print('abs:', abscissas[i_point]) + print('w:',weights[i_point]) + print('found nan in wts or absc') + raise Exception() + + # Boundary conditions + state[0] = projection(weights[-2], abscissas[-2], indices, num_coords, num_nodes) + state[-1] = projection(weights[1], abscissas[1], indices, num_coords, num_nodes) + + # print('get pt0') + # print(state[0]) + xi, wts = chyqmom9(state[0], indices) + abscissas[0] = xi.T + weights[0] = wts + # print('get pt-1') + # print(state[-1]) + xi, wts = chyqmom9(state[-1], indices) + abscissas[-1] = xi.T + weights[-1] = wts + # print(' finish update quad 2d') + + print('max = ', np.max(state[:,:])) + if np.isnan(state).any(): + print('found nan') + ii = np.where(np.isnan(state))[0] + iii = ii[0] + print('state', iii, 'is nan:', state[iii]) + raise Exception() diff --git a/utils/pretty_print_util.py b/utils/pretty_print_util.py index ce5a731..02f89b1 100644 --- a/utils/pretty_print_util.py +++ b/utils/pretty_print_util.py @@ -1,24 +1,28 @@ def i_scalar_pretty_print(message, scalar_name, scalar): - message += scalar_name + ' = %i' - print( message % scalar ) + message += scalar_name + " = %i" + print(message % scalar) + def f_scalar_pretty_print(message, scalar_name, scalar): - message += scalar_name + ' = %.4E' - print( message % scalar ) + message += scalar_name + " = %.4E" + print(message % scalar) + def i_array_pretty_print(message, array_name, array): - message += array_name + ' = [{:s}]' - print( message.format( ', '.join( [ str(a) for a in array ] ) ) ) + message += array_name + " = [{:s}]" + print(message.format(", ".join([str(a) for a in array]))) + def f_array_pretty_print(message, array_name, array): - message += array_name + ' = [{:s}]' - print( message.format( ', '.join( [ '{:.16e}'.format(a) for a in array ] ) ) ) + message += array_name + " = [{:s}]" + print(message.format(", ".join(["{:.16e}".format(a) for a in array]))) + def sym_array_pretty_print(message, array_name, array): - message += array_name + ' = [{:s}]' - print( message.format( ', '.join( [ str(a) for a in array ] ) ) ) + message += array_name + " = [{:s}]" + print(message.format(", ".join([str(a) for a in array]))) diff --git a/utils/pretty_print_util.pyc b/utils/pretty_print_util.pyc deleted file mode 100644 index 41370b5..0000000 Binary files a/utils/pretty_print_util.pyc and /dev/null differ diff --git a/utils/quad.py b/utils/quad.py index 20d5ceb..f4ab38a 100644 --- a/utils/quad.py +++ b/utils/quad.py @@ -1,23 +1,31 @@ import math import numpy as np + def quadrature_1d(weights, abscissas, moment_index): xi_to_idx = abscissas ** moment_index - q = np.dot( weights, xi_to_idx ) + q = np.dot(weights, xi_to_idx) return q + def quadrature_2d(weights, abscissas, moment_index, num_quadrature_nodes): - q = 0. - for i in range( num_quadrature_nodes ): - q += weights[i]*(abscissas[0][i]**moment_index[0]) * \ - (abscissas[1][i]**moment_index[1]) + q = 0.0 + for i in range(num_quadrature_nodes): + q += ( + weights[i] + * (abscissas[0][i] ** moment_index[0]) + * (abscissas[1][i] ** moment_index[1]) + ) return q + def quadrature_3d(weights, abscissas, moment_index, num_quadrature_nodes): - q = 0. - for i in range( num_quadrature_nodes ): - q += weights[i] * \ - abscissas[0,i]**moment_index[0] * \ - abscissas[1,i]**moment_index[1] * \ - abscissas[2,i]**moment_index[2] + q = 0.0 + for i in range(num_quadrature_nodes): + q += ( + weights[i] + * abscissas[0, i] ** moment_index[0] + * abscissas[1, i] ** moment_index[1] + * abscissas[2, i] ** moment_index[2] + ) return q diff --git a/utils/stats_util.py b/utils/stats_util.py index 3d6142b..a7b3048 100644 --- a/utils/stats_util.py +++ b/utils/stats_util.py @@ -4,96 +4,159 @@ import numpy as np from stats_util import * + def raw_gaussian_moments_univar(num_moments, mu, sigma): """ - This function returns raw 1D-Gaussian moments as a function of + This function returns raw 1D-Gaussian moments as a function of mean (mu) and standard deviation (sigma) """ - moments = np.zeros( num_moments ) - for i_moment in range( num_moments ): - moments[ i_moment ] = stats.norm.moment( i_moment, mu, sigma ) + moments = np.zeros(num_moments) + for i_moment in range(num_moments): + moments[i_moment] = stats.norm.moment(i_moment, mu, sigma) return moments - + + def raw_gaussian_moments_bivar(indices, mu1, mu2, sig1, sig2): """ - This function returns raw 2D-Gaussian moments as a function of + This function returns raw 2D-Gaussian moments as a function of means (mu_1,mu_2) and standard deviations (sigma_1,sigma_2) """ num_moments = len(indices) - moments = np.zeros( num_moments) - for i_moment in range( num_moments ): + moments = np.zeros(num_moments) + for i_moment in range(num_moments): i = indices[i_moment][0] j = indices[i_moment][1] - moments[ i_moment ] = (1./math.pi) * 2**((-4.+i+j)/2.) * \ - math.exp( -(mu1**2./(2. * sig1**2))-(mu2**2./(2 * sig2**2.)) ) * \ - sig1**(-1.+i) * sig2**(-1+j) * \ - ( \ - -math.sqrt(2.) * \ - (-1.+(-1.)**i) * \ - mu1 * \ - sc.gamma(1.+i/2.) * \ - sc.hyp1f1(1+i/2.,3./2.,mu1**2./(2. * sig1**2.)) \ - + (1.+(-1.)**i) * \ - sig1 * \ - sc.gamma((1.+i)/2.) * \ - sc.hyp1f1((1.+i)/2.,1./2.,mu1**2./(2. * sig1**2.)) - ) * \ - ( \ - -math.sqrt(2.) * \ - (-1+(-1)**j) * \ - mu2 * \ - sc.gamma(1.+j/2.) * \ - sc.hyp1f1(1.+j/2.,3./2.,mu2**2./(2. * sig2**2.)) + \ - (1.+(-1)**j) * \ - sig2 * \ - sc.gamma((1.+j)/2.) * \ - sc.hyp1f1((1.+j)/2.,1./2.,mu2**2./(2. * sig2**2.)) \ - ) + moments[i_moment] = ( + (1.0 / math.pi) + * 2 ** ((-4.0 + i + j) / 2.0) + * math.exp( + -(mu1 ** 2.0 / (2.0 * sig1 ** 2)) - (mu2 ** 2.0 / (2 * sig2 ** 2.0)) + ) + * sig1 ** (-1.0 + i) + * sig2 ** (-1 + j) + * ( + -math.sqrt(2.0) + * (-1.0 + (-1.0) ** i) + * mu1 + * sc.gamma(1.0 + i / 2.0) + * sc.hyp1f1(1 + i / 2.0, 3.0 / 2.0, mu1 ** 2.0 / (2.0 * sig1 ** 2.0)) + + (1.0 + (-1.0) ** i) + * sig1 + * sc.gamma((1.0 + i) / 2.0) + * sc.hyp1f1( + (1.0 + i) / 2.0, 1.0 / 2.0, mu1 ** 2.0 / (2.0 * sig1 ** 2.0) + ) + ) + * ( + -math.sqrt(2.0) + * (-1 + (-1) ** j) + * mu2 + * sc.gamma(1.0 + j / 2.0) + * sc.hyp1f1(1.0 + j / 2.0, 3.0 / 2.0, mu2 ** 2.0 / (2.0 * sig2 ** 2.0)) + + (1.0 + (-1) ** j) + * sig2 + * sc.gamma((1.0 + j) / 2.0) + * sc.hyp1f1( + (1.0 + j) / 2.0, 1.0 / 2.0, mu2 ** 2.0 / (2.0 * sig2 ** 2.0) + ) + ) + ) return moments -def raw_gaussian_moments_trivar(indices,mu1,mu2,mu3,sig1,sig2,sig3): + +def raw_gaussian_moments_trivar(indices, mu1, mu2, mu3, sig1, sig2, sig3): num_moments = len(indices) - moments = np.zeros( num_moments) + moments = np.zeros(num_moments) i = 0 for idx2 in indices: idx = idx2.tolist() - if idx == [0, 0, 0]: moments[i] = 1 - if idx == [0, 0, 1]: moments[i] = mu3 - if idx == [0, 0, 2]: moments[i] = mu3**2 + sig3**2 - if idx == [0, 0, 3]: moments[i] = mu3**3 + 3*mu3*sig3**2 - if idx == [0, 0, 4]: moments[i] = mu3**4 + 6*mu3**2*sig3**2 + 3*sig3**4 - if idx == [0, 1, 0]: moments[i] = mu2 - if idx == [0, 1, 1]: moments[i] = mu2*mu3 - if idx == [0, 1, 2]: moments[i] = mu2*mu3**2 + mu2*sig3**2 - if idx == [0, 1, 3]: moments[i] = mu2*mu3**3 + 3*mu2*mu3*sig3**2 - if idx == [0, 2, 0]: moments[i] = mu2**2 + sig2**2 - if idx == [0, 2, 1]: moments[i] = mu2**2*mu3 + mu3*sig2**2 - if idx == [0, 2, 2]: moments[i] = mu2**2*mu3**2 + mu3**2*sig2**2 + mu2**2*sig3**2 + sig2**2*sig3**2 - if idx == [0, 3, 0]: moments[i] = mu2**3 + 3*mu2*sig2**2 - if idx == [0, 3, 1]: moments[i] = mu2**3*mu3 + 3*mu2*mu3*sig2**2 - if idx == [0, 4, 0]: moments[i] = mu2**4 + 6*mu2**2*sig2**2 + 3*sig2**4 - if idx == [1, 0, 0]: moments[i] = mu1 - if idx == [1, 0, 1]: moments[i] = mu1*mu3 - if idx == [1, 0, 2]: moments[i] = mu1*mu3**2 + mu1*sig3**2 - if idx == [1, 0, 3]: moments[i] = mu1*mu3**3 + 3*mu1*mu3*sig3**2 - if idx == [1, 1, 0]: moments[i] = mu1*mu2 - if idx == [1, 1, 1]: moments[i] = mu1*mu2*mu3 - if idx == [1, 1, 2]: moments[i] = mu1*mu2*mu3**2 + mu1*mu2*sig3**2 - if idx == [1, 2, 0]: moments[i] = mu1*mu2**2 + mu1*sig2**2 - if idx == [1, 2, 1]: moments[i] = mu1*mu2**2*mu3 + mu1*mu3*sig2**2 - if idx == [1, 3, 0]: moments[i] = mu1*mu2**3 + 3*mu1*mu2*sig2**2 - if idx == [2, 0, 0]: moments[i] = mu1**2 + sig1**2 - if idx == [2, 0, 1]: moments[i] = mu1**2*mu3 + mu3*sig1**2 - if idx == [2, 0, 2]: moments[i] = mu1**2*mu3**2 + mu3**2*sig1**2 + mu1**2*sig3**2 + sig1**2*sig3**2 - if idx == [2, 1, 0]: moments[i] = mu1**2*mu2 + mu2*sig1**2 - if idx == [2, 1, 1]: moments[i] = mu1**2*mu2*mu3 + mu2*mu3*sig1**2 - if idx == [2, 2, 0]: moments[i] = mu1**2*mu2**2 + mu2**2*sig1**2 + mu1**2*sig2**2 + sig1**2*sig2**2 - if idx == [3, 0, 0]: moments[i] = mu1**3 + 3*mu1*sig1**2 - if idx == [3, 0, 1]: moments[i] = mu1**3*mu3 + 3*mu1*mu3*sig1**2 - if idx == [3, 1, 0]: moments[i] = mu1**3*mu2 + 3*mu1*mu2*sig1**2 - if idx == [4, 0, 0]: moments[i] = mu1**4 + 6*mu1**2*sig1**2 + 3*sig1**4 + if idx == [0, 0, 0]: + moments[i] = 1 + if idx == [0, 0, 1]: + moments[i] = mu3 + if idx == [0, 0, 2]: + moments[i] = mu3 ** 2 + sig3 ** 2 + if idx == [0, 0, 3]: + moments[i] = mu3 ** 3 + 3 * mu3 * sig3 ** 2 + if idx == [0, 0, 4]: + moments[i] = mu3 ** 4 + 6 * mu3 ** 2 * sig3 ** 2 + 3 * sig3 ** 4 + if idx == [0, 1, 0]: + moments[i] = mu2 + if idx == [0, 1, 1]: + moments[i] = mu2 * mu3 + if idx == [0, 1, 2]: + moments[i] = mu2 * mu3 ** 2 + mu2 * sig3 ** 2 + if idx == [0, 1, 3]: + moments[i] = mu2 * mu3 ** 3 + 3 * mu2 * mu3 * sig3 ** 2 + if idx == [0, 2, 0]: + moments[i] = mu2 ** 2 + sig2 ** 2 + if idx == [0, 2, 1]: + moments[i] = mu2 ** 2 * mu3 + mu3 * sig2 ** 2 + if idx == [0, 2, 2]: + moments[i] = ( + mu2 ** 2 * mu3 ** 2 + + mu3 ** 2 * sig2 ** 2 + + mu2 ** 2 * sig3 ** 2 + + sig2 ** 2 * sig3 ** 2 + ) + if idx == [0, 3, 0]: + moments[i] = mu2 ** 3 + 3 * mu2 * sig2 ** 2 + if idx == [0, 3, 1]: + moments[i] = mu2 ** 3 * mu3 + 3 * mu2 * mu3 * sig2 ** 2 + if idx == [0, 4, 0]: + moments[i] = mu2 ** 4 + 6 * mu2 ** 2 * sig2 ** 2 + 3 * sig2 ** 4 + if idx == [1, 0, 0]: + moments[i] = mu1 + if idx == [1, 0, 1]: + moments[i] = mu1 * mu3 + if idx == [1, 0, 2]: + moments[i] = mu1 * mu3 ** 2 + mu1 * sig3 ** 2 + if idx == [1, 0, 3]: + moments[i] = mu1 * mu3 ** 3 + 3 * mu1 * mu3 * sig3 ** 2 + if idx == [1, 1, 0]: + moments[i] = mu1 * mu2 + if idx == [1, 1, 1]: + moments[i] = mu1 * mu2 * mu3 + if idx == [1, 1, 2]: + moments[i] = mu1 * mu2 * mu3 ** 2 + mu1 * mu2 * sig3 ** 2 + if idx == [1, 2, 0]: + moments[i] = mu1 * mu2 ** 2 + mu1 * sig2 ** 2 + if idx == [1, 2, 1]: + moments[i] = mu1 * mu2 ** 2 * mu3 + mu1 * mu3 * sig2 ** 2 + if idx == [1, 3, 0]: + moments[i] = mu1 * mu2 ** 3 + 3 * mu1 * mu2 * sig2 ** 2 + if idx == [2, 0, 0]: + moments[i] = mu1 ** 2 + sig1 ** 2 + if idx == [2, 0, 1]: + moments[i] = mu1 ** 2 * mu3 + mu3 * sig1 ** 2 + if idx == [2, 0, 2]: + moments[i] = ( + mu1 ** 2 * mu3 ** 2 + + mu3 ** 2 * sig1 ** 2 + + mu1 ** 2 * sig3 ** 2 + + sig1 ** 2 * sig3 ** 2 + ) + if idx == [2, 1, 0]: + moments[i] = mu1 ** 2 * mu2 + mu2 * sig1 ** 2 + if idx == [2, 1, 1]: + moments[i] = mu1 ** 2 * mu2 * mu3 + mu2 * mu3 * sig1 ** 2 + if idx == [2, 2, 0]: + moments[i] = ( + mu1 ** 2 * mu2 ** 2 + + mu2 ** 2 * sig1 ** 2 + + mu1 ** 2 * sig2 ** 2 + + sig1 ** 2 * sig2 ** 2 + ) + if idx == [3, 0, 0]: + moments[i] = mu1 ** 3 + 3 * mu1 * sig1 ** 2 + if idx == [3, 0, 1]: + moments[i] = mu1 ** 3 * mu3 + 3 * mu1 * mu3 * sig1 ** 2 + if idx == [3, 1, 0]: + moments[i] = mu1 ** 3 * mu2 + 3 * mu1 * mu2 * sig1 ** 2 + if idx == [4, 0, 0]: + moments[i] = mu1 ** 4 + 6 * mu1 ** 2 * sig1 ** 2 + 3 * sig1 ** 4 i = i + 1 return moments