diff --git a/.gitmodules b/.gitmodules index baa7d310..76f0f9ac 100644 --- a/.gitmodules +++ b/.gitmodules @@ -7,3 +7,6 @@ [submodule "external/catch2"] path = external/catch2 url = https://github.com/catchorg/Catch2.git +[submodule "external/HighFive"] + path = external/HighFive + url = https://github.com/BlueBrain/HighFive.git diff --git a/CMakeLists.txt b/CMakeLists.txt index e73422a0..8dfc61c5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -17,6 +17,9 @@ endif() project(athenaPK LANGUAGES CXX) +# Adding HighFive library +add_subdirectory(external/HighFive) + set(CMAKE_CXX_EXTENSIONS OFF) set(CMAKE_EXPORT_COMPILE_COMMANDS ON) set(CMAKE_CXX_STANDARD 17) @@ -83,5 +86,4 @@ add_subdirectory(src) if (AthenaPK_ENABLE_TESTING) include(CTest) add_subdirectory(tst/regression) -endif() - +endif() \ No newline at end of file diff --git a/external/HighFive b/external/HighFive new file mode 160000 index 00000000..7340d1fc --- /dev/null +++ b/external/HighFive @@ -0,0 +1 @@ +Subproject commit 7340d1fcc7c455519c113216d8415e003ed37cc6 diff --git a/inputs/cooling_tables/Visu.ipynb b/inputs/cooling_tables/Visu.ipynb new file mode 100644 index 00000000..ecd18ae4 --- /dev/null +++ b/inputs/cooling_tables/Visu.ipynb @@ -0,0 +1,102 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "9ceae106-4747-4b78-82d9-7a34f3eb8066", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "id": "acdd51e4-bd01-430a-a6f6-96e9ed05b478", + "metadata": { + "tags": [] + }, + "outputs": [ + { + "data": { + "text/plain": [ + "(-25.0, -21.0)" + ] + }, + "execution_count": 26, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAk8AAAG2CAYAAABmsmIiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABVEElEQVR4nO3dd1hUV/4/8PcMZSjCIAxSZJBiAcWC2I0RggXLqjFr25g1EdlNMWvdJKZZslk3kcSfMd9YUjSWrJuNZWNiEls0UZMoWEFBigVnREBhhuZQ5v7+QCaigHORYQrv1/PME+fOuXc+J9eEN+eee65EEAQBRERERGQUqbkLICIiIrImDE9EREREIjA8EREREYnA8EREREQkAsMTERERkQgMT0REREQiMDwRERERicDwRERERCQCwxMRERGRCAxPRERERCJYZXi6fPky4uPjERwcDGdnZ4SGhmLx4sWoqKio027OnDmIioqCTCZDr169jDq2TqfDiy++CIVCAVdXV4wbNw7Xrl0zQS+IiIjIGllleEpLS4Ner8e6deuQmpqKlStXYu3atXj11VfrtBMEATNnzsSUKVOMPvbcuXOxc+dObNu2DUeOHEFJSQnGjh2L6urq5u4GERERWSGJrTwYeMWKFVizZg2ys7Pv+2zJkiXYtWsXTp8+3egxNBoNvL29sXnzZkPgUqvVUCqV2LNnD0aOHGmK0omIiMiK2Ju7gOai0Wjg6en5UMdITk5GZWUlRowYYdjm7++PiIgIHDt2rMHwpNPpoNPpDO/1ej1u3boFLy8vSCSSh6qJiIiIWoYgCCguLoa/vz+k0oYvztlEeMrKysLq1avx3nvvPdRxcnNz4ejoiLZt29bZ7uPjg9zc3Ab3W758OZYuXfpQ301ERESWIScnBwEBAQ1+blHhacmSJQ8MISdOnECfPn0M79VqNeLi4jBp0iTMmjXLJHUJgtDoCNKiRYswf/58w3uNRoPAwEDk5OTA3d3dJDURERFR89JqtVAqlXBzc2u0nUWFp9mzZ2Pq1KmNtgkKCjL8Wa1WIyYmBgMHDsT69esf+vt9fX1RUVGBwsLCOqNPeXl5GDRoUIP7yWQyyGSy+7a7u7szPBEREVmZB025sajwpFAooFAojGqrUqkQExODqKgobNiwodFrk8aKioqCg4MD9u3bh8mTJwMArl+/jpSUFLz77rsPfXwiIiKyfla5VIFarUZ0dDSUSiUSExORn5+P3Nzc++YlZWZm4vTp08jNzUV5eTlOnz6N06dPG9aDUqlUCAsLw/HjxwEAcrkc8fHxWLBgAQ4cOIBTp05h+vTp6N69O4YNG9bi/SQiIiLLY1EjT8bau3cvMjMzkZmZed+ErrtXXpg1axYOHz5seB8ZGQkAuHTpEoKCglBZWYn09HSUlZUZ2qxcuRL29vaYPHkyysvLERsbi40bN8LOzs7EvSIiIiJrYDPrPFkSrVYLuVwOjUbDOU9ERERWwtif31Z52Y6IiIjIXBieiIiIiERgeCIiIiISgeGJiIiISASGJyIiIiIRGJ6IiIiIRGB4IiIiIhKB4YmIiIhIBIYnIiIiIhEYnoiIiIhEYHgiIiIiEoHhiYiIiEgEhiciIiIiERieiIiIiERgeCIiIiISgeGJiIiISASGJyIiIiIRGJ6IiIiIRGB4IiIiIhKB4YmIiIhIBIYnIiIiIhEYnoiIiIhEYHgiIiIiEoHhiYiIiEgEhiciIiIiERieiIiIiERgeCIiIiISgeGJiIiISASGJyIiIiIRGJ6IiIiIRGB4IiIiIhKB4YmIiIhIBIYnIiIiIhEYnoiIiIhEYHgiIiIiEsEqw9Ply5cRHx+P4OBgODs7IzQ0FIsXL0ZFRUWddnPmzEFUVBRkMhl69epl1LGjo6MhkUjqvKZOnWqCXhAREZE1sjd3AU2RlpYGvV6PdevWoWPHjkhJSUFCQgJKS0uRmJhoaCcIAmbOnInffvsNZ8+eNfr4CQkJWLZsmeG9s7Nzs9ZPRERE1ssqw1NcXBzi4uIM70NCQpCeno41a9bUCU8ffPABACA/P19UeHJxcYGvr2/zFUxEREQ2wyov29VHo9HA09OzWY61detWKBQKdOvWDQsXLkRxcXGj7XU6HbRabZ0XERER2SarHHm6V1ZWFlavXo333nvvoY/15JNPIjg4GL6+vkhJScGiRYtw5swZ7Nu3r8F9li9fjqVLlz70dxMREZHls6iRpyVLltw3WfveV1JSUp191Go14uLiMGnSJMyaNeuha0hISMCwYcMQERGBqVOn4quvvsL+/ftx8uTJBvdZtGgRNBqN4ZWTk/PQdRAREZFlsqiRp9mzZz/wzragoCDDn9VqNWJiYjBw4ECsX7/eJDX17t0bDg4OyMjIQO/evettI5PJIJPJTPL9REREZFksKjwpFAooFAqj2qpUKsTExCAqKgobNmyAVGqaQbTU1FRUVlbCz8/PJMcnIiIi62JRl+2MpVarER0dDaVSicTEROTn5yM3Nxe5ubl12mVmZuL06dPIzc1FeXk5Tp8+jdOnTxvWg1KpVAgLC8Px48cB1MydWrZsGZKSknD58mXs2bMHkyZNQmRkJAYPHtzi/SQiIiLLY1EjT8bau3cvMjMzkZmZiYCAgDqfCYJg+POsWbNw+PBhw/vIyEgAwKVLlxAUFITKykqkp6ejrKwMAODo6IgDBw5g1apVKCkpgVKpxJgxY7B48WLY2dm1QM+IiIjI0kmEu9MGNQutVgu5XA6NRgN3d3dzl0NERERGMPbnt1VetiMiIiIyF4YnIiIiIhEYnoiIiIhEYHgiIiIiEoHhiYiIiEgEhiciIiIiERieiIiIiERgeCIiIiISgeGJiIiISASGJyIiIiIRGJ6IiIiIRGB4IiIiIhKB4YmIiIhIBIYnIiIiIhEYnoiIiIhEYHgiIiIiEoHhiQwEQTB3CURERBbP3twFkPll3CjGG/9LQfKVQoR6t0FkYFs8Hx0KpaeLuUsjIiKyOAxPrdyaQ1l4f186KqtrRp3ScouRlluMny7m46vnBsJP7mzmComIiCwLL9u1YslXbuGd79NQWS1gWHg7fPPiI1g7PQohCleoisrx1KfHUVhaYe4yiYiILArDUyv2wYFMAMAfowLw8Z/7IKK9HHERvtgU3w9+cidk5pXgr5uTUa3nXCgiIqJaDE+t1JmcIhy+mA87qQQvPtYREonE8FlAWxdsju+HNjJ7HL98CxuOXjJjpURERJaF4amVWn0wAwAwoVd7dPByve/zju3c8OrocADAih/SkZVf0qL1ERERWSqGp1bownUt9l/Ig1QCvBAT2mC7af2UGNJJAV2VHn//7xnoefmOiIiI4ak1+vqMGgAwoqsvQrzbNNhOIpHgX0/0QBuZPU5eLcKOU6qWKpGIiMhiMTy1QvvP3wAAjOru+8C27T2c8eJjHQEAiT+ko6yiyqS1ERERWTqGp1bmckEpMvJKYC+VILpLO6P2mTEoCAFtnZGrvY2Pf+LkcSIiat0YnlqZ/RdqRp36h3hC7uxg1D5ODnZ4OS4MALDupyzkaW+brD4iIiJLx/DUyuy7c8luWLiPqP3G9vBDZKAHyiqq8eGPmaYojYiIyCowPLUihaUVSLpSCEB8eJJIJPj7yC4AgG3Hc6AuKm/2+oiIiKwBw1Mr8mN6Hqr1AsJ83Zr00N9BoQoMCPFERbUe/8fRJyIiaqUYnlqRwxfzAQDDu4obdbrbvGGdAQBfJuUg51ZZs9RFRERkTRieWpGTV2su2fUP9mryMfqHeOGRjgpUVgscfSIiolaJ4amVuFmiQ86tckgkQA+l/KGONW94JwDAV8nXcPUmR5+IiKh1YXhqJU7nFAEAOnq3gbuTcUsUNCSqgyce7eyNKr1geEYeERFRa2GV4eny5cuIj49HcHAwnJ2dERoaisWLF6OiosLQ5syZM5g2bRqUSiWcnZ0RHh6OVatWPfDYOp0OL774IhQKBVxdXTFu3Dhcu3bNlN1pEaeuFgEAeik9muV484bVjD7tOKXCpYLSZjkmERGRNbDK8JSWlga9Xo9169YhNTUVK1euxNq1a/Hqq68a2iQnJ8Pb2xtbtmxBamoqXnvtNSxatAgffvhho8eeO3cudu7ciW3btuHIkSMoKSnB2LFjUV1dbepumVTtyFOvQI9mOV5kYFs8FtYO1XoBqw9w9ImIiFoPiSAIgrmLaA4rVqzAmjVrkJ2d3WCbF154ARcuXMDBgwfr/Vyj0cDb2xubN2/GlClTAABqtRpKpRJ79uzByJEjjapFq9VCLpdDo9HA3d1dfGeamV4voOfSvSjWVWHP34agq3/z1HT2WhHGfXgUUgmwb/5QhDbykGEiIiJLZ+zPb6sceaqPRqOBp6fnQ7VJTk5GZWUlRowYYdjm7++PiIgIHDt2rNlqbWlZ+SUo1lXB2cEOnX2aL+D0CPDAsHAf6AXgA44+ERFRK2ET4SkrKwurV6/Gs88+22CbX375BV9++SX++te/NtgmNzcXjo6OaNu2bZ3tPj4+yM3NbXA/nU4HrVZb52VJauc7dQ+Qw96ueU/53Dtzn74+o0bGjeJmPTYREZElsqjwtGTJEkgkkkZfSUlJdfZRq9WIi4vDpEmTMGvWrHqPm5qaivHjx+PNN9/E8OHDRdclCAIkEkmDny9fvhxyudzwUiqVor/DlE7dme8U2Uzzne4W0V6Okd18IAjAKo4+ERFRK2Bv7gLuNnv2bEydOrXRNkFBQYY/q9VqxMTEYODAgVi/fn297c+fP4/HHnsMCQkJeP311xs9tq+vLyoqKlBYWFhn9CkvLw+DBg1qcL9FixZh/vz5hvdardaiAtSpO4tjRjbTnXb3mjusM35IvYFvz13H7FwtwnzNP8+LiIjIVCwqPCkUCigUCqPaqlQqxMTEICoqChs2bIBUev8gWmpqKh577DHMmDEDb7/99gOPGRUVBQcHB+zbtw+TJ08GAFy/fh0pKSl49913G9xPJpNBJpMZVXdLq6jSIzOvBADQPcDDJN8R7ueOMd398O2561i1PwNrpkeZ5HuIiIgsgUVdtjOWWq1GdHQ0lEolEhMTkZ+fj9zc3DrzklJTUxETE4Phw4dj/vz5hs/z8/MNbVQqFcLCwnD8+HEAgFwuR3x8PBYsWIADBw7g1KlTmD59Orp3745hw4a1eD+bQ3ZBCar0Atyc7OEvdzLZ98wZ1gkSCfBdSi5S1RqTfQ8REZG5WWV42rt3LzIzM3Hw4EEEBATAz8/P8Kr13//+F/n5+di6dWudz/v27WtoU1lZifT0dJSV/f6IkZUrV2LChAmYPHkyBg8eDBcXF+zevRt2dnYt2sfmkp5bM4m7i49bo/O2HlZnHzeM7eEPAFi1v3XNfdKUV+K8Wou84tvmLoWIiFqAzazzZEksaZ2nd75Pw5pDWfhT/0D88/HuJv2uzLwSjFh5GHoB2D37EXQPeLhn6Fm6H1Jz8eb/UnBDqwMA2EklGB7ug2cGB6F/SNMfvkxERObR6tZ5ovpdvDPyFObrZvLv6tiuDcb3ag8A+Nf3F2CruVyvF7By30X8dXOyITh5uDigWi/g+9RcTFn/KxJ/SLfZ/hMRtXYMTzYu7a7Ldi1h/vDOcLST4mjmTRxKz3/wDlbo7T0XDMsyPDM4CKlLR+L0myPww9xHMblPAADgwx8zMWfbaeiqrPuxPkREdD+GJxtWfLsSqqJyAECXFhh5AgClpwueGRwEoCZkVFXrW+R7W8qJy7fw2dFLAIB3nuiOxX/oBldZzU2rXXzd8O4fe+LdJ3rAXirB12fUmP/lGY5AERHZGIYnG3bxRs0SBT7uMni4OLbY9z4f0xFtXRyQmVeC/yTltNj3mtrtymq8vP0sBAGY0keJKX0D6203ua8Snz7dF/ZSCb49ex3/92NmC1dKRESmxPBkwwx32rXwopVyZwfMia15bMvKfRdRfLuyRb/fVFYfzEB2finaucnw6pjwRtsO7eyNtyZEAAAS917E3tSGH+9DRETWheHJhqXn1jxjryUmi9/ryQEdEKxwRUFJBdYezmrx729uecW38fFPNZfrlo2PgNzZ4YH7TOsXiBkDOwAAFnx5Bjm3yh6wBxERWQOGJxuWfudBvZ1baLL43RzspHhlVBgA4JOfL0F9Z+6Vtdp07AoqqvXoHeiBuAhfo/d7fWxX9A70QLGuCn/bdgqVNjYHjIioNWJ4slGCIBgu25lj5AkARnT1Qb9gT+iq9Ej8Id0sNTSHUl0VNv96BQDwl0dDRe3rYCfFqqmRcHOyx6mrRVi576IpSiQiohbE8GSj8ot1KCyrhFRSs/6SOUgkErx+Z27QjlMqJF8pNEsdD+vLpBxoyisR5OWC4V19RO+v9HTBO0/0AACsOZyFpMu3mrtEIiJqQQxPNqr2YcAdvFzh5GC+R8v0CPDApKiatY/e2JVidUsXVFXr8emRmrlOs4aEwE7atEfcjO7uh0lRARAE4KWvzuJ2Jdd/IiKyVgxPNiqroBQAEKJwNXMlwCujwiB3dsD561rD5S9rcSg9H9cKy+Hp6og/3gmBTfX6mK5o5yZDdkEp/l8re/4fEZEtYXiyUdn5NSNPId7mD09ebWR4Ka4LAOD9vReRp7WeB+huP3kNADAxsv1Dj+DJXRzw9p3nC67/KQvnrmkeuj4iImp5DE82Kjv/zsiTt3nmO91rat9A9AyQo1hXhX/uuWDucoxSWFqBAxfyAABPPOSoU63hXX3wh57+0AvA67vOoVrP1ceJiKwNw5ONyi64M/JkAZftAMBOKsFbEyIgkQC7TqtxLKvA3CU90O6zalRU69HVzx3hfs230OgbY8LhJrPHmWsa/Pv41WY7LhERtQyGJxt0u7Ia1wpr1lWylJEnoGby+PT+NYtGvvm/VFRUWfbk8e3JNZfsmmvUqVY7dyfMH9EZALDih3QUlOia9fhERGRaDE826MrNMggC4OZkD0WblnumnTEWjugCL1dHZOaV4JMj2eYup0GZecU4c00De6kE43v5N/vxnxrQAV393KEpr8S/vktr9uMTEZHpMDzZoN8ni7eBRNK0W+tNRe7igFdH16z9tGp/Bq7cLDVzRfXbdUoNAIju4g1FG1mzH9/eTop/PF7z7Luvkq/hBNd+IiKyGgxPNij7zjIFoRYy3+leE3u3x+COXtBV6fHazhQIgmVNmhYEAXvOXQcA/KFn84861eod2BZT+yoBWOcaWERErRXDkw36/U47ywxPEokEb0/oDpm9FEcyC7DjpMrcJdWRfqMY2QWlcLSXIjZc/IriYrwcF4a2Lg5Iyy3GxmOXTfpdRETUPBiebJDhTjsLmix+ryCFK+YM6wQA+Me353HTgiZN7zlbM+o0tLM32sjsTfpdbV0d8XJczQOUV+67iFyN9ayBRUTUWjE82RhBECx+5KlWwpAQhPm6obCsEv/41nLWftqTkgsAGN3dt0W+b3IfJSIDPVBaUY23vj3fIt9JRERNx/BkY26VVkBTXgmJBAjysuzw5GAnxb+e6AGJBNh5SoWfLuabuyRk3ChGZl4JHO1Mf8mullQqwT8mREAqAb49ex0/Z5j/3wMRETWM4cnG1E4W95c7m/WBwMbqpfTA04OCAACv7TqH8grzPjD32zsTxYd0UsDdyaHFvrebvxx/HhgEoGYNLF0VHxxMRGSpGJ5sjCU9085YC0Z0gb/cCTm3yrHqgHkfmPvduZpLdqO6+7X4d88f0RnebjJcKijFxz9Z7hpYREStHcOTjblUUAbAch7LYow2MnssG1+z5tEnP2cjLVdrljoy80qQfqMYDnYSDG+hS3Z3c3dywOtjatbAWn0wEzm3ylq8BiIiejCGJxtz9VbNZbsOFj7f6V7DuvogrpsvqvQCFu04B70ZHpj73Z1LdoM7KiB3ablLdncb19MfA0Nq1sBa/HWqxa2BRUREDE8258rNmtGKDl4uZq5EvCXjuqGNzB6nrhbhCzM8MNdwl11Ey1+yqyWRSPDWhG5wsJPgYFoe9p2/YbZaiIiofgxPNkQQBFy14vDkK3fCwjsPzH3n+zTkaVtuzaNLBaW4cF0Le6kEI7q1/CW7u3Vs54aEISEAgKW7z6Ososqs9RARUV0MTzbkVmkFinVVkEiAgLbWF54A4KmBQegZIEfx7Sos/abl1jyqfRzLwFAveLiY/2HKLz7WCe09nKEqKseHBzPNXQ4REd2F4cmGXLkzwdjX3ckqlimoj51Ugn9O7A47qQTfnr2OH9PyWuR7v0upCU9jzHCXXX2cHe2wZFw3AMDHP2cjM6/YzBUREVEthicbUnvJLtDTOkedanXzl2Pm4CAAwOu7UlCiM+1lq+z8EqSotLCTSjC8q3kv2d1teFcfxIa1Q2W1gNd2pphlEj0REd2P4cmGWPNk8XvNG94ZSs+ay1bL95j20S27TtU8mHhIJwW82shM+l1iLRnXDc4Odvjt0i38JynH3OUQEREYnmzKFStdpqA+Lo72eOeJHgCArb9dxbHMApN8jyAI2Hm6Jjw9HtneJN/xMJSeLlhwZxL9P/dcwI0WnERPRET1Y3iyIbZy2a7WoFAFnuwfCAB4ecdZlJrg8l3ylULk3CqHq6MdRnRtmQcBi/XM4GDDJPo3/5di7nKIiFo9hicbUjth3BYu29VaNDoc7T2ckXOrHCt+SG/24++4c8kuLsIPzo6WOcneTirBv57oAXupBD+k3jAs5klERObB8GQjyiqqkF+sAwB08LT+y3a12sjssXxidwDAxmOXcfzSrWY7tq6qGt+erQkiE3tb3iW7u4X7ueO56FAAwJtfp0JTVmnmioiIWi+rDE+XL19GfHw8goOD4ezsjNDQUCxevBgVFRWGNmfOnMG0adOgVCrh7OyM8PBwrFq16oHHjo6OhkQiqfOaOnWqKbvTLGoni8udHcz2aBFTebSzNyb3CQAAvPTVGZRXVDfLcfedvwFNeSV83Z0wIMSrWY5pSi/EdESItyvyi3X4p4kn0RMRUcOsMjylpaVBr9dj3bp1SE1NxcqVK7F27Vq8+uqrhjbJycnw9vbGli1bkJqaitdeew2LFi3Chx9++MDjJyQk4Pr164bXunXrTNmdZmFLd9rV57UxXeHr7oTLN8vwr++aJzh8duQSAGByXyXspJJmOaYpOTnYGSbR/ycpB0cyTDOJnoiIGmdv7gKaIi4uDnFxcYb3ISEhSE9Px5o1a5CYmAgAmDlzZp19QkJC8Msvv2DHjh2YPXt2o8d3cXGBr69lTh5uiLU+ENhYcmcHvPPHHpjx2XF8/ssVxIS1Q3SXdk0+3smrhTh5tQiOdlI8NaBDM1ZqWn2DPPHngR2w6ZcreOmrM/h+3qNwd7KtkUYiIktnlSNP9dFoNPD09HzoNgCwdetWKBQKdOvWDQsXLkRxceOrO+t0Omi12jqvlmYYebKRO+3qM7SzN54eFAQA+PtXZ3GrtKLxHRrx6Z1Rp3G9/OHtZllrOz3Iy3Fh6ODlArXmNt7a3XKPsCEioho2EZ6ysrKwevVqPPvssw22+eWXX/Dll1/ir3/9a6PHevLJJ/Hvf/8bhw4dwhtvvIHt27dj4sSJje6zfPlyyOVyw0upVDapHw/j6p077QJt9LJdrVdGhaFjuzbIL9Zh7n9Oo7oJq26risrxfUouAGDm4ODmLtHkXGX2SJzUExIJ8N/ka9h//oa5SyIialUsKjwtWbLkvsna976SkpLq7KNWqxEXF4dJkyZh1qxZ9R43NTUV48ePx5tvvonhw4c3WkNCQgKGDRuGiIgITJ06FV999RX279+PkydPNrjPokWLoNFoDK+cnJZfCfpaYTkAQGmlDwQ2lpODHT6YGgknByl+upjfpOUL1h7KQrVewKBQL3T1dzdBlabXN8gTCUNCAACv7DiHwocYhSMiInEsas7T7NmzH3hnW1BQkOHParUaMTExGDhwINavX19v+/Pnz+Oxxx5DQkICXn/9ddE19e7dGw4ODsjIyEDv3r3rbSOTySCTme/Sj14vQFUbnjydzVZHS+nq7453/9gTf/v3Kaw9nIVwPzeM72XcUgOnc4qw5bcrAIDZMR1NWabJzR/eGQfT8pCZV4I3/peCD/9U/99PIiJqXhYVnhQKBRQKhVFtVSoVYmJiEBUVhQ0bNkAqvX8QLTU1FY899hhmzJiBt99+u0k1paamorKyEn5+fk3avyXkFetQUa2HnVQCX3cnc5fTIsb19EeqWoN1h7Ox8L9n4OxghxHdGp/kX1Wtx6Id5yAINY9iGdTRuL9rlsrJwQ7vTeqJiWuO4Zuz1xEXocbYHv7mLouIyOZZ1GU7Y6nVakRHR0OpVCIxMRH5+fnIzc1Fbm6uoU1qaipiYmIwfPhwzJ8/3/B5fn6+oY1KpUJYWBiOHz8OoGbu1LJly5CUlITLly9jz549mDRpEiIjIzF48OAW76exrhXWzHfykzvB3s4qT2mTvDQyDH/o6Y/KagEvfHESe1NzG23/yZFLuHBdCw8XB7w+JryFqjStnkoPvHBn8czXd6Ugj8++IyIyOav8Sbt3715kZmbi4MGDCAgIgJ+fn+FV67///S/y8/OxdevWOp/37dvX0KayshLp6ekoK6sJH46Ojjhw4ABGjhyJLl264G9/+xtGjBiB/fv3w87OMh/dAQA5d8KTrc93upedVIKVk3saAtRzW0/iw4MZ9U4i3/zrFbzzfRoA4NVR4fBqY1132DVm9mOd0M3fHUVllXhp+1kIgvhJ9EREZDyJwP/TNjutVgu5XA6NRgN3d9NPSF59IAPv7buISVEBWDGpp8m/z9LUXo77b/I1AED/YE9MH9AB/YM9kau9jW/PXse6n7IBAE8N6IBl47tBIrH8RTHFyLhRjDGrj6CiSo9/TIjAdCtau4qIyFIY+/PbouY8UdMY7rSz4TWeGmNvJ8W7f+yB/iFeeGNXCn67dAu/1fMMvL891hHzhne2ueAEAJ183PByXBje+uY83v72AvoHe6KTj5u5yyIisklWedmO6qq9bBfQ1vbvtGuIRCLBH6MCsGfOEMx6JBhhvjXBwcPFAUM6KbBySk/MH9HFJoNTrWcGBeGRjgqUV1bj+a0nUVZRZe6SiIhsEkeebEDtyFNAK5vzVJ9ghSteH9sVAFBWUQVnBzubDkx3k0olWDmlF0Z/8DMy8krw+q4UvDepZ6vpPxFRS+HIk5Wr1gtQF7WeNZ7EcHG0b3XBwdtNhtXTIiGVADtOqvDfpGvmLomIyOYwPFm5XO1tVOkFONhJ0M6tdazxRI0bEOKFBSO6AADe+F8KLlxv+WctEhHZMoYnK3ftzjPt/D2cYSdtXaMs1LDnhoYiuos3dFV6PL/1JIpvV5q7JCIim8HwZOVayzPtSBypVIKVk3vBX+6ESwWld1ZW56okRETNgeHJyvFOO2pIW1dHrP5Tb9hLJfjm7HVs+fWKuUsiIrIJDE9W7vc77Rie6H5RHdrilVFhAIC3vrmAs9eKzFsQEZENYHiycrXPtWutC2TSg8U/EowRXX1QUV0z/6mwtMLcJRERWTWGJyuXc4sjT9Q4iUSCFZN6ooOXC64VluOFL06iqlpv7rKIiJqkslqPtFwtNGXmuxGGi2RasapqPXK1twFwgUxqnNzZAeuf6oPHPzqKY1k3sfy7NLxxZzFRIiJLVFhagVztbeQV65BxoxgXrhfjwnUtMvNKUFGtx+ppkfhDT3+z1MbwZMVytbdRrRfgaCeFdxuZucshC9fF1w3vT+6JZ7ecxKdHLqFjuzaY1i/Q3GUREUEQBFy9VYZUtRZJlwvxU0Y+MvNKGmzvJrOHppwjT9QEqjuTxf08nCDlGk9khLgIP8yJ7YRVBzLw+q4U+MmdEN2lnbnLIqJWRK8XkJFXghSVBqlqLVLVGpy/rkXx7fufx+nl6ghFGxmCFC4I93NHuJ87uvq5I6Cts1mfIMHwZMVUdx7L0t6D853IeHOHdULOrTLsOKXCC1tP4stnB6Kbv9zcZRGRjarWC8jKrwlLv2bfxI/p+cgv1t3XztFOis6+bdC9vQce6ajAIx0VkLs4mKHiB2N4smJqhidqAolEgn890QO52ts4lnUTMzeewM7nB8Off4+IqBkIggBVUTlOXi3Cj2l5OJiWd98lNhdHO0T4y9HV3x3d/N3RzV+OTj5t4GBnHfexMTxZsdqRJ/7QI7Ec7aVYMz0Kk9Yew8UbJXhmwwn897mBcHeyzN/yiMiyleqq8NPFfPyQmoufMwpw854lUVwc7dDN3x09AzwQ3aUd+ga3hczezkzVPjyGJytWu0Bmey5TQE0gd3bAhmf64fH/O4r0G8V4bksyNjzdD4721vGbHxGZT1W1HhdvlOBUTiF+TMvDTxkFqKj6fQkUe6kEXXzdMCjUC8O7+qJ3oAfsrWRUyRgMT1asduQpgCNP1ETtPZzx2dN9MXndLziaeROLdpxD4qQeZp2ISUSW6XZlNX66mI9vz13HgQt5KNHVneAd6OmCkd18MLyrL3oEyOHkYL0jSw/C8GSlBEH4fc4TR57oIUS0l+P/nuyNWZ8nYfvJawho64x5wzubuywisgC6qmr8fLEA3567jn3nb9QJTG1k9uiplKN/sBdGdvNFZ582reYXL4YnK3WrtAK3K2uGSH3lTmauhqxdTJd2eGt8BF7deQ6rDmRA7uyAmY8Em7ssIjKDiio9jmYW4Juz17H3fG6dJQR83Z0wpocfRnf3Qy+lB+xa6TI5DE9WqvaSXTs3mVVPuiPL8af+gbihvY1VBzKw7JvzcLSXYvqADuYui4hagF4vIOlKIXaeUmHPuet17o5r5ybD6O5++ENPP0Qq23JdQTA8WS0VJ4uTCcwd1gm3q6qx7nA2Xt+VAkd7KSb3UZq7LCIykcy8Yuw8pcKuU2rDL+UA4O0mw+gIX4zp4Y8+HRiY7mWS8HT69Gn06tXLFIemO7hAJpmCRCLBK3Fh0FXqsfHYZby8/Sxk9lKM79Xe3KURUTO5VVpxJzCpcE6lMWxvI7PHqAhfTIhsjwEhXq32kpwxmi08aTQabN26FZ988gnOnDmD6urq5jo01YPhiUxFIpFg8R+6oqJajy9+u4r5X56Bo50Uo7r7mbs0Imqiar2AI5kF+M+Jq9h3/gYqqwUANUsKRHfxxoTI9hgW7mPTd8g1p4cOTwcPHsRnn32GHTt2wNnZGVOmTMGaNWuaozZqBC/bkSlJJBL8Y3wEKqr0+Cr5Gl789ymstZNiWFcfc5dGRCLk3CrDf5Ov4aukHKg1tw3bewTI8ceoAIzt4Q9PV0czVmidmhSerl27ho0bN2LDhg3IycnByJEjsWnTJowbNw6OjjwJLYEjT2RqUqkE7zzRA7oqPXafUeP5rSfxyYw+eLSzt7lLI6JG3K6sxt7zN/DliRwczSqAUDPIBLmzAx6PbI/JfZTo6u9u3iKtnOjwNHr0aOzbtw8RERGYPXs2nnzySbRrx6eytzQ1H81CLcBOKsH7k3uioqoaP6TeQMKmJGx8ph8GhnqZuzQiuouqqBw/X8zHTxn5OJJRAO1dyws80lGByX2VGNGVl+Wai+jw9P3336Nbt25YunQpxowZAzs7noiWVlZRhcKymttIedmOTM3BTorV03rj2S3JOJiWh/jPT2DTzH7oE+Rp7tKIWrXrmnJs/fUqvku5jqz80jqf+cmdMCkqAJP6KKH0dDFThbZLdHg6evQoPvvsM0yfPh1OTk6YNm0aZsyYgd69e5uiPqpH7XwnNyd7PsiVWoSjvRQf3VmF/EhmAZ7ZcAL75g/lAq1ELeyG9jYOp+dj34UbOJiWh2p9zTU5qQSIDGyLIZ0UGNJJgV7KtrxbzoQkglB7NVScsrIybNu2DZ999hmOHTuGrl27YsaMGXjyySfh7+/f3HVaFa1WC7lcDo1GA3f35r+u/GN6Hp7ZcAJhvm74fu6jzX58ooaUV1Rj8rpfcE6lQfwjwXhjbFdzl0Rk87LyS7DzpAo/puchVa2t81n/YE88OaADhnb2htyZv0w/LGN/fjc5PN0tPT0dn376KTZv3oyCggLExsbi+++/f9jDWi1Th6etv13BaztTEBvWDp8+3bfZj0/UmMMX8zHjs+NwdrDD0Vce4506RCagq6rGsayb2PLLFRxIyzNsl0iAHgEeiOnijVERfuji62bGKm2PsT+/m2Wdpy5duuDdd9/F8uXLsXv3bnz22WfNcVhqAJcpIHN6tJMC3dvLcU6lwYajl7BgRBdzl0RkEzTllTiUnoe952/gcHq+4SG8EgnwWJd2GNPDD4929oaijczMlVKzrjBuZ2eHCRMmYMKECc15WLoHlykgc5JIJHghJhTPbjmJjccu4y+PhsCNc++ImuRWaQW+PavG3vM38EvWTVTpf78Y1M5NhlERvpgxKAgh3m3MWCXdi8+2s0IceSJzG9HVF6HersjKL8XL289i1dRIONhJzV0WkdW4ob2NT37OxpZfr6K88vcncnRs1wYjuvpgeFcf9Azw4DPlLJRV/t/u8uXLiI+PR3BwMJydnREaGorFixejoqLC0ObmzZuIi4uDv78/ZDIZlEolZs+eDa1W28iRAZ1OhxdffBEKhQKurq4YN24crl27ZuouicI1nsjcpFIJlozrBgc7Cfacy8VzW5Jxu5KPZCJqTH6xDpt/vYKp63/BgOUH8PHPl1BeWY2ufu54dXQYflwYjf3zh+KluDBEBvJhvJbMKkee0tLSoNfrsW7dOnTs2BEpKSlISEhAaWkpEhMTAQBSqRTjx4/HP/7xD3h7eyMzMxMvvPACbt26hS+++KLBY8+dOxe7d+/Gtm3b4OXlhQULFmDs2LFITk62iDWtKqv1yNXWLLEfwPBEZjSkkzc+/nMf/HVzMvZfyEPCpiSsf6oPnB3N/98JkaXQlFfi+5Tr+N9pNX7Nvom7rsqhX5AnnosJRXRnb0gkDErWpFnutrMEK1aswJo1a5Cdnd1gmw8++AArVqxATk5OvZ9rNBp4e3tj8+bNmDJlCgBArVZDqVRiz549GDlypFG1mPJuu5xbZRjy7o9wtJMi7a04/mZCZncsswCzNiWhrKIa/YI98dnTfdFGZpW/lxE1i4oqPQ5fzMeuUyrsu3ADFVV6w2c9A+QY08MPoyL8uHilBTL53Xbz58+vd7tEIoGTkxM6duyI8ePHw9OzZVYh1mg0jX6XWq3Gjh07MHTo0AbbJCcno7KyEiNGjDBs8/f3R0REBI4dO9ZgeNLpdNDpdIb3D7o0+DBUhkt2TgxOZBEGdVRg08x+eGbDCRy/dAtPfvIbPpvRB168I4haEUEQcPJqEXadUuGbs2rDUyAAoFO7NpgQ2R7jevozMNmIJoenU6dO4eTJk6iurkaXLl0gCAIyMjJgZ2eHsLAwfPTRR1iwYAGOHDmCrl1Nu5BeVlYWVq9ejffee+++z6ZNm4b//e9/KC8vxx/+8Ad88sknDR4nNzcXjo6OaNu2bZ3tPj4+yM3NbXC/5cuXY+nSpU3vgAic70SWqE+QJ7Ym9MefPzuOMzlFmLjmGDY+0w/BCldzl0ZkUpcKSrHrlAq7Tqtw5WaZYbu3mwzje/pjQmR7dPN352U5G9PkCePjx4/HsGHDoFarkZycjJMnT0KlUmH48OGYNm0aVCoVHn30UcybN8/oYy5ZsgQSiaTRV1JSUp191Go14uLiMGnSJMyaNeu+Y65cuRInT57Erl27kJWV1eCIWWMEQWj0L/6iRYug0WgMr4YuCzYHw512DE9kYXoEeOCrZwchoK0zrtwsw8SPjiL5yi1zl0XU7K5ryrHh6CU8/tFRxCQewqoDGbhyswwujnaYGNkem2b2w6+LYvH62K6IaC9ncLJBTZ7z1L59e+zbt+++UaXU1FSMGDECKpUKJ0+exIgRI1BQUGDUMQsKCh7YNigoCE5ONc/TUqvViImJQf/+/bFx40ZIpY1nwSNHjmDIkCFQq9Xw8/O77/ODBw8iNjYWt27dqjP61LNnT0yYMMHo0SVTznl6ZftZbDuRg7nDOmHusM7Nemyi5pBfrMOsz0/gzDUNZPZS/L8pvTCq+/3/vRFZk2uFZfg+JRd7zl3HyatFhu1SSc3NE49HtseIbj5wceR8P2tm8jlPGo0GeXl594Wn/Px8w5wfDw+POssHPIhCoYBCoTCqrUqlQkxMDKKiorBhw4YHBiegZgQJQJ35SXeLioqCg4MD9u3bh8mTJwMArl+/jpSUFLz77rtG9sK0VLxsRxbO202Gf/9lAP7271PYfyEPz39xEgtHdMFzQ0M5T4+sivZ2Jb47dx3bk1U4fvn3UVSJBOjToS1GRfhhbE8/tHPjA7JbmyaHp/Hjx2PmzJl477330LdvX0gkEhw/fhwLFy40rDB+/PhxdO7c/KMjarUa0dHRCAwMRGJiIvLz8w2f+fr6AgD27NmDGzduoG/fvmjTpg3Onz+Pl156CYMHD0ZQUBCAmgAWGxuLTZs2oV+/fpDL5YiPj8eCBQvg5eUFT09PLFy4EN27d8ewYcOavR9NURueuEwBWTIXR3use6oPlu1Oxee/XMGKH9Jx6moR3pvckw8vJYtWVa3Hz5kF2HlShR9Sc6G7c6ecRFKztMDo7n6Ii/CFjzsDU2vW5PC0bt06zJs3D1OnTkVVVc3zd+zt7TFjxgysXLkSABAWFtboBO2m2rt3LzIzM5GZmYmAgIA6n9WOLjk7O+Pjjz/GvHnzoNPpoFQqMXHiRLzyyiuGtpWVlUhPT0dZ2e+T/FauXAl7e3tMnjwZ5eXliI2NxcaNGy1ijSdBEAwTxrm6OFk6O6kES8dHINzPHW9+nYr9F25g/IdHsGZ6FML9mv+B2URNJQgCUlRa7Dh1DbvPqFFQ8vsVk47t2uCJ3gGYEOkPPzn/v0s1Hnqdp5KSEmRnZ0MQBISGhqJNGz5/x1RzngpKdOjzj/2QSIC0t+Igszd/oCMyxtlrRXhuy0moisrh5CDFPx/vjom9Ax68I5EJXSssw/9Oq7HzlAqZeSWG7Z6ujvhDDz88ERWA7pzw3aqYfM5TrTZt2qBHjx4PexgyQu2ddt5tZAxOZFV6BHjgmxcfwdz/nMbhi/mY/+UZHMkswLLxEVxQk1qUpqwS36Vcx45TKhy/9Ps8Jpm9FMO7+uDxyPZ4tLM3n9VIjXqo/2sVFRXh008/xYULFyCRSBAeHo74+HjI5fLmqo/uwkt2ZM3aujris6f74sODmVh14CJ2nFQh+Uoh3n2iB/qHeJm7PLJh2fklOHAhDwfSbuDE5UJU33lGikQCDAj2wuO92yMuwhfuTpyPR8ZpcnhKSkrCyJEj4ezsjH79+kEQBKxcuRL//Oc/sXfvXvTu3bs56yT8PlmcazyRtbKTSjBnWCcM6uiFudtO48rNMkxZ/yv+1D8Qr4wK4w8vahaV1XqcuHwLBy7k4WBaHi4VlNb5vIuPGyZEtsf4Xv68c5mapMnhad68eRg3bhw+/vhj2NvXHKaqqgqzZs3C3Llz8dNPPzVbkVTjWiFHnsg29A3yxHdzh2D5njT8+/hVfPHbVRy4cAPLxkdgZDdfc5dHVqiwtAKHLuZh/4U8/JSej2JdleEzBzsJBoR44bGwdngsrB06eHHle3o4DzXydHdwAmrutnvppZfQp0+fZimO6nJxtEN7D2co2/LZSGT93J0csHxid4zv5Y9FO87hUkEp/ro5GaMifLF0XDe0463g1AhBEJCRd+dy3IUbOHm1EPq7bn/ycnVETFg7xIa1wyOdFHDjqCY1oyaHJ3d3d1y9ehVhYWF1tufk5MDNze2hC6P7vRQXhpfiwh7ckMiKDAjxwndzhuCDAxlY91M2vkvJxZHMArw6OhxT+ii5sCYZVFTp8dulm4b5Szm3yut8HubrhmHhPngsvB16BnjAjn93yESaHJ6mTJmC+Ph4JCYmYtCgQZBIJDhy5Aj+/ve/Y9q0ac1ZIxHZOCcHO7wUF4axPfzx8vazOKfSYNGOc/jit6t48w9d0TfI09wlkpkUllbgx/Q87L9wAz9dLEDJXZfjHO2lGBTqhdhwHzwW1o7zQanFNHmdp4qKCvz973/H2rVrUVVVBUEQ4OjoiOeeew7/+te/IJPJmrtWq2HKZ9sR2bqqaj02HruM/7c/w/CDckx3P7wyKgxKT16ybg2y8ktw4MIN7D+fh6Qrt+pcjlO0kSE2rB1iw2sux/FZctScjP35/dCLZJaVlSErKwuCIKBjx45wceH/3BieiB5efrEO7++7iP+cuAq9UDPKMOuRYDwf05FrQ9kYXVU1jl+6hUPp+fgxLQ/Z99wdV3s5blhXH/RoL+elXDIZk4Sn+fPnG13A+++/b3RbW8PwRNR8zqu1eOub8/gl+yaAmpGHv4/sjD9GKTmnxYqpi8pxIC0Ph9PzcDTzJsorqw2f1d4dNyzcB7Hh7RDAm2SohZgkPMXExBjVTiKR4ODBg8Ye1uYwPBE1L0EQsO/8Dby95wKu3Kx5FmVXP3e8MbYrBoZygU1roNcLOKfSYP+FG9h/IQ8XrmvrfN7OTYahnb0R3aUdHu3Mu+PIPFrssh3dj+GJyDR0VdXYdOwKPjiQYVjHJ6aLNxaM6IKI9nyygaUp1VXhWNZNHLhwAwfS8pBfrDN8JpUAvQPbIiasHaK7eKOrnzufIUdmx/BkRgxPRKZ1s6RmPtS2EzmGR22M6eGH+cM7I9SbDyc3l9q1lw6l5+HwxXycuFSIimq94XNXRzsM7eKN2DAfxIS1g6eroxmrJbofw5MZMTwRtYxLBaVYue8ivj6jBlDz+JfHI9vjhZiOCFZwFemWUHy7Ekczb+LwxTwcTs+HWnO7zudKT2fEdGmHYeE+6B/iyYeak0VjeDIjhieilnXhuhbv7U3H/gt5AGouCY3u7ocn+3fAgBBPXg5qRtV6ASkqDY5kFuCni/lIvlKIqrvWEpDZSzEgxOvO/CVvBCtc+e+frAbDkxkxPBGZx6mrhfjwYCYOpOUZtgUrXDGlrxJP9A6At1vrXX/uYVwrLMPPGQU4klGAo1kFKCqrrPN5sMIVQzt7Y2gXbwwI9oKzI0eXyDoxPJkRwxOReaWqNdjy6xV8fVqN0oqaW+DtpRJEd/HGH3r6Y1i4D1y5VlSD8opv49TVIhzNLMDPGQW4dM+6S24yewwM9cKQzt54tJOCD9olm8HwZEYMT0SWoVRXhd1n1Nh2Igenc4oM250cpIgN98HY7n4Y3EkB91Z8W3zx7UqcU2lwJkeDMzlFOHut6L55S3ZSCSKVHnikkwJDOnmjZ4Ac9nZSM1VMZDoMT2bE8ERkeS7eKMbuM2rsPqPG5TtrRQE1waCX0gND7gSDHgFyONhgMBAEATe0OmTnlyAzv6QmLF0rQlZ+Ce79KSCVAJ3auaF/iCce6ajAgFCvVh0wqfVgeDIjhiciyyUIAlJUWnx9RoUDF+5/FIijvRRd/dzRvb0c3QPk6BEgR0fvNlYz0lJeUY1LBaXILihBVl7NP7PzS5GdX2K4hHmv9h7O6KX0QE+lHD0DPBDRXs7LmtQqMTyZEcMTkfW4VliGIxk1c3uOZBZAU155Xxsnh98DVQcvV3i7yaBoI4O3mwxero5wc7J/6HBVrRegKa9EUVkFCstq/llUVonCu/9553NteRX85E6IDGwLV5kdsvNLkZVfE5JUReUNfoedVIJATxeEKFzRrb0cPQPk6BHgwYn0RHcwPJkRwxORddLrBVy9VYazKg3OXSvCOZUGKSotSu6sZt4YJwcp2sgc4OZkD5m9FPZ2EthJpXCQSmAnlUAvCKjSC9DrBVQLAqr1Nd93u6oaRWWV0N6uvO/yWVN5uDggROGKUO82CPFugxBvV4R6uyLQ0xWO9tYxgkZkDgxPZsTwRGQ79HoBl26WIkWlwblrGlzX3kZ+sQ4FJTrkF+tQfPvBwUoMN5k9PFwd4OHsCA8XB7R1qfmnh4sj2ro4wMPFAW4yB1wqKMXpnCLoqvQI9Xa9E5BqwhJX7iZqGmN/fvOiNhFRI6RSCUK92yDUuw3G92p/3+cVVXqU6qpQoqtC8e2af1ZU6VGp16OqWkBVtR7VggA7iQRSqQR2kpqRqNqXo730TihyhNzZwSYnqxPZGoYnIqKH4GgvhaO9I9pytIeo1eCvOEREREQiMDwRERERicDwRERERCQCwxMRERGRCAxPRERERCIwPBERERGJwPBEREREJALDExEREZEIDE9EREREIjA8EREREYlgleHp8uXLiI+PR3BwMJydnREaGorFixejoqLC0ObmzZuIi4uDv78/ZDIZlEolZs+eDa1W2+ixo6OjIZFI6rymTp1q6i4RERGRlbDKZ9ulpaVBr9dj3bp16NixI1JSUpCQkIDS0lIkJiYCAKRSKcaPH49//OMf8Pb2RmZmJl544QXcunULX3zxRaPHT0hIwLJlywzvnZ2dTdofIiIish5WGZ7i4uIQFxdneB8SEoL09HSsWbPGEJ7atm2L5557ztCmQ4cOeP7557FixYoHHt/FxQW+vr7NXzgRERFZPau8bFcfjUYDT0/PBj9Xq9XYsWMHhg4d+sBjbd26FQqFAt26dcPChQtRXFzcaHudTgetVlvnRURERLbJJsJTVlYWVq9ejWefffa+z6ZNmwYXFxe0b98e7u7u+OSTTxo91pNPPol///vfOHToEN544w1s374dEydObHSf5cuXQy6XG15KpfKh+kNERESWSyIIgmDuImotWbIES5cubbTNiRMn0KdPH8N7tVqNoUOHYujQofUGo9zcXBQVFSE9PR2vvvoqhg4dio8++sjompKTk9GnTx8kJyejd+/e9bbR6XTQ6XSG91qtFkqlEhqNBu7u7kZ/FxEREZmPVquFXC5/4M9viwpPBQUFKCgoaLRNUFAQnJycANQEp5iYGPTv3x8bN26EVNr4QNqRI0cwZMgQqNVq+Pn5GVWTIAiQyWTYvHkzpkyZYtQ+xv7LJyIiIsth7M9vi5owrlAooFAojGqrUqkQExODqKgobNiw4YHBCagJQgDqjBI9SGpqKiorK40OW0RERGTbLCo8GUutViM6OhqBgYFITExEfn6+4bPau+T27NmDGzduoG/fvmjTpg3Onz+Pl156CYMHD0ZQUBCAmgAWGxuLTZs2oV+/fsjKysLWrVsxevRoKBQKnD9/HgsWLEBkZCQGDx5sjq4SERGRhbHK8LR3715kZmYiMzMTAQEBdT6rHV1ydnbGxx9/jHnz5kGn00GpVGLixIl45ZVXDG0rKyuRnp6OsrIyAICjoyMOHDiAVatWoaSkBEqlEmPGjMHixYthZ2fXch0kIiIii2VRc55sBec8ERERWR9jf37bxFIFRERERC2F4YmIiIhIBIYnIiIiIhEYnoiIiIhEYHgiIiIiEoHhiYiIiEgEhiciIiIiERieiIiIiERgeCIiIiISgeGJiIiISASGJyIiIiIRGJ6IiIiIRGB4IiIiIhKB4YmIiIhIBIYnIiIiIhEYnoiIiIhEYHgiIiIiEoHhiYiIiEgEhiciIiIiERieiIiIiERgeCIiIiISgeGJiIiISASGJyIiIiIRGJ6IiIiIRGB4IiIiIhKB4YmIiIhIBIYnIiIiIhEYnoiIiIhEYHgiIiIiEoHhiYiIiEgEhiciIiIiERieiIiIiERgeCIiIiISgeGJiIiISASGJyIiIiIRrDI8Xb58GfHx8QgODoazszNCQ0OxePFiVFRU1Nv+5s2bCAgIgEQiQVFRUaPH1ul0ePHFF6FQKODq6opx48bh2rVrJugFERERWSOrDE9paWnQ6/VYt24dUlNTsXLlSqxduxavvvpqve3j4+PRo0cPo449d+5c7Ny5E9u2bcORI0dQUlKCsWPHorq6ujm7QERERFZKIgiCYO4imsOKFSuwZs0aZGdn19m+Zs0a/Oc//8Gbb76J2NhYFBYWwsPDo95jaDQaeHt7Y/PmzZgyZQoAQK1WQ6lUYs+ePRg5cqRRtWi1Wsjlcmg0Gri7uz9Uv4iIiKhlGPvz2ypHnuqj0Wjg6elZZ9v58+exbNkybNq0CVLpg7uanJyMyspKjBgxwrDN398fEREROHbsWIP76XQ6aLXaOi8iIiKyTTYRnrKysrB69Wo8++yzhm06nQ7Tpk3DihUrEBgYaNRxcnNz4ejoiLZt29bZ7uPjg9zc3Ab3W758OeRyueGlVCqb1hEiIiKyeBYVnpYsWQKJRNLoKykpqc4+arUacXFxmDRpEmbNmmXYvmjRIoSHh2P69OkPXZcgCJBIJA1+vmjRImg0GsMrJyfnob+TiIiILJO9uQu42+zZszF16tRG2wQFBRn+rFarERMTg4EDB2L9+vV12h08eBDnzp3DV199BaAmAAGAQqHAa6+9hqVLl953bF9fX1RUVKCwsLDO6FNeXh4GDRrUYE0ymQwymeyB/SMiIiLrZ1HhSaFQQKFQGNVWpVIhJiYGUVFR2LBhw31zmrZv347y8nLD+xMnTmDmzJn4+eefERoaWu8xo6Ki4ODggH379mHy5MkAgOvXryMlJQXvvvtuE3tFREREtsSiwpOx1Go1oqOjERgYiMTEROTn5xs+8/X1BYD7AlJBQQEAIDw83HC3nUqlQmxsLDZt2oR+/fpBLpcjPj4eCxYsgJeXFzw9PbFw4UJ0794dw4YNa5nOERERkUWzyvC0d+9eZGZmIjMzEwEBAXU+E7PyQmVlJdLT01FWVmbYtnLlStjb22Py5MkoLy9HbGwsNm7cCDs7u2arn4iIiKyXzazzZEm4zhMREZH1aXXrPBERERG1BIYnIiIiIhEYnoiIiIhEYHgiIiIiEoHhiYiIiEgEhiciIiIiERieiIiIiERgeCIiIiISgeGJiIiISASGJyIiIiIRGJ6IiIiIRGB4IiIiIhKB4YmIiIhIBIYnIiIiIhEYnoiIiIhEYHgiIiIiEoHhiYiIiEgEhiciIiIiERieiIiIiERgeCIiIiISgeGJiIiISASGJyIiIiIRGJ6IiIiIRGB4IiIiIhKB4YmIiIhIBIYnIiIiIhEYnoiIiIhEYHgiIiIiEoHhiYiIiEgEhiciIiIiERieiIiIiERgeCIiIiISgeGJiIiISASGJyIiIiIRGJ6IiIiIRLDK8HT58mXEx8cjODgYzs7OCA0NxeLFi1FRUVFv+5s3byIgIAASiQRFRUWNHjs6OhoSiaTOa+rUqSboBREREVkje3MX0BRpaWnQ6/VYt24dOnbsiJSUFCQkJKC0tBSJiYn3tY+Pj0ePHj2gUqmMOn5CQgKWLVtmeO/s7NxstRMREZF1s8rwFBcXh7i4OMP7kJAQpKenY82aNfeFpzVr1qCoqAhvvvkmvvvuO6OO7+LiAl9f32atmYiIiGyDVV62q49Go4Gnp2edbefPn8eyZcuwadMmSKXGd3Xr1q1QKBTo1q0bFi5ciOLi4uYul4iIiKyUVY483SsrKwurV6/Ge++9Z9im0+kwbdo0rFixAoGBgcjOzjbqWE8++SSCg4Ph6+uLlJQULFq0CGfOnMG+ffsa3Een00Gn0xnea7XapneGiIiILJpFjTwtWbLkvsna976SkpLq7KNWqxEXF4dJkyZh1qxZhu2LFi1CeHg4pk+fLqqGhIQEDBs2DBEREZg6dSq++uor7N+/HydPnmxwn+XLl0MulxteSqVSXMeJiIjIakgEQRDMXUStgoICFBQUNNomKCgITk5OAGqCU0xMDPr374+NGzfWuTTXq1cvnDt3DhKJBAAgCAL0ej3s7Ozw2muvYenSpUbVJAgCZDIZNm/ejClTptTbpr6RJ6VSCY1GA3d3d6O+h4iIiMxLq9VCLpc/8Oe3RV22UygUUCgURrVVqVSIiYlBVFQUNmzYcN+cpu3bt6O8vNzw/sSJE5g5cyZ+/vlnhIaGGl1TamoqKisr4efn12AbmUwGmUxm9DGJiIjIellUeDKWWq1GdHQ0AgMDkZiYiPz8fMNntXfJ3RuQake0wsPD4eHhAaAmgMXGxmLTpk3o168fsrKysHXrVowePRoKhQLnz5/HggULEBkZicGDB7dM54iIiMiiWWV42rt3LzIzM5GZmYmAgIA6n4m5CllZWYn09HSUlZUBABwdHXHgwAGsWrUKJSUlUCqVGDNmDBYvXgw7O7tm7QMRERFZJ4ua82QrjL1mSkRERJbD2J/fFnW3HREREZGlY3giIiIiEoHhiYiIiEgEhiciIiIiERieiIiIiERgeCIiIiISgeGJiIiISASGJyIiIiIRGJ6IiIiIRGB4IiIiIhKB4YmIiIhIBIYnIiIiIhEYnoiIiIhEYHgiIiIiEoHhiYiIiEgEhiciIiIiERieiIiIiERgeCIiIiISgeGJiIiISASGJyIiIiIRGJ6IiIiIRGB4IiIiIhKB4YmIiIhIBIYnIiIiIhEYnoiIiIhEYHgiIiIiEoHhiYiIiEgEhiciIiIiERieiIiIiERgeCIiIiISgeGJiIiISASGJyIiIiIRGJ6IiIiIRGB4IiIiIhKB4YmIiIhIBIYnIiIiIhGsMjxdvnwZ8fHxCA4OhrOzM0JDQ7F48WJUVFTUaSeRSO57rV27ttFj63Q6vPjii1AoFHB1dcW4ceNw7do1U3aHiIiIrIi9uQtoirS0NOj1eqxbtw4dO3ZESkoKEhISUFpaisTExDptN2zYgLi4OMN7uVze6LHnzp2L3bt3Y9u2bfDy8sKCBQswduxYJCcnw87OziT9ISIiIushEQRBMHcRzWHFihVYs2YNsrOzDdskEgl27tyJCRMmGHUMjUYDb29vbN68GVOmTAEAqNVqKJVK7NmzByNHjjTqOFqtFnK5HBqNBu7u7qL7QkRERC3P2J/fVjnyVB+NRgNPT8/7ts+ePRuzZs1CcHAw4uPj8Ze//AVSaf1XK5OTk1FZWYkRI0YYtvn7+yMiIgLHjh1rMDzpdDrodLo6tQA1J4GIiIisQ+3P7QeNK9lEeMrKysLq1avx3nvv1dn+1ltvITY2Fs7Ozjhw4AAWLFiAgoICvP766/UeJzc3F46Ojmjbtm2d7T4+PsjNzW3w+5cvX46lS5fet12pVDahN0RERGROxcXFjU7zsajLdkuWLKk3hNztxIkT6NOnj+G9Wq3G0KFDMXToUHzyySeN7vvee+9h2bJlhpGhe33xxRd45pln6owiAcDw4cMRGhra4GTze0ee9Ho9bt26BS8vL0gkkkZrMoZWq4VSqUROTo7NXga09T7aev8A9tEW2Hr/APbRFpiyf4IgoLi4GP7+/g1epQIsbORp9uzZmDp1aqNtgoKCDH9Wq9WIiYnBwIEDsX79+gcef8CAAdBqtbhx4wZ8fHzu+9zX1xcVFRUoLCysM/qUl5eHQYMGNXhcmUwGmUxWZ5uHh8cD6xHL3d3dJv9DuJut99HW+wewj7bA1vsHsI+2wFT9e9CNZYCFhSeFQgGFQmFUW5VKhZiYGERFRWHDhg2NJsRap06dgpOTU4PBJioqCg4ODti3bx8mT54MALh+/TpSUlLw7rvvGt0PIiIisl0WFZ6MpVarER0djcDAQCQmJiI/P9/wma+vLwBg9+7dyM3NxcCBA+Hs7Iwff/wRr732Gv7yl78YRolUKhViY2OxadMm9OvXD3K5HPHx8ViwYAG8vLzg6emJhQsXonv37hg2bJhZ+kpERESWxSrD0969e5GZmYnMzEwEBATU+ax2CpeDgwM++ugjzJ8/H3q9HiEhIVi2bBleeOEFQ9vKykqkp6ejrKzMsG3lypWwt7fH5MmTUV5ejtjYWGzcuNGsazzJZDIsXrz4vkuDtsTW+2jr/QPYR1tg6/0D2EdbYAn9s6gJ40RERESWziofz0JERERkLgxPRERERCIwPBERERGJwPBEREREJALDk4VZvnw5JBIJ5s6d22i7w4cPIyoqCk5OTggJCWlw9XNLY0z/Dh06BIlEct8rLS2t5QoVYcmSJffVWrtkRkOs7fyJ7aO1ncNaKpUK06dPh5eXF1xcXNCrVy8kJyc3uo81nUux/bO28xgUFFRvvXffZX0vazp/gPg+Wts5rKqqwuuvv47g4GA4Ozsb7pTX6/WN7tfS59EqlyqwVSdOnMD69evRo0ePRttdunQJo0ePRkJCArZs2YKjR4/i+eefh7e3N5544okWqlY8Y/tXKz09vc7qsd7e3qYq7aF169YN+/fvN7xvbGkLaz1/YvpYy5rOYWFhIQYPHoyYmBh89913aNeuHbKyshp9WoA1ncum9K+WtZzHEydOoLq62vA+JSUFw4cPx6RJk+ptb03nr5bYPtaylnP4zjvvYO3atfj888/RrVs3JCUl4ZlnnoFcLsecOXPq3ccs51Egi1BcXCx06tRJ2LdvnzB06FBhzpw5DbZ96aWXhLCwsDrb/vrXvwoDBgwwcZVNJ6Z/P/74owBAKCwsbLH6HsbixYuFnj17Gt3eGs+f2D5a2zkUBEF4+eWXhUceeUTUPtZ0LpvSP2s8j3ebM2eOEBoaKuj1+no/t6bz15AH9dHazuGYMWOEmTNn1tk2ceJEYfr06Q3uY47zyMt2FuKFF17AmDFjjFrJ/JdffsGIESPqbBs5ciSSkpJQWVlpqhIfipj+1YqMjISfnx9iY2Px448/mrC6h5eRkQF/f38EBwdj6tSpyM7ObrCtNZ4/QFwfa1nTOfz666/Rp08fTJo0Ce3atUNkZCQ+/vjjRvexpnPZlP7VsqbzWKuiogJbtmzBzJkzG3xAuzWdv/oY08da1nIOH3nkERw4cAAXL14EAJw5cwZHjhzB6NGjG9zHHOeR4ckCbNu2DSdPnsTy5cuNap+bm3vfg419fHxQVVWFgoICU5T4UMT2z8/PD+vXr8f27duxY8cOdOnSBbGxsfjpp59MXGnT9O/fH5s2bcIPP/yAjz/+GLm5uRg0aBBu3rxZb3trO3+A+D5a2zkEgOzsbKxZswadOnXCDz/8gGeffRZ/+9vfsGnTpgb3saZz2ZT+WeN5rLVr1y4UFRXh6aefbrCNNZ2/+hjTR2s7hy+//DKmTZuGsLAwODg4IDIyEnPnzsW0adMa3Mcc55FznswsJycHc+bMwd69e+Hk5GT0fvf+liHcWSj+Qb99tLSm9K9Lly7o0qWL4f3AgQORk5ODxMREPProo6YqtclGjRpl+HP37t0xcOBAhIaG4vPPP8f8+fPr3cdazl8tsX20tnMIAHq9Hn369ME///lPADW/qaempmLNmjX485//3OB+1nIum9I/azyPtT799FOMGjUK/v7+jbazlvNXH2P6aG3n8D//+Q+2bNmCL774At26dcPp06cxd+5c+Pv7Y8aMGQ3u19LnkSNPZpacnIy8vDxERUXB3t4e9vb2OHz4MD744APY29vXmRhYy9fXF7m5uXW25eXlwd7eHl5eXi1VulGa0r/6DBgwABkZGSautnm4urqie/fuDdZrTeevIQ/qY30s/Rz6+fmha9eudbaFh4fj6tWrDe5jTeeyKf2rj6WfRwC4cuUK9u/fj1mzZjXazprO372M7WN9LPkc/v3vf8crr7yCqVOnonv37njqqacwb968Rq9cmOM8cuTJzGJjY3Hu3Lk625555hmEhYXh5ZdfrveOpoEDB2L37t11tu3duxd9+vSBg4ODSesVqyn9q8+pU6fg5+dnihKbnU6nw4ULFzBkyJB6P7em89eQB/WxPpZ+DgcPHoz09PQ62y5evIgOHTo0uI81ncum9K8+ln4eAWDDhg1o164dxowZ02g7azp/9zK2j/Wx5HNYVlYGqbTuuI6dnV2jSxWY5TyabCo6Ndm9d6O98sorwlNPPWV4n52dLbi4uAjz5s0Tzp8/L3z66aeCg4OD8NVXX5mhWvEe1L+VK1cKO3fuFC5evCikpKQIr7zyigBA2L59uxmqfbAFCxYIhw4dErKzs4Vff/1VGDt2rODm5iZcvnxZEATbOH9i+2ht51AQBOH48eOCvb298PbbbwsZGRnC1q1bBRcXF2HLli2GNtZ8LpvSP2s8j9XV1UJgYKDw8ssv3/eZNZ+/u4npo7WdwxkzZgjt27cXvvnmG+HSpUvCjh07BIVCIbz00kuGNpZwHhmeLNC94WLGjBnC0KFD67Q5dOiQEBkZKTg6OgpBQUHCmjVrWrbIh/Cg/r3zzjtCaGio4OTkJLRt21Z45JFHhG+//bblCzXSlClTBD8/P8HBwUHw9/cXJk6cKKSmpho+t4XzJ7aP1nYOa+3evVuIiIgQZDKZEBYWJqxfv77O59Z+LsX2zxrP4w8//CAAENLT0+/7zNrPXy0xfbS2c6jVaoU5c+YIgYGBgpOTkxASEiK89tprgk6nM7SxhPMoEYQ7s6qIiIiI6IE4YZyIiIhIBIYnIiIiIhEYnoiIiIhEYHgiIiIiEoHhiYiIiEgEhiciIiIiERieiIiIiERgeCKiVic6Ohpz5841dxlEZKUYnoiImtnGjRshkUgafR06dMjcZRJRE/HBwEREzWzKlCmIi4szvJ84cSIiIiKwbNkywzZPT09zlEZEzYAjT0TU6hUWFuLPf/4z2rZtCxcXF4waNQoZGRl12nz88cdQKpVwcXHB448/jvfffx8eHh71Hs/Z2Rm+vr6Gl6OjI1xcXO7bRkTWieGJiFq9p59+GklJSfj666/xyy+/QBAEjB49GpWVlQCAo0eP4tlnn8WcOXNw+vRpDB8+HG+//baZqyYic+FlOyJq1TIyMvD111/j6NGjGDRoEABg69atUCqV2LVrFyZNmoTVq1dj1KhRWLhwIQCgc+fOOHbsGL755htzlk5EZsKRJyJq1S5cuAB7e3v079/fsM3LywtdunTBhQsXAADp6eno169fnf3ufU9ErQfDExG1aoIgNLhdIpHc9+cH7UdEto/hiYhata5du6Kqqgq//fabYdvNmzdx8eJFhIeHAwDCwsJw/PjxOvslJSW1aJ1EZDkYnoioVevUqRPGjx+PhIQEHDlyBGfOnMH06dPRvn17jB8/HgDw4osvYs+ePXj//feRkZGBdevW4bvvvrtvNIqIWgeGJyJq9TZs2ICoqCiMHTsWAwcOhCAI2LNnDxwcHAAAgwcPxtq1a/H++++jZ8+e+P777zFv3jw4OTmZuXIiMgeJwAv3RESiJSQkIC0tDT///LO5SyGiFsalCoiIjJCYmIjhw4fD1dUV3333HT7//HN89NFH5i6LiMyAI09EREaYPHkyDh06hOLiYoSEhODFF1/Es88+a+6yiMgMGJ6IiIiIROCEcSIiIiIRGJ6IiIiIRGB4IiIiIhKB4YmIiIhIBIYnIiIiIhEYnoiIiIhEYHgiIiIiEoHhiYiIiEgEhiciIiIiEf4/R8VOIcUZBIYAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# File\n", + "table = open('gnat-sternberg.cooling','r')\n", + "lines = table.readlines()\n", + "\n", + "# Lists\n", + "logT = []\n", + "logL = []\n", + "\n", + "for line in lines[8:]:\n", + " \n", + " split = line.split()\n", + " \n", + " logT.append(float(split[0]))\n", + " logL.append(float(split[1]))\n", + "\n", + "# Close\n", + "table.close()\n", + "\n", + "# Plot\n", + "fig,ax = plt.subplots()\n", + "ax.plot(logT,logL)\n", + "ax.set_xlabel('log T')\n", + "ax.set_ylabel('log $\\Lambda$')\n", + "ax.set_ylim(-25,-21)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6c23cebb-cb61-41fc-91fd-1fe35faad040", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "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.11.3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/src/eos/adiabatic_glmmhd.cpp b/src/eos/adiabatic_glmmhd.cpp index 24b8a155..124fd2a7 100644 --- a/src/eos/adiabatic_glmmhd.cpp +++ b/src/eos/adiabatic_glmmhd.cpp @@ -41,6 +41,7 @@ void AdiabaticGLMMHDEOS::ConservedToPrimitive(MeshData *md) const { const auto nhydro = pkg->Param("nhydro"); const auto nscalars = pkg->Param("nscalars"); + auto this_on_device = (*this); parthenon::par_for( @@ -50,7 +51,23 @@ void AdiabaticGLMMHDEOS::ConservedToPrimitive(MeshData *md) const { const auto &cons = cons_pack(b); auto &prim = prim_pack(b); // auto &nu = entropy_pack(b); + + // Adding the gks gjs gis + + // Getting the global indexing + + auto pmb = md->GetBlockData(b)->GetBlockPointer(); + auto pm = pmb->pmy_mesh; + auto hydro_pkg = pmb->packages.Get("Hydro"); + + const auto gis = pmb->loc.lx1 * pmb->block_size.nx1; + const auto gjs = pmb->loc.lx2 * pmb->block_size.nx2; + const auto gks = pmb->loc.lx3 * pmb->block_size.nx3; + + // ... + + // - return this_on_device.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i); + return this_on_device.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i, gks, gjs,gis); }); } diff --git a/src/eos/adiabatic_glmmhd.hpp b/src/eos/adiabatic_glmmhd.hpp index baebc681..9b59f3e2 100644 --- a/src/eos/adiabatic_glmmhd.hpp +++ b/src/eos/adiabatic_glmmhd.hpp @@ -30,10 +30,10 @@ class AdiabaticGLMMHDEOS : public EquationOfState { gamma_{gamma} {} void ConservedToPrimitive(MeshData *md) const override; - + KOKKOS_INLINE_FUNCTION Real GetGamma() const { return gamma_; } - + //---------------------------------------------------------------------------------------- // \!fn Real EquationOfState::SoundSpeed(Real prim[NHYDRO]) // \brief returns adiabatic sound speed given vector of primitive variables @@ -53,7 +53,7 @@ class AdiabaticGLMMHDEOS : public EquationOfState { return std::sqrt(0.5 * (qsq + std::sqrt(tmp * tmp + 4.0 * asq * ct2)) / d); } // - + //---------------------------------------------------------------------------------------- // \!fn Real EquationOfState::ConsToPrim(View4D cons, View4D prim, const int& k, const // int& j, const int& i) \brief Fills an array of primitives given an array of @@ -61,7 +61,7 @@ class AdiabaticGLMMHDEOS : public EquationOfState { template KOKKOS_INLINE_FUNCTION void ConsToPrim(View4D cons, View4D prim, const int &nhydro, const int &nscalars, const int &k, const int &j, - const int &i) const { + const int &i, const int &gks, const int &gjs, const int &gis) const { auto gam = GetGamma(); auto gm1 = gam - 1.0; auto density_floor_ = GetDensityFloor(); diff --git a/src/eos/adiabatic_hydro.cpp b/src/eos/adiabatic_hydro.cpp index 20e53a75..475a754c 100644 --- a/src/eos/adiabatic_hydro.cpp +++ b/src/eos/adiabatic_hydro.cpp @@ -30,26 +30,40 @@ using parthenon::ParArray4D; // Container &rc, // int il, int iu, int jl, int ju, int kl, int ku) // \brief Converts conserved into primitive variables in adiabatic hydro. + void AdiabaticHydroEOS::ConservedToPrimitive(MeshData *md) const { auto const cons_pack = md->PackVariables(std::vector{"cons"}); auto prim_pack = md->PackVariables(std::vector{"prim"}); auto ib = md->GetBlockData(0)->GetBoundsI(IndexDomain::entire); auto jb = md->GetBlockData(0)->GetBoundsJ(IndexDomain::entire); auto kb = md->GetBlockData(0)->GetBoundsK(IndexDomain::entire); - + auto pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro"); const auto nhydro = pkg->Param("nhydro"); const auto nscalars = pkg->Param("nscalars"); - + auto this_on_device = (*this); - + parthenon::par_for( DEFAULT_LOOP_PATTERN, "ConservedToPrimitive", parthenon::DevExecSpace(), 0, cons_pack.GetDim(5) - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) { + + // Getting the global indexing + + auto pmb = md->GetBlockData(b)->GetBlockPointer(); + auto pm = pmb->pmy_mesh; + auto hydro_pkg = pmb->packages.Get("Hydro"); + + const auto gis = pmb->loc.lx1 * pmb->block_size.nx1; + const auto gjs = pmb->loc.lx2 * pmb->block_size.nx2; + const auto gks = pmb->loc.lx3 * pmb->block_size.nx3; + + // ... + const auto &cons = cons_pack(b); auto &prim = prim_pack(b); - - return this_on_device.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i); + + return this_on_device.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i, gks, gjs ,gis); }); } diff --git a/src/eos/adiabatic_hydro.hpp b/src/eos/adiabatic_hydro.hpp index bf317923..0addb9ee 100644 --- a/src/eos/adiabatic_hydro.hpp +++ b/src/eos/adiabatic_hydro.hpp @@ -51,15 +51,16 @@ class AdiabaticHydroEOS : public EquationOfState { template KOKKOS_INLINE_FUNCTION void ConsToPrim(View4D cons, View4D prim, const int &nhydro, const int &nscalars, const int &k, const int &j, - const int &i) const { + const int &i, const int &gks, const int &gjs, const int &gis) const { + Real gm1 = GetGamma() - 1.0; auto density_floor_ = GetDensityFloor(); auto pressure_floor_ = GetPressureFloor(); auto e_floor_ = GetInternalEFloor(); - + auto velocity_ceiling_ = GetVelocityCeiling(); auto e_ceiling_ = GetInternalECeiling(); - + Real &u_d = cons(IDN, k, j, i); Real &u_m1 = cons(IM1, k, j, i); Real &u_m2 = cons(IM2, k, j, i); @@ -71,7 +72,11 @@ class AdiabaticHydroEOS : public EquationOfState { Real &w_vy = prim(IV2, k, j, i); Real &w_vz = prim(IV3, k, j, i); Real &w_p = prim(IPR, k, j, i); - + + if (u_d <= 0.0) {std::cout << "(k,j,i)=(" << k << "," << j << "," << i << ") \n" ; + std::cout << "(gks,gjs,gis)=(" << gks << "," << gjs << "," << gis << ") \n" ; + } + // Let's apply floors explicitly, i.e., by default floor will be disabled (<=0) // and the code will fail if a negative density is encountered. PARTHENON_REQUIRE(u_d > 0.0 || density_floor_ > 0.0, @@ -88,7 +93,7 @@ class AdiabaticHydroEOS : public EquationOfState { Real e_k = 0.5 * di * (SQR(u_m1) + SQR(u_m2) + SQR(u_m3)); w_p = gm1 * (u_e - e_k); - + // apply velocity ceiling. By default ceiling is std::numeric_limits::infinity() const Real w_v2 = SQR(w_vx) + SQR(w_vy) + SQR(w_vz); if (w_v2 > SQR(velocity_ceiling_)) { @@ -105,11 +110,29 @@ class AdiabaticHydroEOS : public EquationOfState { u_e -= e_k - e_k_new; e_k = e_k_new; } - + // Let's apply floors explicitly, i.e., by default floor will be disabled (<=0) // and the code will fail if a negative pressure is encountered. + + // Trying to understand why negative pressure appear in tests + + if (w_p <= 0.0) {std::cout << "w_p = gm1 * (u_e - e_k)"; + std::cout << "w_p=" << w_p << "\n" ; + std::cout << "gm1=" << gm1 << "\n" ; + std::cout << "u_e=" << u_e << "\n" ; + std::cout << "e_k=" << e_k << "\n" ; + std::cout << "(k,j,i)=(" << k << "," << j << "," << i << ") \n" ; + std::cout << "(gks,gjs,gis)=(" << gks << "," << gjs << "," << gis << ") \n" ; + } + + //if (pressure_floor_ < 0.0) {std::cout << "pressure_floor_" << pressure_floor_ << "\n" ;} + //if (e_floor_ < 0.0) {std::cout << "e_floor_" << e_floor_ << "\n" ;} + + // The first argument check whether one of the conditions is true + // By default, floors are deactivated, ie. pressure floor and e_floor + // are -1. The code will eventually crash when w_p turns < 0. PARTHENON_REQUIRE(w_p > 0.0 || pressure_floor_ > 0.0 || e_floor_ > 0.0, - "Got negative pressure. Consider enabling first-order flux " + "Wow. Got negative pressure. Consider enabling first-order flux " "correction or setting a reasonble pressure or temperature floor."); // Pressure floor (if present) takes precedence over temperature floor diff --git a/src/hydro/srcterms/old/tabular_cooling_edit.cpp b/src/hydro/srcterms/old/tabular_cooling_edit.cpp new file mode 100644 index 00000000..b8d7a0fd --- /dev/null +++ b/src/hydro/srcterms/old/tabular_cooling_edit.cpp @@ -0,0 +1,743 @@ +//======================================================================================== +// AthenaPK - a performance portable block structured AMR astrophysical MHD code. +// Copyright (c) 2021-2023, Athena-Parthenon Collaboration. All rights reserved. +// Licensed under the 3-clause BSD License, see LICENSE file for details +//======================================================================================== +//! \file tabular_cooling.cpp +// \brief Applies tabular cooling +// +//======================================================================================== + +// C++ headers +#include +#include +#include + +// Parthenon headers +#include +#include +#include +#include + +// AthenaPK headers +#include "../../units.hpp" +#include "tabular_cooling.hpp" +#include "utils/error_checking.hpp" + +namespace cooling { +using namespace parthenon; + +TabularCooling::TabularCooling(ParameterInput *pin, + std::shared_ptr hydro_pkg) { + auto units = hydro_pkg->Param("units"); + + const std::string table_filename = pin->GetString("cooling", "table_filename"); + + const int log_temp_col = pin->GetOrAddInteger("cooling", "log_temp_col", 0); + const int log_lambda_col = pin->GetOrAddInteger("cooling", "log_lambda_col", 1); + + // Heating function as defined in Wladimir et al. 2021, ie. insure thermal eq. at n=0.1cm-3 + // and T = 1e4. By default, set to zero (heating deactivated) + // Provided in cgs units, ie. erg / s, then rescaled to code units + const Real heating_const_cgs = pin->GetOrAddReal("cooling", "heating_const", 0.0); + const Real gamma_units = heating_const_cgs * (units.erg() / units.s()); + + // Just for verification, check what units.s() is + //std::cout << "units.s()=" << units.s(); + + // Convert erg cm^3/s to code units + const Real lambda_units_cgs = pin->GetReal("cooling", "lambda_units_cgs"); + const Real lambda_units = + lambda_units_cgs / (units.erg() * pow(units.cm(), 3) / units.s()); + + // Heating parameters + const Real + + std::cout << "Entering lambda cout loop\n" + + const auto integrator_str = pin->GetOrAddString("cooling", "integrator", "rk12"); + if (integrator_str == "rk12") { + //std::cout << "Cooling integrator is rk12"; + integrator_ = CoolIntegrator::rk12; + } else if (integrator_str == "rk45") + //std::cout << "Cooling integrator is rk45"; + integrator_ = CoolIntegrator::rk45; + } else if (integrator_str == "townsend") { + //std::cout << "Cooling integrator is Townsend"; + integrator_ = CoolIntegrator::townsend; + } else { + integrator_ = CoolIntegrator::undefined; + } + + max_iter_ = pin->GetOrAddInteger("cooling", "max_iter", 100); + cooling_time_cfl_ = pin->GetOrAddReal("cooling", "cfl", 0.1); + d_log_temp_tol_ = pin->GetOrAddReal("cooling", "d_log_temp_tol", 1e-8); + d_e_tol_ = pin->GetOrAddReal("cooling", "d_e_tol", 1e-8); + // negative means disabled + T_floor_ = pin->GetOrAddReal("hydro", "Tfloor", -1.0); + + std::stringstream msg; + + /**************************************** + * Read tab file with IOWrapper + ****************************************/ + IOWrapper input; + input.Open(table_filename.c_str(), IOWrapper::FileMode::read); + + /**************************************** + * Read tab file from IOWrapper into a stringstream tab + ****************************************/ + std::stringstream tab_ss; + const int bufsize = 4096; + char *buf = new char[bufsize]; + std::ptrdiff_t ret; + parthenon::IOWrapperSizeT word_size = sizeof(char); + + do { + if (Globals::my_rank == 0) { // only the master process reads the cooling table + ret = input.Read(buf, word_size, bufsize); + } +#ifdef MPI_PARALLEL + // then broadcasts it + // no need for fence as cooling table is independent of execution/memory space + MPI_Bcast(&ret, sizeof(std::ptrdiff_t), MPI_BYTE, 0, MPI_COMM_WORLD); + MPI_Bcast(buf, ret, MPI_BYTE, 0, MPI_COMM_WORLD); +#endif + tab_ss.write(buf, ret); // add the buffer into the stream + } while (ret == bufsize); // till EOF (or par_end is found) + + delete[] buf; + input.Close(); + + /**************************************** + * Determine log_temps and and log_lambdas vectors + ****************************************/ + std::vector log_temps, log_lambdas; + std::string line; + std::size_t first_char; + while (tab_ss.good()) { + getline(tab_ss, line); + if (line.empty()) continue; // skip blank line + first_char = line.find_first_not_of(" "); // skip white space + if (first_char == std::string::npos) continue; // line is all white space + if (line.compare(first_char, 1, "#") == 0) continue; // skip comments + + // Parse the numbers on the line + std::istringstream iss(line); + std::vector line_data{std::istream_iterator{iss}, + std::istream_iterator{}}; + // Check size + if (line_data.empty() || line_data.size() <= std::max(log_temp_col, log_lambda_col)) { + msg << "### FATAL ERROR in function [TabularCooling::TabularCooling]" << std::endl + << "Index " << std::max(log_temp_col, log_lambda_col) << " out of range on \"" + << line << "\"" << std::endl; + PARTHENON_FAIL(msg); + } + + try { + const Real log_temp = std::stod(line_data[log_temp_col]); + const Real log_lambda = std::stod(line_data[log_lambda_col]); + + // Add to growing list + log_temps.push_back(log_temp); + log_lambdas.push_back(log_lambda - std::log10(lambda_units)); + + } catch (const std::invalid_argument &ia) { + msg << "### FATAL ERROR in function [TabularCooling::TabularCooling]" << std::endl + << "Number: \"" << ia.what() << "\" could not be parsed as double" << std::endl; + PARTHENON_FAIL(msg); + } + } + + /**************************************** + * Check some assumtions about the cooling table + ****************************************/ + + // Reconstructing the cooling function in new coordinate system + std::cout << "Entering lambda cout loop\n" + for (int i = 0; i < log_temps.size(); i++) { + std::cout << "log_temps" << log_temps[i] << ", log_lambdas" << log_lambdas[i] << "\n"; + } + + // Ensure at least two data points in the table to interpolate from + if (log_temps.size() < 2 || log_lambdas.size() < 2) { + msg << "### FATAL ERROR in function [TabularCooling::TabularCooling]" << std::endl + << "Not enough data to interpolate cooling" << std::endl; + PARTHENON_FAIL(msg); + } + + // Ensure that the first log_temp is increasing + const Real log_temp_start = log_temps[0]; + const Real d_log_temp = log_temps[1] - log_temp_start; + + if (d_log_temp <= 0) { + msg << "### FATAL ERROR in function [TabularCooling::TabularCooling]" << std::endl + << "second log_temp in table is descreasing" << std::endl; + PARTHENON_FAIL(msg); + } + + // Ensure that log_temps is evenly spaced + for (size_t i = 1; i < log_temps.size(); i++) { + const Real d_log_temp_i = log_temps[i] - log_temps[i - 1]; + + if (d_log_temp_i < 0) { + msg << "### FATAL ERROR in function [TabularCooling::TabularCooling]" << std::endl + << "log_temp in table is descreasing at i= " << i + << " log_temp= " << log_temps[i] << std::endl; + PARTHENON_FAIL(msg); + } + + // The Dedt() function currently relies on an equally spaced cooling table for faster + // lookup (direct indexing). It is used in the subcycling method and in order to + // restrict `dt` by a cooling cfl. + if (((integrator_ != CoolIntegrator::townsend) || (cooling_time_cfl_ > 0.0)) && + (fabs(d_log_temp_i - d_log_temp) / d_log_temp > d_log_temp_tol_)) { + msg << "### FATAL ERROR in function [TabularCooling::TabularCooling]" << std::endl + << "d_log_temp in table is uneven at i=" << i << " log_temp=" << log_temps[i] + << " d_log_temp= " << d_log_temp << " d_log_temp_i= " << d_log_temp_i + << " diff= " << d_log_temp_i - d_log_temp + << " rel_diff= " << fabs(d_log_temp_i - d_log_temp) / d_log_temp + << " tol= " << d_log_temp_tol_ << std::endl; + PARTHENON_FAIL(msg); + } + } + + /**************************************** + * Move values read into the data table + ****************************************/ + + n_temp_ = log_temps.size(); + log_temp_start_ = log_temps[0]; + log_temp_final_ = log_temps[n_temp_ - 1]; + d_log_temp_ = d_log_temp; + lambda_final_ = std::pow(10.0, log_lambdas[n_temp_ - 1]); + + // Setup log_lambdas_ used in Dedt() + { + // log_lambdas is used if the integrator isn't Townsend, if the cooling CFL + // is set, or if cooling time is a extra derived field. Since we don't have + // a good way to check the last condition we always initialize log_lambdas_ + log_lambdas_ = ParArray1D("log_lambdas_", n_temp_); + + // Read log_lambdas in host_log_lambdas, changing to code units along the way + auto host_log_lambdas = Kokkos::create_mirror_view(log_lambdas_); + for (unsigned int i = 0; i < n_temp_; i++) { + host_log_lambdas(i) = log_lambdas[i]; + } + // Copy host_log_lambdas into device memory + Kokkos::deep_copy(log_lambdas_, host_log_lambdas); + } + // Setup Townsend cooling, i.e., precalulcate piecewise powerlaw approx. + if (integrator_ == CoolIntegrator::townsend) { + lambdas_ = ParArray1D("lambdas_", n_temp_); + temps_ = ParArray1D("temps_", n_temp_); + + // Read log_lambdas in host_lambdas, changing to code units along the way + auto host_lambdas = Kokkos::create_mirror_view(lambdas_); + auto host_temps = Kokkos::create_mirror_view(temps_); + for (unsigned int i = 0; i < n_temp_; i++) { + host_lambdas(i) = std::pow(10.0, log_lambdas[i]); + host_temps(i) = std::pow(10.0, log_temps[i]); + } + // Copy host_lambdas into device memory + Kokkos::deep_copy(lambdas_, host_lambdas); + Kokkos::deep_copy(temps_, host_temps); + + // Coeffs are for intervals, i.e., only n_temp_ - 1 entries + const auto n_bins = n_temp_ - 1; + townsend_Y_k_ = ParArray1D("townsend_Y_k_", n_bins); + townsend_alpha_k_ = ParArray1D("townsend_alpha_k_", n_bins); + + // Initialize on host (make *this* recursion simpler) + auto host_townsend_Y_k = Kokkos::create_mirror_view(townsend_Y_k_); + auto host_townsend_alpha_k = Kokkos::create_mirror_view(townsend_alpha_k_); + + // Initialize piecewise power law indices + for (unsigned int i = 0; i < n_bins; i++) { + // Could use log_lambdas_ here, but using lambdas_ instead as they've already been + // converted to code units. + host_townsend_alpha_k(i) = + (std::log10(host_lambdas(i + 1)) - std::log10(host_lambdas(i))) / + (log_temps[i + 1] - log_temps[i]); + PARTHENON_REQUIRE(host_townsend_alpha_k(i) != 1.0, + "Need to implement special case for Townsend piecewise fits."); + } + + // Calculate TEF (temporal evolution functions Y_k recursively), (Eq. A6) + host_townsend_Y_k(n_bins - 1) = 0.0; // Last Y_N = Y(T_ref) = 0 + + for (int i = n_bins - 2; i >= 0; i--) { + const auto alpha_k_m1 = host_townsend_alpha_k(i) - 1.0; + const auto step = (host_lambdas(n_bins) / host_lambdas(i)) * + (host_temps(i) / host_temps(n_bins)) * + (std::pow(host_temps(i) / host_temps(i + 1), alpha_k_m1) - 1.0) / + alpha_k_m1; + + host_townsend_Y_k(i) = host_townsend_Y_k(i + 1) - step; + } + + Kokkos::deep_copy(townsend_alpha_k_, host_townsend_alpha_k); + Kokkos::deep_copy(townsend_Y_k_, host_townsend_Y_k); + } + + // Create a lightweight object for computing cooling rates within kernels + const auto mbar_over_kb = hydro_pkg->Param("mbar_over_kb"); + const auto adiabatic_index = hydro_pkg->Param("AdiabaticIndex"); + const auto He_mass_fraction = hydro_pkg->Param("He_mass_fraction"); + + cooling_table_obj_ = CoolingTableObj(log_lambdas_, log_temp_start_, log_temp_final_, + d_log_temp_, n_temp_, mbar_over_kb, + adiabatic_index, 1.0 - He_mass_fraction, units); +} + +void TabularCooling::SrcTerm(MeshData *md, const Real dt) const { + if (integrator_ == CoolIntegrator::rk12) { + SubcyclingFixedIntSrcTerm(md, dt, RK12Stepper()); + } else if (integrator_ == CoolIntegrator::rk45) { + SubcyclingFixedIntSrcTerm(md, dt, RK45Stepper()); + } else if (integrator_ == CoolIntegrator::townsend) { + TownsendSrcTerm(md, dt); + } else { + PARTHENON_FAIL("Unknown cooling integrator."); + } +} + +template +void TabularCooling::SubcyclingFixedIntSrcTerm(MeshData *md, const Real dt_, + const RKStepper rk_stepper) const { + + const auto dt = dt_; // HACK capturing parameters still broken with Cuda 11.6 ... + auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro"); + const bool mhd_enabled = hydro_pkg->Param("fluid") == Fluid::glmmhd; + // Grab member variables for compiler + + const CoolingTableObj cooling_table_obj = cooling_table_obj_; + const auto gm1 = (hydro_pkg->Param("AdiabaticIndex") - 1.0); + const auto mbar_gm1_over_kb = hydro_pkg->Param("mbar_over_kb") * gm1; + + const unsigned int max_iter = max_iter_; + + const Real min_sub_dt = dt / max_iter; + + const Real d_e_tol = d_e_tol_; + + // Determine the cooling floor, whichever is higher of the cooling table floor + // or fluid solver floor + const auto temp_cool_floor = std::pow(10.0, log_temp_start_); // low end of cool table + const Real temp_floor = (T_floor_ > temp_cool_floor) ? T_floor_ : temp_cool_floor; + + const Real internal_e_floor = temp_floor / mbar_gm1_over_kb; // specific internal en. + + // Grab some necessary variables + const auto &prim_pack = md->PackVariables(std::vector{"prim"}); + const auto &cons_pack = md->PackVariables(std::vector{"cons"}); + // need to include ghost zones as this source is called prior to the other fluxes when + // split + IndexRange ib = md->GetBlockData(0)->GetBoundsI(IndexDomain::entire); + IndexRange jb = md->GetBlockData(0)->GetBoundsJ(IndexDomain::entire); + IndexRange kb = md->GetBlockData(0)->GetBoundsK(IndexDomain::entire); + + par_for( + DEFAULT_LOOP_PATTERN, "TabularCooling::SubcyclingSplitSrcTerm", DevExecSpace(), 0, + cons_pack.GetDim(5) - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i) { + auto &cons = cons_pack(b); + auto &prim = prim_pack(b); + // Need to use `cons` here as prim may still contain state at t_0; + const Real rho = cons(IDN, k, j, i); + // TODO(pgrete) with potentially more EOS, a separate get_pressure (or similar) + // function could be useful. + Real internal_e = + cons(IEN, k, j, i) - 0.5 * + (SQR(cons(IM1, k, j, i)) + SQR(cons(IM2, k, j, i)) + + SQR(cons(IM3, k, j, i))) / + rho; + if (mhd_enabled) { + internal_e -= 0.5 * (SQR(cons(IB1, k, j, i)) + SQR(cons(IB2, k, j, i)) + + SQR(cons(IB3, k, j, i))); + } + internal_e /= rho; + const Real internal_e_initial = internal_e; + + bool dedt_valid = true; + + // Wrap DeDt into a functor for the RKStepper + auto DeDt_wrapper = [&](const Real t, const Real e, bool &valid) { + return cooling_table_obj.DeDt(e, rho, valid); + }; + + Real sub_t = 0; // current subcycle time + // Try full dt. If error is too large adaptive timestepping will reduce sub_dt + Real sub_dt = dt; + + // Check if cooling is actually happening, e.g., when T below T_cool_min or if + // temperature is already below floor. + const Real dedt_initial = DeDt_wrapper(0.0, internal_e_initial, dedt_valid); + if (dedt_initial == 0.0 || internal_e_initial <= internal_e_floor) { + return; // Function end here if no cooling needed, ie. if initially edot = 0 + } + + // Use minumum subcycle timestep when d_e_tol == 0 + if (d_e_tol == 0) { + sub_dt = min_sub_dt; + } + + unsigned int sub_iter = 0; + // check for dedt != 0.0 required in case cooling floor it hit during subcycling + while ((sub_t * (1 + KEpsilon_) < dt) && + (DeDt_wrapper(sub_t, internal_e, dedt_valid) != 0.0)) { + + if (sub_iter > max_iter) { + // Due to sub_dt >= min_dt, this error should never happen + PARTHENON_FAIL( + "FATAL ERROR in [TabularCooling::SubcyclingFixedIntSrcTerm]: Sub " + "cycles exceed max_iter (This should be impossible)"); + } + + // Next higher order estimate + Real internal_e_next_h; + // Error in estimate of higher order + Real d_e_err; + // Number of attempts on this subcycle + unsigned int sub_attempt = 0; + // Whether to reattempt this subcycle + bool reattempt_sub = true; + do { + // Next lower order estimate + Real internal_e_next_l; + // Do one dual order RK step + dedt_valid = true; + RKStepper::Step(sub_t, sub_dt, internal_e, DeDt_wrapper, internal_e_next_h, + internal_e_next_l, dedt_valid); + + sub_attempt++; + + if (!dedt_valid) { + if (sub_dt == min_sub_dt) { + // Cooling is so fast that even the minimum subcycle dt would lead to + // negative internal energy -- so just cool to the floor of the cooling + // table + sub_dt = (dt - sub_t); + internal_e_next_h = internal_e_floor; + reattempt_sub = false; + } else { + reattempt_sub = true; + sub_dt = min_sub_dt; + } + } else { + + // Compute error + d_e_err = fabs((internal_e_next_h - internal_e_next_l) / internal_e_next_h); + + reattempt_sub = false; + // Accepting or reattempting the subcycle: + // + // -If the error is small, accept the subcycle + // + // -If the error on the subcycle is too high, compute a new time + // step to reattempt the subcycle + // -But if the new time step is smaller than the minimum subcycle + // time step (total step duration/ max iterations), just use the + // minimum subcycle time step instead + + if (std::isnan(d_e_err)) { + reattempt_sub = true; + sub_dt = min_sub_dt; + } else if (d_e_err >= d_e_tol && sub_dt > min_sub_dt) { + // Reattempt this subcycle + reattempt_sub = true; + // Error was too high, shrink the timestep + if (d_e_tol == 0) { + sub_dt = min_sub_dt; + } else { + sub_dt = RKStepper::OptimalStep(sub_dt, d_e_err, d_e_tol); + } + // Don't drop timestep under maximum iteration count + if (sub_dt < min_sub_dt || sub_attempt >= max_iter) { + sub_dt = min_sub_dt; + } + } + } + + } while (reattempt_sub); + + // Accept this subcycle + sub_t += sub_dt; + + internal_e = internal_e_next_h; // !!!! internal_e updated here !!!! + + // skip to the end of subcycling if error is 0 (very unlikely) + if (d_e_err == 0) { + sub_dt = dt - sub_t; + } else { + // Grow the timestep + // (or shrink in case d_e_err >= d_e_tol and sub_dt is already at min_sub_dt) + sub_dt = RKStepper::OptimalStep(sub_dt, d_e_err, d_e_tol); + } + + if (d_e_tol == 0) { + sub_dt = min_sub_dt; + } + + // Don't drop timestep under the minimum step size + sub_dt = std::max(sub_dt, min_sub_dt); + + // Limit by end time + sub_dt = std::min(sub_dt, dt - sub_t); + + sub_iter++; + } + + // If cooled below floor, reset to floor value. + // This could happen if the floor value is larger than the lower end of the + // cooling table or if they are close and the last subcycle in the cooling above + // the lower end pushed the temperature below the lower end (and the floor). + // Literally do: if internal_e > internal_e_floor, internal_e remain unchanged + // otherwise, internal_e get internal_e_floor. + // internal_e is already updated before (see above) + internal_e = (internal_e > internal_e_floor) ? internal_e : internal_e_floor; + + // Remove the cooling from the total energy density + cons(IEN, k, j, i) += rho * (internal_e - internal_e_initial); + // Latter technically not required if no other tasks follows before + // ConservedToPrim conversion, but keeping it for now (better safe than sorry). + prim(IPR, k, j, i) = rho * internal_e * gm1; + }); +} + +void TabularCooling::TownsendSrcTerm(parthenon::MeshData *md, + const parthenon::Real dt_) const { + auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro"); + const bool mhd_enabled = hydro_pkg->Param("fluid") == Fluid::glmmhd; + + // Grab member variables for compiler + const auto dt = dt_; // HACK capturing parameters still broken with Cuda 11.6 ... + + const auto units = hydro_pkg->Param("units"); + const auto gm1 = (hydro_pkg->Param("AdiabaticIndex") - 1.0); + const auto mbar_gm1_over_kb = hydro_pkg->Param("mbar_over_kb") * gm1; + const Real X_by_mh2 = + std::pow((1 - hydro_pkg->Param("He_mass_fraction")) / units.mh(), 2); + + const auto lambdas = lambdas_; + const auto temps = temps_; + const auto alpha_k = townsend_alpha_k_; + const auto Y_k = townsend_Y_k_; + + const auto internal_e_floor = T_floor_ / mbar_gm1_over_kb; + const auto temp_cool_floor = std::pow(10.0, log_temp_start_); // low end of cool table + + // Grab some necessary variables + const auto &prim_pack = md->PackVariables(std::vector{"prim"}); + const auto &cons_pack = md->PackVariables(std::vector{"cons"}); + // need to include ghost zones as this source is called prior to the other fluxes when + // split + IndexRange ib = md->GetBlockData(0)->GetBoundsI(IndexDomain::entire); + IndexRange jb = md->GetBlockData(0)->GetBoundsJ(IndexDomain::entire); + IndexRange kb = md->GetBlockData(0)->GetBoundsK(IndexDomain::entire); + + const auto nbins = alpha_k.extent_int(0); + + // Get reference values + const auto temp_final = std::pow(10.0, log_temp_final_); + const auto lambda_final = lambda_final_; + + par_for( + DEFAULT_LOOP_PATTERN, "TabularCooling::TownsendSrcTerm", DevExecSpace(), 0, + cons_pack.GetDim(5) - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i) { + auto &cons = cons_pack(b); + auto &prim = prim_pack(b); + // Need to use `cons` here as prim may still contain state at t_0; + const auto rho = cons(IDN, k, j, i); + // TODO(pgrete) with potentially more EOS, a separate get_pressure (or similar) + // function could be useful. + auto internal_e = + cons(IEN, k, j, i) - 0.5 * + (SQR(cons(IM1, k, j, i)) + SQR(cons(IM2, k, j, i)) + + SQR(cons(IM3, k, j, i))) / + rho; + if (mhd_enabled) { + internal_e -= 0.5 * (SQR(cons(IB1, k, j, i)) + SQR(cons(IB2, k, j, i)) + + SQR(cons(IB3, k, j, i))); + } + internal_e /= rho; + + // If temp is below floor, reset and return + if (internal_e <= internal_e_floor) { + // Remove the cooling from the total energy density + cons(IEN, k, j, i) += rho * (internal_e_floor - internal_e); + // Latter technically not required if no other tasks follows before + // ConservedToPrim conversion, but keeping it for now (better safe than sorry). + prim(IPR, k, j, i) = rho * internal_e_floor * gm1; + return; + } + + auto temp = mbar_gm1_over_kb * internal_e; + // Temperature is above floor (see conditional above) but below cooling table: + // -> no cooling + if (temp < temp_cool_floor) { + return; + } + const Real n_h2_by_rho = rho * X_by_mh2; + + // Get the index of the right temperature bin + // TODO(?) this could be optimized for using a binary search + auto idx = 0; + while ((idx < nbins - 1) && (temps(idx + 1) < temp)) { + idx += 1; + } + + // Compute the Temporal Evolution Function Y(T) (Eq. A5) + const auto alpha_k_m1 = alpha_k(idx) - 1.0; + const auto tef = + Y_k(idx) + (lambda_final / lambdas(idx)) * (temps(idx) / temp_final) * + (std::pow(temps(idx) / temp, alpha_k_m1) - 1.0) / alpha_k_m1; + + // Compute the adjusted TEF for new timestep (Eqn. 26) (term in brackets) + const auto tef_adj = + tef + lambda_final * dt / temp_final * mbar_gm1_over_kb * n_h2_by_rho; + + // TEF is a strictly decreasing function and new_tef > tef + // Check if the new TEF falls into a lower bin, i.e., find the right bin for A7 + // If so, update slopes and coefficients + while ((idx > 0) && (tef_adj > Y_k(idx))) { + idx -= 1; + } + + // Compute the Inverse Temporal Evolution Function Y^{-1}(Y) (Eq. A7) + const auto temp_new = + temps(idx) * + std::pow(1 - (1.0 - alpha_k(idx)) * (lambdas(idx) / lambda_final) * + (temp_final / temps(idx)) * (tef_adj - Y_k(idx)), + 1.0 / (1.0 - alpha_k(idx))); + // Set new temp (at the lowest to the lower end of the cooling table) + const auto internal_e_new = temp_new > temp_cool_floor + ? temp_new / mbar_gm1_over_kb + : temp_cool_floor / mbar_gm1_over_kb; + cons(IEN, k, j, i) += rho * (internal_e_new - internal_e); + // Latter technically not required if no other tasks follows before + // ConservedToPrim conversion, but keeping it for now (better safe than sorry). + prim(IPR, k, j, i) = rho * internal_e_new * gm1; + }); +} + +Real TabularCooling::EstimateTimeStep(MeshData *md) const { + if (cooling_time_cfl_ <= 0.0) { + return std::numeric_limits::max(); + } + + if (isnan(cooling_time_cfl_) || isinf(cooling_time_cfl_)) { + return std::numeric_limits::infinity(); + } + + auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro"); + + const CoolingTableObj cooling_table_obj = cooling_table_obj_; + const auto gm1 = (hydro_pkg->Param("AdiabaticIndex") - 1.0); + const auto mbar_gm1_over_kb = hydro_pkg->Param("mbar_over_kb") * gm1; + + // Determine the cooling floor, whichever is higher of the cooling table floor + // or fluid solver floor + const auto temp_cool_floor = std::pow(10.0, log_temp_start_); // low end of cool table + const Real temp_floor = (T_floor_ > temp_cool_floor) ? T_floor_ : temp_cool_floor; + + const Real internal_e_floor = temp_floor / mbar_gm1_over_kb; // specific internal en. + + // Grab some necessary variables + const auto &prim_pack = md->PackVariables(std::vector{"prim"}); + IndexRange ib = md->GetBlockData(0)->GetBoundsI(IndexDomain::interior); + IndexRange jb = md->GetBlockData(0)->GetBoundsJ(IndexDomain::interior); + IndexRange kb = md->GetBlockData(0)->GetBoundsK(IndexDomain::interior); + + Real min_cooling_time = std::numeric_limits::infinity(); + Kokkos::Min reducer_min(min_cooling_time); + + Kokkos::parallel_reduce( + "TabularCooling::TimeStep", + Kokkos::MDRangePolicy>( + {0, kb.s, jb.s, ib.s}, {prim_pack.GetDim(5), kb.e + 1, jb.e + 1, ib.e + 1}, + {1, 1, 1, ib.e + 1 - ib.s}), + KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i, + Real &thread_min_cooling_time) { + auto &prim = prim_pack(b); + + const Real rho = prim(IDN, k, j, i); + const Real pres = prim(IPR, k, j, i); + + const Real internal_e = pres / (rho * gm1); + + const Real de_dt = cooling_table_obj.DeDt(internal_e, rho); + + // Compute cooling time + // If de_dt is zero (temperature is smaller than lower end of cooling table) or + // current temp is below floor, use infinite cooling time + const Real cooling_time = ((de_dt == 0) || (internal_e < internal_e_floor)) + ? std::numeric_limits::infinity() + : fabs(internal_e / de_dt); + + thread_min_cooling_time = std::min(cooling_time, thread_min_cooling_time); + }, + reducer_min); + + return cooling_time_cfl_ * min_cooling_time; +} + +void TabularCooling::TestCoolingTable(ParameterInput *pin) const { + + const std::string test_filename = pin->GetString("cooling", "test_filename"); + + const auto rho0 = pin->GetReal("cooling", "test_rho0"); + const auto rho1 = pin->GetReal("cooling", "test_rho1"); + const auto n_rho = pin->GetInteger("cooling", "test_n_rho"); + + const auto pres0 = pin->GetReal("cooling", "test_pres0"); + const auto pres1 = pin->GetReal("cooling", "test_pres1"); + const auto n_pres = pin->GetInteger("cooling", "test_n_pres"); + + // Grab member variables for compiler + + // Everything needed by DeDt + const CoolingTableObj cooling_table_obj = cooling_table_obj_; + const auto gm1 = pin->GetReal("hydro", "gamma") - 1.0; + + // Make some device arrays to store the test data + ParArray2D d_rho("d_rho", n_rho, n_pres), d_pres("d_pres", n_rho, n_pres), + d_internal_e("d_internal_e", n_rho, n_pres), d_de_dt("d_de_dt", n_rho, n_pres); + + par_for( + loop_pattern_mdrange_tag, "TabularCooling::TestCoolingTable", DevExecSpace(), 0, + n_rho - 1, 0, n_pres - 1, KOKKOS_LAMBDA(const int &j, const int &i) { + const Real rho = rho0 * pow(rho1 / rho0, static_cast(j) / (n_rho - 1)); + const Real pres = pres0 * pow(pres1 / pres0, static_cast(i) / (n_pres - 1)); + + d_rho(j, i) = rho; + d_pres(j, i) = pres; + + const Real internal_e = pres / (rho * gm1); + + d_internal_e(j, i) = internal_e; + + const Real de_dt = cooling_table_obj.DeDt(internal_e, rho); + + d_de_dt(j, i) = de_dt; + }); + + // Copy Device arrays to host + auto h_rho = Kokkos::create_mirror_view_and_copy(HostMemSpace(), d_rho); + auto h_pres = Kokkos::create_mirror_view_and_copy(HostMemSpace(), d_pres); + auto h_internal_e = Kokkos::create_mirror_view_and_copy(HostMemSpace(), d_internal_e); + auto h_de_dt = Kokkos::create_mirror_view_and_copy(HostMemSpace(), d_de_dt); + + // Write to file + std::ofstream file(test_filename); + file << "#rho pres internal_e de_dt" << std::endl; + for (int j = 0; j < n_rho; j++) { + for (int i = 0; i < n_pres; i++) { + file << h_rho(j, i) << " " << h_pres(j, i) << " " << h_internal_e(j, i) << " " + << h_de_dt(j, i) << " " << std::endl; + } + } +} + +} // namespace cooling diff --git a/src/hydro/srcterms/old/tabular_cooling_edit.hpp b/src/hydro/srcterms/old/tabular_cooling_edit.hpp new file mode 100644 index 00000000..5caea98e --- /dev/null +++ b/src/hydro/srcterms/old/tabular_cooling_edit.hpp @@ -0,0 +1,314 @@ +//======================================================================================== +// AthenaPK - a performance portable block structured AMR astrophysical MHD code. +// Copyright (c) 2021, Athena-Parthenon Collaboration. All rights reserved. +// Licensed under the 3-clause BSD License, see LICENSE file for details +//======================================================================================== +//! \file tabular_cooling.hpp +// \brief Class for defining tabular cooling +#ifndef HYDRO_SRCTERMS_TABULAR_COOLING_HPP_ +#define HYDRO_SRCTERMS_TABULAR_COOLING_HPP_ + +// C++ headers +#include // stringstream +#include // istream_iterator +#include // stringstream +#include // runtime_error +#include // string +#include // vector + +// Parthenon headers +#include +#include +#include +#include +#include + +// AthenaPK headers +#include "../../main.hpp" +#include "../../units.hpp" + +#ifdef MPI_PARALLEL +#include +#endif + +namespace cooling { + +// Struct to take one RK step using heun's method to compute 2nd and 1st order estimations +// in y1_h and y1_l +struct RK12Stepper { + template + static KOKKOS_INLINE_FUNCTION void + Step(const parthenon::Real t0, const parthenon::Real h, const parthenon::Real y0, + Function f, parthenon::Real &y1_h, parthenon::Real &y1_l, bool &valid) { + const parthenon::Real f_t0_y0 = f(t0, y0, valid); + y1_l = y0 + h * f_t0_y0; // 1st order + y1_h = y0 + h / 2. * (f_t0_y0 + f(t0 + h, y1_l, valid)); // 2nd order + } + static KOKKOS_INLINE_FUNCTION parthenon::Real OptimalStep(const parthenon::Real h, + const parthenon::Real err, + const parthenon::Real tol) { + return 0.95 * h * pow(tol / err, 2); + } +}; + +// Struct to take a 5th and 4th order RK step to compute 5th and 4th order estimations in +// y1_h and y1_l +struct RK45Stepper { + template + static KOKKOS_INLINE_FUNCTION void + Step(const parthenon::Real t0, const parthenon::Real h, const parthenon::Real y0, + Function f, parthenon::Real &y1_h, parthenon::Real &y1_l, bool &valid) { + const parthenon::Real k1 = h * f(t0, y0, valid); + const parthenon::Real k2 = h * f(t0 + 1. / 4. * h, y0 + 1. / 4. * k1, valid); + const parthenon::Real k3 = + h * f(t0 + 3. / 8. * h, y0 + 3. / 32. * k1 + 9. / 32. * k2, valid); + const parthenon::Real k4 = + h * f(t0 + 12. / 13. * h, + y0 + 1932. / 2197. * k1 - 7200. / 2197. * k2 + 7296. / 2197. * k3, valid); + const parthenon::Real k5 = + h * f(t0 + h, + y0 + 439. / 216. * k1 - 8. * k2 + 3680. / 513. * k3 - 845. / 4104. * k4, + valid); + const parthenon::Real k6 = h * f(t0 + 1. / 2. * h, + y0 - 8. / 27. * k1 + 2. * k2 - 3544. / 2565. * k3 + + 1859. / 4104. * k4 - 11. / 40. * k5, + valid); // TODO(forrestglines): Check k2? + y1_l = y0 + 25. / 216. * k1 + 1408. / 2565. * k3 + 2197. / 4104. * k4 - + 1. / 5. * k5; // 4th order + y1_h = y0 + 16. / 135. * k1 + 6656. / 12825. * k3 + 28561. / 56430. * k4 - + 9. / 50. * k5 + 2. / 55. * k6; // 5th order + } + static KOKKOS_INLINE_FUNCTION parthenon::Real OptimalStep(const parthenon::Real h, + const parthenon::Real err, + const parthenon::Real tol) { + return 0.95 * h * pow(tol / err, 5); + } +}; + +enum class CoolIntegrator { undefined, rk12, rk45, townsend }; + +class CoolingTableObj { + /************************************************************ + * Cooling Table Object, for interpolating a cooling rate out of a cooling + * table. Currently assumes evenly space log_temperatures in cooling table + * + * Lightweight object intended for inlined computation within kernels + ************************************************************/ + private: + // Log cooling rate/ne^3 + parthenon::ParArray1D log_lambdas_; + + // Spacing of cooling table + // TODO: assumes evenly spaced cooling table + parthenon::Real log_temp_start_, log_temp_final_, d_log_temp_; + unsigned int n_temp_; + + // Mean molecular mass * ( adiabatic_index -1) / boltzmann_constant + parthenon::Real mbar_gm1_over_k_B_; + + // (Hydrogen mass fraction / hydrogen atomic mass)^2 + parthenon::Real x_H_over_m_h2_; + + public: + CoolingTableObj() + : log_lambdas_(), log_temp_start_(NAN), log_temp_final_(NAN), d_log_temp_(NAN), + n_temp_(0), mbar_gm1_over_k_B_(NAN), x_H_over_m_h2_(NAN) {} + CoolingTableObj(const parthenon::ParArray1D log_lambdas, + const parthenon::Real log_temp_start, + const parthenon::Real log_temp_final, const parthenon::Real d_log_temp, + const unsigned int n_temp, const parthenon::Real mbar_over_kb, + const parthenon::Real adiabatic_index, const parthenon::Real x_H, + const Units units) + : log_lambdas_(log_lambdas), log_temp_start_(log_temp_start), + log_temp_final_(log_temp_final), d_log_temp_(d_log_temp), n_temp_(n_temp), + mbar_gm1_over_k_B_(mbar_over_kb * (adiabatic_index - 1)), + x_H_over_m_h2_(SQR(x_H / units.mh())) {} + + // Interpolate a cooling rate from the table + // from internal energy density and density + KOKKOS_INLINE_FUNCTION parthenon::Real + DeDt(const parthenon::Real &e, const parthenon::Real &rho, bool &is_valid) const { + using namespace parthenon; + + if (e < 0 || std::isnan(e)) { + is_valid = false; + return 0; + } + + const Real temp = mbar_gm1_over_k_B_ * e; + const Real log_temp = log10(temp); + Real log_lambda; + if (log_temp < log_temp_start_) { + return 0; // Return no variation + + } else if (log_temp > log_temp_final_) { + // Cooling insured by free-free cooling ie. Bremsstrahlung + // Above table + // Return de/dt + // TODO(forrestglines):Currently free-free cooling is used for + // temperatures above the table. This behavior could be generalized via + // templates + log_lambda = 0.5 * log_temp - 0.5 * log_temp_final_ + log_lambdas_(n_temp_ - 1); + + // NEED TO BE MODIFIED HERE: ADD THE HEATING FUNCTION AND CHECK FOR CONTINUITY + + } else { + // Inside table, interpolate assuming log spaced temperatures + + // Determine where temp is in the table + const unsigned int i_temp = + static_cast((log_temp - log_temp_start_) / d_log_temp_); + const Real log_temp_i = log_temp_start_ + d_log_temp_ * i_temp; + + // log_temp should be between log_temps[i_temp] and log_temps[i_temp+1] + PARTHENON_REQUIRE(log_temp >= log_temp_i && log_temp <= log_temp_i + d_log_temp_, + "FATAL ERROR in [CoolingTable::DeDt]: Failed to find log_temp"); + + const Real log_lambda_i = log_lambdas_(i_temp); + const Real log_lambda_ip1 = log_lambdas_(i_temp + 1); + + // Linearly interpolate lambda at log_temp + log_lambda = log_lambda_i + (log_temp - log_temp_i) * + (log_lambda_ip1 - log_lambda_i) / d_log_temp_; + } + // Return de/dt + const Real lambda = pow(10., log_lambda); + const Real de_dt = -lambda * x_H_over_m_h2_ * rho; + return de_dt; + } + + // Interpolate a cooling rate from the table + // from internal energy density and density + + KOKKOS_INLINE_FUNCTION parthenon::Real + DeDt(const parthenon::Real &e, const parthenon::Real &rho, bool &is_valid) const { + using namespace parthenon; + + if (e < 0 || std::isnan(e)) { + is_valid = false; + return 0; + } + + const Real temp = mbar_gm1_over_k_B_ * e; + const Real log_temp = log10(temp); + Real log_lambda; + if (log_temp < log_temp_start_) { + return 0; // Return no variation + + } else if (log_temp > log_temp_final_) { + // Cooling insured by free-free cooling ie. Bremsstrahlung + // Above table + // Return de/dt + // TODO(forrestglines):Currently free-free cooling is used for + // temperatures above the table. This behavior could be generalized via + // templates + log_lambda = 0.5 * log_temp - 0.5 * log_temp_final_ + log_lambdas_(n_temp_ - 1); + + // NEED TO BE MODIFIED HERE: ADD THE HEATING FUNCTION AND CHECK FOR CONTINUITY + + } else { + // Inside table, interpolate assuming log spaced temperatures + + // Determine where temp is in the table + const unsigned int i_temp = + static_cast((log_temp - log_temp_start_) / d_log_temp_); + const Real log_temp_i = log_temp_start_ + d_log_temp_ * i_temp; + + // log_temp should be between log_temps[i_temp] and log_temps[i_temp+1] + PARTHENON_REQUIRE(log_temp >= log_temp_i && log_temp <= log_temp_i + d_log_temp_, + "FATAL ERROR in [CoolingTable::DeDt]: Failed to find log_temp"); + + const Real log_lambda_i = log_lambdas_(i_temp); + const Real log_lambda_ip1 = log_lambdas_(i_temp + 1); + + // Linearly interpolate lambda at log_temp + log_lambda = log_lambda_i + (log_temp - log_temp_i) * + (log_lambda_ip1 - log_lambda_i) / d_log_temp_; + } + // Return de/dt + const Real lambda = pow(10., log_lambda); + const Real de_dt = -lambda * x_H_over_m_h2_ * rho; + return de_dt; + } + + + + KOKKOS_INLINE_FUNCTION parthenon::Real DeDt(const parthenon::Real &e, + const parthenon::Real &rho) const { + bool is_valid = true; + return DeDt(e, rho, is_valid); + } +}; + +class TabularCooling { + private: + // Defines uniformly spaced log temperature range of the table + unsigned int n_temp_; + parthenon::Real log_temp_start_, log_temp_final_, d_log_temp_, lambda_final_; + + // Table of log cooling rates + // TODO(forrestglines): Make log_lambdas_ explicitly a texture cache array, use CUDA to + // interpolate directly + // Log versions are used in subcyling cooling where cooling rates are interpolated + // Non-log versions are used for Townsend cooling + parthenon::ParArray1D log_lambdas_; + parthenon::ParArray1D lambdas_; + parthenon::ParArray1D temps_; + // Townsend cooling temporal evolution function + parthenon::ParArray1D townsend_Y_k_; + // Townsend cooling power law indices + parthenon::ParArray1D townsend_alpha_k_; + + CoolIntegrator integrator_; + + // Temperature floor (assumed in Kelvin and only used in cooling function) + // This is either the temperature floor used by the hydro method or the + // lowest temperature in the cooling table (assuming zero cooling below the + // table), whichever temperature is higher + parthenon::Real T_floor_; + + // Maximum number of iterations/subcycles + unsigned int max_iter_; + + // Cooling CFL + parthenon::Real cooling_time_cfl_; + + // Minimum timestep that the cooling may limit the simulation timestep + // Use nonpositive values to disable + parthenon::Real min_cooling_timestep_; + + // Tolerances + parthenon::Real d_log_temp_tol_, d_e_tol_; + + // Used for roundoff as subcycle approaches end of timestep + static constexpr parthenon::Real KEpsilon_ = 1e-12; + + CoolingTableObj cooling_table_obj_; + + public: + TabularCooling(parthenon::ParameterInput *pin, + std::shared_ptr hydro_pkg); + + void SrcTerm(parthenon::MeshData *md, const parthenon::Real dt) const; + + // Townsend 2009 exact integration scheme + void TownsendSrcTerm(parthenon::MeshData *md, + const parthenon::Real dt) const; + + // (Adaptive) subcyling using a fixed integration scheme + template + void SubcyclingFixedIntSrcTerm(parthenon::MeshData *md, + const parthenon::Real dt, + const RKStepper rk_stepper) const; + + parthenon::Real EstimateTimeStep(parthenon::MeshData *md) const; + + // Get a lightweight object for computing cooling rate from the cooling table + const CoolingTableObj GetCoolingTableObj() const { return cooling_table_obj_; } + + void TestCoolingTable(parthenon::ParameterInput *pin) const; +}; + +} // namespace cooling + +#endif // HYDRO_SRCTERMS_TABULAR_COOLING_HPP_ diff --git a/src/hydro/srcterms/tabular_cooling.cpp b/src/hydro/srcterms/tabular_cooling.cpp index 9164f4cc..ee4645b5 100644 --- a/src/hydro/srcterms/tabular_cooling.cpp +++ b/src/hydro/srcterms/tabular_cooling.cpp @@ -36,17 +36,29 @@ TabularCooling::TabularCooling(ParameterInput *pin, const int log_temp_col = pin->GetOrAddInteger("cooling", "log_temp_col", 0); const int log_lambda_col = pin->GetOrAddInteger("cooling", "log_lambda_col", 1); + // Heating constants + const Real heating_const_cgs = pin->GetOrAddReal("cooling", "heating_const", 0.0); + const Real gamma_units = heating_const_cgs * (units.erg() / units.s()); + + // Checkpoint + std::cout << "heating_const_cgs=" << heating_const_cgs << "\n"; + std::cout << "units.s()=" << units.s() << "\n"; + std::cout << "gamma_units=" << gamma_units << "\n"; + const Real lambda_units_cgs = pin->GetReal("cooling", "lambda_units_cgs"); // Convert erg cm^3/s to code units const Real lambda_units = lambda_units_cgs / (units.erg() * pow(units.cm(), 3) / units.s()); - + const auto integrator_str = pin->GetOrAddString("cooling", "integrator", "rk12"); if (integrator_str == "rk12") { + std::cout << "Cooling integrator is rk12\n"; integrator_ = CoolIntegrator::rk12; } else if (integrator_str == "rk45") { + std::cout << "Cooling integrator is rk45\n"; integrator_ = CoolIntegrator::rk45; - } else if (integrator_str == "townsend") { + } else if (integrator_str == "townsend\n") { + std::cout << "Cooling integrator is Townsend"; integrator_ = CoolIntegrator::townsend; } else { integrator_ = CoolIntegrator::undefined; @@ -115,43 +127,49 @@ TabularCooling::TabularCooling(ParameterInput *pin, << line << "\"" << std::endl; PARTHENON_FAIL(msg); } - + try { const Real log_temp = std::stod(line_data[log_temp_col]); const Real log_lambda = std::stod(line_data[log_lambda_col]); - + // Add to growing list log_temps.push_back(log_temp); log_lambdas.push_back(log_lambda - std::log10(lambda_units)); - + } catch (const std::invalid_argument &ia) { msg << "### FATAL ERROR in function [TabularCooling::TabularCooling]" << std::endl << "Number: \"" << ia.what() << "\" could not be parsed as double" << std::endl; PARTHENON_FAIL(msg); } } - + /**************************************** * Check some assumtions about the cooling table ****************************************/ - + + // Reconstructing the cooling function in new coordinate system + //std::cout << "Entering lambda cout loop\n"; + //for (int i = 0; i < log_temps.size(); i++) { + //std::cout << "log_temps" << log_temps[i] << ", log_lambdas" << log_lambdas[i] << "\n"; + //} + // Ensure at least two data points in the table to interpolate from if (log_temps.size() < 2 || log_lambdas.size() < 2) { msg << "### FATAL ERROR in function [TabularCooling::TabularCooling]" << std::endl << "Not enough data to interpolate cooling" << std::endl; PARTHENON_FAIL(msg); } - + // Ensure that the first log_temp is increasing const Real log_temp_start = log_temps[0]; const Real d_log_temp = log_temps[1] - log_temp_start; - + if (d_log_temp <= 0) { msg << "### FATAL ERROR in function [TabularCooling::TabularCooling]" << std::endl << "second log_temp in table is descreasing" << std::endl; PARTHENON_FAIL(msg); } - + // Ensure that log_temps is evenly spaced for (size_t i = 1; i < log_temps.size(); i++) { const Real d_log_temp_i = log_temps[i] - log_temps[i - 1]; @@ -162,7 +180,7 @@ TabularCooling::TabularCooling(ParameterInput *pin, << " log_temp= " << log_temps[i] << std::endl; PARTHENON_FAIL(msg); } - + // The Dedt() function currently relies on an equally spaced cooling table for faster // lookup (direct indexing). It is used in the subcycling method and in order to // restrict `dt` by a cooling cfl. @@ -177,17 +195,18 @@ TabularCooling::TabularCooling(ParameterInput *pin, PARTHENON_FAIL(msg); } } - + /**************************************** * Move values read into the data table ****************************************/ - + n_temp_ = log_temps.size(); log_temp_start_ = log_temps[0]; log_temp_final_ = log_temps[n_temp_ - 1]; d_log_temp_ = d_log_temp; lambda_final_ = std::pow(10.0, log_lambdas[n_temp_ - 1]); - + gamma_units_ = gamma_units; + // Setup log_lambdas_ used in Dedt() { // log_lambdas is used if the integrator isn't Townsend, if the cooling CFL @@ -215,6 +234,7 @@ TabularCooling::TabularCooling(ParameterInput *pin, host_lambdas(i) = std::pow(10.0, log_lambdas[i]); host_temps(i) = std::pow(10.0, log_temps[i]); } + // Copy host_lambdas into device memory Kokkos::deep_copy(lambdas_, host_lambdas); Kokkos::deep_copy(temps_, host_temps); @@ -223,11 +243,11 @@ TabularCooling::TabularCooling(ParameterInput *pin, const auto n_bins = n_temp_ - 1; townsend_Y_k_ = ParArray1D("townsend_Y_k_", n_bins); townsend_alpha_k_ = ParArray1D("townsend_alpha_k_", n_bins); - + // Initialize on host (make *this* recursion simpler) auto host_townsend_Y_k = Kokkos::create_mirror_view(townsend_Y_k_); auto host_townsend_alpha_k = Kokkos::create_mirror_view(townsend_alpha_k_); - + // Initialize piecewise power law indices for (unsigned int i = 0; i < n_bins; i++) { // Could use log_lambdas_ here, but using lambdas_ instead as they've already been @@ -262,7 +282,7 @@ TabularCooling::TabularCooling(ParameterInput *pin, const auto He_mass_fraction = hydro_pkg->Param("He_mass_fraction"); cooling_table_obj_ = CoolingTableObj(log_lambdas_, log_temp_start_, log_temp_final_, - d_log_temp_, n_temp_, mbar_over_kb, + d_log_temp_, n_temp_, gamma_units_, mbar_over_kb, adiabatic_index, 1.0 - He_mass_fraction, units); } @@ -290,20 +310,20 @@ void TabularCooling::SubcyclingFixedIntSrcTerm(MeshData *md, const Real dt const CoolingTableObj cooling_table_obj = cooling_table_obj_; const auto gm1 = (hydro_pkg->Param("AdiabaticIndex") - 1.0); const auto mbar_gm1_over_kb = hydro_pkg->Param("mbar_over_kb") * gm1; - + const unsigned int max_iter = max_iter_; - + const Real min_sub_dt = dt / max_iter; - + const Real d_e_tol = d_e_tol_; - + // Determine the cooling floor, whichever is higher of the cooling table floor // or fluid solver floor const auto temp_cool_floor = std::pow(10.0, log_temp_start_); // low end of cool table const Real temp_floor = (T_floor_ > temp_cool_floor) ? T_floor_ : temp_cool_floor; - + const Real internal_e_floor = temp_floor / mbar_gm1_over_kb; // specific internal en. - + // Grab some necessary variables const auto &prim_pack = md->PackVariables(std::vector{"prim"}); const auto &cons_pack = md->PackVariables(std::vector{"cons"}); @@ -312,7 +332,7 @@ void TabularCooling::SubcyclingFixedIntSrcTerm(MeshData *md, const Real dt IndexRange ib = md->GetBlockData(0)->GetBoundsI(IndexDomain::entire); IndexRange jb = md->GetBlockData(0)->GetBoundsJ(IndexDomain::entire); IndexRange kb = md->GetBlockData(0)->GetBoundsK(IndexDomain::entire); - + par_for( DEFAULT_LOOP_PATTERN, "TabularCooling::SubcyclingSplitSrcTerm", DevExecSpace(), 0, cons_pack.GetDim(5) - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, diff --git a/src/hydro/srcterms/tabular_cooling.hpp b/src/hydro/srcterms/tabular_cooling.hpp index c1ffa8b2..bea2499b 100644 --- a/src/hydro/srcterms/tabular_cooling.hpp +++ b/src/hydro/srcterms/tabular_cooling.hpp @@ -108,19 +108,25 @@ class CoolingTableObj { // (Hydrogen mass fraction / hydrogen atomic mass)^2 parthenon::Real x_H_over_m_h2_; + + // Heating constant unit + parthenon::Real gamma_units_; public: CoolingTableObj() : log_lambdas_(), log_temp_start_(NAN), log_temp_final_(NAN), d_log_temp_(NAN), - n_temp_(0), mbar_gm1_over_k_B_(NAN), x_H_over_m_h2_(NAN) {} + n_temp_(0), mbar_gm1_over_k_B_(NAN), x_H_over_m_h2_(NAN), gamma_units_(NAN) {} CoolingTableObj(const parthenon::ParArray1D log_lambdas, const parthenon::Real log_temp_start, const parthenon::Real log_temp_final, const parthenon::Real d_log_temp, - const unsigned int n_temp, const parthenon::Real mbar_over_kb, + const unsigned int n_temp, + const parthenon::Real gamma_units, + const parthenon::Real mbar_over_kb, const parthenon::Real adiabatic_index, const parthenon::Real x_H, const Units units) : log_lambdas_(log_lambdas), log_temp_start_(log_temp_start), log_temp_final_(log_temp_final), d_log_temp_(d_log_temp), n_temp_(n_temp), + gamma_units_(gamma_units), mbar_gm1_over_k_B_(mbar_over_kb * (adiabatic_index - 1)), x_H_over_m_h2_(SQR(x_H / units.mh())) {} @@ -129,12 +135,12 @@ class CoolingTableObj { KOKKOS_INLINE_FUNCTION parthenon::Real DeDt(const parthenon::Real &e, const parthenon::Real &rho, bool &is_valid) const { using namespace parthenon; - + if (e < 0 || std::isnan(e)) { is_valid = false; return 0; } - + const Real temp = mbar_gm1_over_k_B_ * e; const Real log_temp = log10(temp); Real log_lambda; @@ -168,7 +174,11 @@ class CoolingTableObj { } // Return de/dt const Real lambda = pow(10., log_lambda); - const Real de_dt = -lambda * x_H_over_m_h2_ * rho; + const Real gamma = gamma_units_; + const Real de_dt = -lambda * x_H_over_m_h2_ * rho + gamma*pow(x_H_over_m_h2_,0.5); + + //std::cout << "cool = " << -lambda * x_H_over_m_h2_ * rho << "\t heat = " << gamma*pow(x_H_over_m_h2_,0.5) << "\t e = " << e << "\t rho = " << rho << "\n" ; + return de_dt; } @@ -183,7 +193,7 @@ class TabularCooling { private: // Defines uniformly spaced log temperature range of the table unsigned int n_temp_; - parthenon::Real log_temp_start_, log_temp_final_, d_log_temp_, lambda_final_; + parthenon::Real log_temp_start_, log_temp_final_, d_log_temp_, lambda_final_, gamma_units_; // Table of log cooling rates // TODO(forrestglines): Make log_lambdas_ explicitly a texture cache array, use CUDA to diff --git a/src/hydro/srcterms/tabular_cooling_working.cpp b/src/hydro/srcterms/tabular_cooling_working.cpp new file mode 100644 index 00000000..31ed7f3e --- /dev/null +++ b/src/hydro/srcterms/tabular_cooling_working.cpp @@ -0,0 +1,733 @@ +//======================================================================================== +// AthenaPK - a performance portable block structured AMR astrophysical MHD code. +// Copyright (c) 2021-2023, Athena-Parthenon Collaboration. All rights reserved. +// Licensed under the 3-clause BSD License, see LICENSE file for details +//======================================================================================== +//! \file tabular_cooling.cpp +// \brief Applies tabular cooling +// +//======================================================================================== + +// C++ headers +#include +#include +#include + +// Parthenon headers +#include +#include +#include +#include + +// AthenaPK headers +#include "../../units.hpp" +#include "tabular_cooling.hpp" +#include "utils/error_checking.hpp" + +namespace cooling { +using namespace parthenon; + +TabularCooling::TabularCooling(ParameterInput *pin, + std::shared_ptr hydro_pkg) { + auto units = hydro_pkg->Param("units"); + + const std::string table_filename = pin->GetString("cooling", "table_filename"); + + const int log_temp_col = pin->GetOrAddInteger("cooling", "log_temp_col", 0); + const int log_lambda_col = pin->GetOrAddInteger("cooling", "log_lambda_col", 1); + + // Heating constants + const Real heating_const_cgs = pin->GetOrAddReal("cooling", "heating_const", 0.0); + const Real gamma_units = heating_const_cgs * (units.erg() / units.s()); + + // Checkpoint + std::cout << "heating_const_cgs=" << heating_const_cgs << "\n"; + std::cout << "units.s()=" << units.s() << "\n"; + std::cout << "gamma_units=" << gamma_units << "\n"; + + const Real lambda_units_cgs = pin->GetReal("cooling", "lambda_units_cgs"); + // Convert erg cm^3/s to code units + const Real lambda_units = + lambda_units_cgs / (units.erg() * pow(units.cm(), 3) / units.s()); + + const auto integrator_str = pin->GetOrAddString("cooling", "integrator", "rk12"); + if (integrator_str == "rk12") { + std::cout << "Cooling integrator is rk12\n"; + integrator_ = CoolIntegrator::rk12; + } else if (integrator_str == "rk45") { + std::cout << "Cooling integrator is rk45\n"; + integrator_ = CoolIntegrator::rk45; + } else if (integrator_str == "townsend\n") { + std::cout << "Cooling integrator is Townsend"; + integrator_ = CoolIntegrator::townsend; + } else { + integrator_ = CoolIntegrator::undefined; + } + max_iter_ = pin->GetOrAddInteger("cooling", "max_iter", 100); + cooling_time_cfl_ = pin->GetOrAddReal("cooling", "cfl", 0.1); + d_log_temp_tol_ = pin->GetOrAddReal("cooling", "d_log_temp_tol", 1e-8); + d_e_tol_ = pin->GetOrAddReal("cooling", "d_e_tol", 1e-8); + // negative means disabled + T_floor_ = pin->GetOrAddReal("hydro", "Tfloor", -1.0); + + std::stringstream msg; + + /**************************************** + * Read tab file with IOWrapper + ****************************************/ + IOWrapper input; + input.Open(table_filename.c_str(), IOWrapper::FileMode::read); + + /**************************************** + * Read tab file from IOWrapper into a stringstream tab + ****************************************/ + std::stringstream tab_ss; + const int bufsize = 4096; + char *buf = new char[bufsize]; + std::ptrdiff_t ret; + parthenon::IOWrapperSizeT word_size = sizeof(char); + + do { + if (Globals::my_rank == 0) { // only the master process reads the cooling table + ret = input.Read(buf, word_size, bufsize); + } +#ifdef MPI_PARALLEL + // then broadcasts it + // no need for fence as cooling table is independent of execution/memory space + MPI_Bcast(&ret, sizeof(std::ptrdiff_t), MPI_BYTE, 0, MPI_COMM_WORLD); + MPI_Bcast(buf, ret, MPI_BYTE, 0, MPI_COMM_WORLD); +#endif + tab_ss.write(buf, ret); // add the buffer into the stream + } while (ret == bufsize); // till EOF (or par_end is found) + + delete[] buf; + input.Close(); + + /**************************************** + * Determine log_temps and and log_lambdas vectors + ****************************************/ + std::vector log_temps, log_lambdas; + std::string line; + std::size_t first_char; + while (tab_ss.good()) { + getline(tab_ss, line); + if (line.empty()) continue; // skip blank line + first_char = line.find_first_not_of(" "); // skip white space + if (first_char == std::string::npos) continue; // line is all white space + if (line.compare(first_char, 1, "#") == 0) continue; // skip comments + + // Parse the numbers on the line + std::istringstream iss(line); + std::vector line_data{std::istream_iterator{iss}, + std::istream_iterator{}}; + // Check size + if (line_data.empty() || line_data.size() <= std::max(log_temp_col, log_lambda_col)) { + msg << "### FATAL ERROR in function [TabularCooling::TabularCooling]" << std::endl + << "Index " << std::max(log_temp_col, log_lambda_col) << " out of range on \"" + << line << "\"" << std::endl; + PARTHENON_FAIL(msg); + } + + try { + const Real log_temp = std::stod(line_data[log_temp_col]); + const Real log_lambda = std::stod(line_data[log_lambda_col]); + + // Add to growing list + log_temps.push_back(log_temp); + log_lambdas.push_back(log_lambda - std::log10(lambda_units)); + + } catch (const std::invalid_argument &ia) { + msg << "### FATAL ERROR in function [TabularCooling::TabularCooling]" << std::endl + << "Number: \"" << ia.what() << "\" could not be parsed as double" << std::endl; + PARTHENON_FAIL(msg); + } + } + + /**************************************** + * Check some assumtions about the cooling table + ****************************************/ + + // Reconstructing the cooling function in new coordinate system + //std::cout << "Entering lambda cout loop\n"; + //for (int i = 0; i < log_temps.size(); i++) { + //std::cout << "log_temps" << log_temps[i] << ", log_lambdas" << log_lambdas[i] << "\n"; + //} + + // Ensure at least two data points in the table to interpolate from + if (log_temps.size() < 2 || log_lambdas.size() < 2) { + msg << "### FATAL ERROR in function [TabularCooling::TabularCooling]" << std::endl + << "Not enough data to interpolate cooling" << std::endl; + PARTHENON_FAIL(msg); + } + + // Ensure that the first log_temp is increasing + const Real log_temp_start = log_temps[0]; + const Real d_log_temp = log_temps[1] - log_temp_start; + + if (d_log_temp <= 0) { + msg << "### FATAL ERROR in function [TabularCooling::TabularCooling]" << std::endl + << "second log_temp in table is descreasing" << std::endl; + PARTHENON_FAIL(msg); + } + + // Ensure that log_temps is evenly spaced + for (size_t i = 1; i < log_temps.size(); i++) { + const Real d_log_temp_i = log_temps[i] - log_temps[i - 1]; + + if (d_log_temp_i < 0) { + msg << "### FATAL ERROR in function [TabularCooling::TabularCooling]" << std::endl + << "log_temp in table is descreasing at i= " << i + << " log_temp= " << log_temps[i] << std::endl; + PARTHENON_FAIL(msg); + } + + // The Dedt() function currently relies on an equally spaced cooling table for faster + // lookup (direct indexing). It is used in the subcycling method and in order to + // restrict `dt` by a cooling cfl. + if (((integrator_ != CoolIntegrator::townsend) || (cooling_time_cfl_ > 0.0)) && + (fabs(d_log_temp_i - d_log_temp) / d_log_temp > d_log_temp_tol_)) { + msg << "### FATAL ERROR in function [TabularCooling::TabularCooling]" << std::endl + << "d_log_temp in table is uneven at i=" << i << " log_temp=" << log_temps[i] + << " d_log_temp= " << d_log_temp << " d_log_temp_i= " << d_log_temp_i + << " diff= " << d_log_temp_i - d_log_temp + << " rel_diff= " << fabs(d_log_temp_i - d_log_temp) / d_log_temp + << " tol= " << d_log_temp_tol_ << std::endl; + PARTHENON_FAIL(msg); + } + } + + /**************************************** + * Move values read into the data table + ****************************************/ + + n_temp_ = log_temps.size(); + log_temp_start_ = log_temps[0]; + log_temp_final_ = log_temps[n_temp_ - 1]; + d_log_temp_ = d_log_temp; + lambda_final_ = std::pow(10.0, log_lambdas[n_temp_ - 1]); + + // Setup log_lambdas_ used in Dedt() + { + // log_lambdas is used if the integrator isn't Townsend, if the cooling CFL + // is set, or if cooling time is a extra derived field. Since we don't have + // a good way to check the last condition we always initialize log_lambdas_ + log_lambdas_ = ParArray1D("log_lambdas_", n_temp_); + + // Read log_lambdas in host_log_lambdas, changing to code units along the way + auto host_log_lambdas = Kokkos::create_mirror_view(log_lambdas_); + for (unsigned int i = 0; i < n_temp_; i++) { + host_log_lambdas(i) = log_lambdas[i]; + } + // Copy host_log_lambdas into device memory + Kokkos::deep_copy(log_lambdas_, host_log_lambdas); + } + // Setup Townsend cooling, i.e., precalulcate piecewise powerlaw approx. + if (integrator_ == CoolIntegrator::townsend) { + lambdas_ = ParArray1D("lambdas_", n_temp_); + temps_ = ParArray1D("temps_", n_temp_); + + // Read log_lambdas in host_lambdas, changing to code units along the way + auto host_lambdas = Kokkos::create_mirror_view(lambdas_); + auto host_temps = Kokkos::create_mirror_view(temps_); + for (unsigned int i = 0; i < n_temp_; i++) { + host_lambdas(i) = std::pow(10.0, log_lambdas[i]); + host_temps(i) = std::pow(10.0, log_temps[i]); + } + // Copy host_lambdas into device memory + Kokkos::deep_copy(lambdas_, host_lambdas); + Kokkos::deep_copy(temps_, host_temps); + + // Coeffs are for intervals, i.e., only n_temp_ - 1 entries + const auto n_bins = n_temp_ - 1; + townsend_Y_k_ = ParArray1D("townsend_Y_k_", n_bins); + townsend_alpha_k_ = ParArray1D("townsend_alpha_k_", n_bins); + + // Initialize on host (make *this* recursion simpler) + auto host_townsend_Y_k = Kokkos::create_mirror_view(townsend_Y_k_); + auto host_townsend_alpha_k = Kokkos::create_mirror_view(townsend_alpha_k_); + + // Initialize piecewise power law indices + for (unsigned int i = 0; i < n_bins; i++) { + // Could use log_lambdas_ here, but using lambdas_ instead as they've already been + // converted to code units. + host_townsend_alpha_k(i) = + (std::log10(host_lambdas(i + 1)) - std::log10(host_lambdas(i))) / + (log_temps[i + 1] - log_temps[i]); + PARTHENON_REQUIRE(host_townsend_alpha_k(i) != 1.0, + "Need to implement special case for Townsend piecewise fits."); + } + + // Calculate TEF (temporal evolution functions Y_k recursively), (Eq. A6) + host_townsend_Y_k(n_bins - 1) = 0.0; // Last Y_N = Y(T_ref) = 0 + + for (int i = n_bins - 2; i >= 0; i--) { + const auto alpha_k_m1 = host_townsend_alpha_k(i) - 1.0; + const auto step = (host_lambdas(n_bins) / host_lambdas(i)) * + (host_temps(i) / host_temps(n_bins)) * + (std::pow(host_temps(i) / host_temps(i + 1), alpha_k_m1) - 1.0) / + alpha_k_m1; + + host_townsend_Y_k(i) = host_townsend_Y_k(i + 1) - step; + } + + Kokkos::deep_copy(townsend_alpha_k_, host_townsend_alpha_k); + Kokkos::deep_copy(townsend_Y_k_, host_townsend_Y_k); + } + + // Create a lightweight object for computing cooling rates within kernels + const auto mbar_over_kb = hydro_pkg->Param("mbar_over_kb"); + const auto adiabatic_index = hydro_pkg->Param("AdiabaticIndex"); + const auto He_mass_fraction = hydro_pkg->Param("He_mass_fraction"); + + cooling_table_obj_ = CoolingTableObj(log_lambdas_, log_temp_start_, log_temp_final_, + d_log_temp_, n_temp_, mbar_over_kb, + adiabatic_index, 1.0 - He_mass_fraction, units); +} + +void TabularCooling::SrcTerm(MeshData *md, const Real dt) const { + if (integrator_ == CoolIntegrator::rk12) { + SubcyclingFixedIntSrcTerm(md, dt, RK12Stepper()); + } else if (integrator_ == CoolIntegrator::rk45) { + SubcyclingFixedIntSrcTerm(md, dt, RK45Stepper()); + } else if (integrator_ == CoolIntegrator::townsend) { + TownsendSrcTerm(md, dt); + } else { + PARTHENON_FAIL("Unknown cooling integrator."); + } +} + +template +void TabularCooling::SubcyclingFixedIntSrcTerm(MeshData *md, const Real dt_, + const RKStepper rk_stepper) const { + + const auto dt = dt_; // HACK capturing parameters still broken with Cuda 11.6 ... + auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro"); + const bool mhd_enabled = hydro_pkg->Param("fluid") == Fluid::glmmhd; + // Grab member variables for compiler + + const CoolingTableObj cooling_table_obj = cooling_table_obj_; + const auto gm1 = (hydro_pkg->Param("AdiabaticIndex") - 1.0); + const auto mbar_gm1_over_kb = hydro_pkg->Param("mbar_over_kb") * gm1; + + const unsigned int max_iter = max_iter_; + + const Real min_sub_dt = dt / max_iter; + + const Real d_e_tol = d_e_tol_; + + // Determine the cooling floor, whichever is higher of the cooling table floor + // or fluid solver floor + const auto temp_cool_floor = std::pow(10.0, log_temp_start_); // low end of cool table + const Real temp_floor = (T_floor_ > temp_cool_floor) ? T_floor_ : temp_cool_floor; + + const Real internal_e_floor = temp_floor / mbar_gm1_over_kb; // specific internal en. + + // Grab some necessary variables + const auto &prim_pack = md->PackVariables(std::vector{"prim"}); + const auto &cons_pack = md->PackVariables(std::vector{"cons"}); + // need to include ghost zones as this source is called prior to the other fluxes when + // split + IndexRange ib = md->GetBlockData(0)->GetBoundsI(IndexDomain::entire); + IndexRange jb = md->GetBlockData(0)->GetBoundsJ(IndexDomain::entire); + IndexRange kb = md->GetBlockData(0)->GetBoundsK(IndexDomain::entire); + + par_for( + DEFAULT_LOOP_PATTERN, "TabularCooling::SubcyclingSplitSrcTerm", DevExecSpace(), 0, + cons_pack.GetDim(5) - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i) { + auto &cons = cons_pack(b); + auto &prim = prim_pack(b); + // Need to use `cons` here as prim may still contain state at t_0; + const Real rho = cons(IDN, k, j, i); + // TODO(pgrete) with potentially more EOS, a separate get_pressure (or similar) + // function could be useful. + Real internal_e = + cons(IEN, k, j, i) - 0.5 * + (SQR(cons(IM1, k, j, i)) + SQR(cons(IM2, k, j, i)) + + SQR(cons(IM3, k, j, i))) / + rho; + if (mhd_enabled) { + internal_e -= 0.5 * (SQR(cons(IB1, k, j, i)) + SQR(cons(IB2, k, j, i)) + + SQR(cons(IB3, k, j, i))); + } + internal_e /= rho; + const Real internal_e_initial = internal_e; + + bool dedt_valid = true; + + // Wrap DeDt into a functor for the RKStepper + auto DeDt_wrapper = [&](const Real t, const Real e, bool &valid) { + return cooling_table_obj.DeDt(e, rho, valid); + }; + + Real sub_t = 0; // current subcycle time + // Try full dt. If error is too large adaptive timestepping will reduce sub_dt + Real sub_dt = dt; + + // Check if cooling is actually happening, e.g., when T below T_cool_min or if + // temperature is already below floor. + const Real dedt_initial = DeDt_wrapper(0.0, internal_e_initial, dedt_valid); + if (dedt_initial == 0.0 || internal_e_initial <= internal_e_floor) { + return; + } + + // Use minumum subcycle timestep when d_e_tol == 0 + if (d_e_tol == 0) { + sub_dt = min_sub_dt; + } + + unsigned int sub_iter = 0; + // check for dedt != 0.0 required in case cooling floor it hit during subcycling + while ((sub_t * (1 + KEpsilon_) < dt) && + (DeDt_wrapper(sub_t, internal_e, dedt_valid) != 0.0)) { + + if (sub_iter > max_iter) { + // Due to sub_dt >= min_dt, this error should never happen + PARTHENON_FAIL( + "FATAL ERROR in [TabularCooling::SubcyclingFixedIntSrcTerm]: Sub " + "cycles exceed max_iter (This should be impossible)"); + } + + // Next higher order estimate + Real internal_e_next_h; + // Error in estimate of higher order + Real d_e_err; + // Number of attempts on this subcycle + unsigned int sub_attempt = 0; + // Whether to reattempt this subcycle + bool reattempt_sub = true; + do { + // Next lower order estimate + Real internal_e_next_l; + // Do one dual order RK step + dedt_valid = true; + RKStepper::Step(sub_t, sub_dt, internal_e, DeDt_wrapper, internal_e_next_h, + internal_e_next_l, dedt_valid); + + sub_attempt++; + + if (!dedt_valid) { + if (sub_dt == min_sub_dt) { + // Cooling is so fast that even the minimum subcycle dt would lead to + // negative internal energy -- so just cool to the floor of the cooling + // table + sub_dt = (dt - sub_t); + internal_e_next_h = internal_e_floor; + reattempt_sub = false; + } else { + reattempt_sub = true; + sub_dt = min_sub_dt; + } + } else { + + // Compute error + d_e_err = fabs((internal_e_next_h - internal_e_next_l) / internal_e_next_h); + + reattempt_sub = false; + // Accepting or reattempting the subcycle: + // + // -If the error is small, accept the subcycle + // + // -If the error on the subcycle is too high, compute a new time + // step to reattempt the subcycle + // -But if the new time step is smaller than the minimum subcycle + // time step (total step duration/ max iterations), just use the + // minimum subcycle time step instead + + if (std::isnan(d_e_err)) { + reattempt_sub = true; + sub_dt = min_sub_dt; + } else if (d_e_err >= d_e_tol && sub_dt > min_sub_dt) { + // Reattempt this subcycle + reattempt_sub = true; + // Error was too high, shrink the timestep + if (d_e_tol == 0) { + sub_dt = min_sub_dt; + } else { + sub_dt = RKStepper::OptimalStep(sub_dt, d_e_err, d_e_tol); + } + // Don't drop timestep under maximum iteration count + if (sub_dt < min_sub_dt || sub_attempt >= max_iter) { + sub_dt = min_sub_dt; + } + } + } + + } while (reattempt_sub); + // Accept this subcycle + sub_t += sub_dt; + + internal_e = internal_e_next_h; + + // skip to the end of subcycling if error is 0 (very unlikely) + if (d_e_err == 0) { + sub_dt = dt - sub_t; + } else { + // Grow the timestep + // (or shrink in case d_e_err >= d_e_tol and sub_dt is already at min_sub_dt) + sub_dt = RKStepper::OptimalStep(sub_dt, d_e_err, d_e_tol); + } + + if (d_e_tol == 0) { + sub_dt = min_sub_dt; + } + + // Don't drop timestep under the minimum step size + sub_dt = std::max(sub_dt, min_sub_dt); + + // Limit by end time + sub_dt = std::min(sub_dt, dt - sub_t); + + sub_iter++; + } + + // If cooled below floor, reset to floor value. + // This could happen if the floor value is larger than the lower end of the + // cooling table or if they are close and the last subcycle in the cooling above + // the lower end pushed the temperature below the lower end (and the floor). + internal_e = (internal_e > internal_e_floor) ? internal_e : internal_e_floor; + + // Remove the cooling from the total energy density + cons(IEN, k, j, i) += rho * (internal_e - internal_e_initial); + // Latter technically not required if no other tasks follows before + // ConservedToPrim conversion, but keeping it for now (better safe than sorry). + prim(IPR, k, j, i) = rho * internal_e * gm1; + }); +} + +void TabularCooling::TownsendSrcTerm(parthenon::MeshData *md, + const parthenon::Real dt_) const { + auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro"); + const bool mhd_enabled = hydro_pkg->Param("fluid") == Fluid::glmmhd; + + // Grab member variables for compiler + const auto dt = dt_; // HACK capturing parameters still broken with Cuda 11.6 ... + + const auto units = hydro_pkg->Param("units"); + const auto gm1 = (hydro_pkg->Param("AdiabaticIndex") - 1.0); + const auto mbar_gm1_over_kb = hydro_pkg->Param("mbar_over_kb") * gm1; + const Real X_by_mh2 = + std::pow((1 - hydro_pkg->Param("He_mass_fraction")) / units.mh(), 2); + + const auto lambdas = lambdas_; + const auto temps = temps_; + const auto alpha_k = townsend_alpha_k_; + const auto Y_k = townsend_Y_k_; + + const auto internal_e_floor = T_floor_ / mbar_gm1_over_kb; + const auto temp_cool_floor = std::pow(10.0, log_temp_start_); // low end of cool table + + // Grab some necessary variables + const auto &prim_pack = md->PackVariables(std::vector{"prim"}); + const auto &cons_pack = md->PackVariables(std::vector{"cons"}); + // need to include ghost zones as this source is called prior to the other fluxes when + // split + IndexRange ib = md->GetBlockData(0)->GetBoundsI(IndexDomain::entire); + IndexRange jb = md->GetBlockData(0)->GetBoundsJ(IndexDomain::entire); + IndexRange kb = md->GetBlockData(0)->GetBoundsK(IndexDomain::entire); + + const auto nbins = alpha_k.extent_int(0); + + // Get reference values + const auto temp_final = std::pow(10.0, log_temp_final_); + const auto lambda_final = lambda_final_; + + par_for( + DEFAULT_LOOP_PATTERN, "TabularCooling::TownsendSrcTerm", DevExecSpace(), 0, + cons_pack.GetDim(5) - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i) { + auto &cons = cons_pack(b); + auto &prim = prim_pack(b); + // Need to use `cons` here as prim may still contain state at t_0; + const auto rho = cons(IDN, k, j, i); + // TODO(pgrete) with potentially more EOS, a separate get_pressure (or similar) + // function could be useful. + auto internal_e = + cons(IEN, k, j, i) - 0.5 * + (SQR(cons(IM1, k, j, i)) + SQR(cons(IM2, k, j, i)) + + SQR(cons(IM3, k, j, i))) / + rho; + if (mhd_enabled) { + internal_e -= 0.5 * (SQR(cons(IB1, k, j, i)) + SQR(cons(IB2, k, j, i)) + + SQR(cons(IB3, k, j, i))); + } + internal_e /= rho; + + // If temp is below floor, reset and return + if (internal_e <= internal_e_floor) { + // Remove the cooling from the total energy density + cons(IEN, k, j, i) += rho * (internal_e_floor - internal_e); + // Latter technically not required if no other tasks follows before + // ConservedToPrim conversion, but keeping it for now (better safe than sorry). + prim(IPR, k, j, i) = rho * internal_e_floor * gm1; + return; + } + + auto temp = mbar_gm1_over_kb * internal_e; + // Temperature is above floor (see conditional above) but below cooling table: + // -> no cooling + if (temp < temp_cool_floor) { + return; + } + const Real n_h2_by_rho = rho * X_by_mh2; + + // Get the index of the right temperature bin + // TODO(?) this could be optimized for using a binary search + auto idx = 0; + while ((idx < nbins - 1) && (temps(idx + 1) < temp)) { + idx += 1; + } + + // Compute the Temporal Evolution Function Y(T) (Eq. A5) + const auto alpha_k_m1 = alpha_k(idx) - 1.0; + const auto tef = + Y_k(idx) + (lambda_final / lambdas(idx)) * (temps(idx) / temp_final) * + (std::pow(temps(idx) / temp, alpha_k_m1) - 1.0) / alpha_k_m1; + + // Compute the adjusted TEF for new timestep (Eqn. 26) (term in brackets) + const auto tef_adj = + tef + lambda_final * dt / temp_final * mbar_gm1_over_kb * n_h2_by_rho; + + // TEF is a strictly decreasing function and new_tef > tef + // Check if the new TEF falls into a lower bin, i.e., find the right bin for A7 + // If so, update slopes and coefficients + while ((idx > 0) && (tef_adj > Y_k(idx))) { + idx -= 1; + } + + // Compute the Inverse Temporal Evolution Function Y^{-1}(Y) (Eq. A7) + const auto temp_new = + temps(idx) * + std::pow(1 - (1.0 - alpha_k(idx)) * (lambdas(idx) / lambda_final) * + (temp_final / temps(idx)) * (tef_adj - Y_k(idx)), + 1.0 / (1.0 - alpha_k(idx))); + // Set new temp (at the lowest to the lower end of the cooling table) + const auto internal_e_new = temp_new > temp_cool_floor + ? temp_new / mbar_gm1_over_kb + : temp_cool_floor / mbar_gm1_over_kb; + cons(IEN, k, j, i) += rho * (internal_e_new - internal_e); + // Latter technically not required if no other tasks follows before + // ConservedToPrim conversion, but keeping it for now (better safe than sorry). + prim(IPR, k, j, i) = rho * internal_e_new * gm1; + }); +} + +Real TabularCooling::EstimateTimeStep(MeshData *md) const { + if (cooling_time_cfl_ <= 0.0) { + return std::numeric_limits::max(); + } + + if (isnan(cooling_time_cfl_) || isinf(cooling_time_cfl_)) { + return std::numeric_limits::infinity(); + } + + auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro"); + + const CoolingTableObj cooling_table_obj = cooling_table_obj_; + const auto gm1 = (hydro_pkg->Param("AdiabaticIndex") - 1.0); + const auto mbar_gm1_over_kb = hydro_pkg->Param("mbar_over_kb") * gm1; + + // Determine the cooling floor, whichever is higher of the cooling table floor + // or fluid solver floor + const auto temp_cool_floor = std::pow(10.0, log_temp_start_); // low end of cool table + const Real temp_floor = (T_floor_ > temp_cool_floor) ? T_floor_ : temp_cool_floor; + + const Real internal_e_floor = temp_floor / mbar_gm1_over_kb; // specific internal en. + + // Grab some necessary variables + const auto &prim_pack = md->PackVariables(std::vector{"prim"}); + IndexRange ib = md->GetBlockData(0)->GetBoundsI(IndexDomain::interior); + IndexRange jb = md->GetBlockData(0)->GetBoundsJ(IndexDomain::interior); + IndexRange kb = md->GetBlockData(0)->GetBoundsK(IndexDomain::interior); + + Real min_cooling_time = std::numeric_limits::infinity(); + Kokkos::Min reducer_min(min_cooling_time); + + Kokkos::parallel_reduce( + "TabularCooling::TimeStep", + Kokkos::MDRangePolicy>( + {0, kb.s, jb.s, ib.s}, {prim_pack.GetDim(5), kb.e + 1, jb.e + 1, ib.e + 1}, + {1, 1, 1, ib.e + 1 - ib.s}), + KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i, + Real &thread_min_cooling_time) { + auto &prim = prim_pack(b); + + const Real rho = prim(IDN, k, j, i); + const Real pres = prim(IPR, k, j, i); + + const Real internal_e = pres / (rho * gm1); + + const Real de_dt = cooling_table_obj.DeDt(internal_e, rho); + + // Compute cooling time + // If de_dt is zero (temperature is smaller than lower end of cooling table) or + // current temp is below floor, use infinite cooling time + const Real cooling_time = ((de_dt == 0) || (internal_e < internal_e_floor)) + ? std::numeric_limits::infinity() + : fabs(internal_e / de_dt); + + thread_min_cooling_time = std::min(cooling_time, thread_min_cooling_time); + }, + reducer_min); + + return cooling_time_cfl_ * min_cooling_time; +} + +void TabularCooling::TestCoolingTable(ParameterInput *pin) const { + + const std::string test_filename = pin->GetString("cooling", "test_filename"); + + const auto rho0 = pin->GetReal("cooling", "test_rho0"); + const auto rho1 = pin->GetReal("cooling", "test_rho1"); + const auto n_rho = pin->GetInteger("cooling", "test_n_rho"); + + const auto pres0 = pin->GetReal("cooling", "test_pres0"); + const auto pres1 = pin->GetReal("cooling", "test_pres1"); + const auto n_pres = pin->GetInteger("cooling", "test_n_pres"); + + // Grab member variables for compiler + + // Everything needed by DeDt + const CoolingTableObj cooling_table_obj = cooling_table_obj_; + const auto gm1 = pin->GetReal("hydro", "gamma") - 1.0; + + // Make some device arrays to store the test data + ParArray2D d_rho("d_rho", n_rho, n_pres), d_pres("d_pres", n_rho, n_pres), + d_internal_e("d_internal_e", n_rho, n_pres), d_de_dt("d_de_dt", n_rho, n_pres); + + par_for( + loop_pattern_mdrange_tag, "TabularCooling::TestCoolingTable", DevExecSpace(), 0, + n_rho - 1, 0, n_pres - 1, KOKKOS_LAMBDA(const int &j, const int &i) { + const Real rho = rho0 * pow(rho1 / rho0, static_cast(j) / (n_rho - 1)); + const Real pres = pres0 * pow(pres1 / pres0, static_cast(i) / (n_pres - 1)); + + d_rho(j, i) = rho; + d_pres(j, i) = pres; + + const Real internal_e = pres / (rho * gm1); + + d_internal_e(j, i) = internal_e; + + const Real de_dt = cooling_table_obj.DeDt(internal_e, rho); + + d_de_dt(j, i) = de_dt; + }); + + // Copy Device arrays to host + auto h_rho = Kokkos::create_mirror_view_and_copy(HostMemSpace(), d_rho); + auto h_pres = Kokkos::create_mirror_view_and_copy(HostMemSpace(), d_pres); + auto h_internal_e = Kokkos::create_mirror_view_and_copy(HostMemSpace(), d_internal_e); + auto h_de_dt = Kokkos::create_mirror_view_and_copy(HostMemSpace(), d_de_dt); + + // Write to file + std::ofstream file(test_filename); + file << "#rho pres internal_e de_dt" << std::endl; + for (int j = 0; j < n_rho; j++) { + for (int i = 0; i < n_pres; i++) { + file << h_rho(j, i) << " " << h_pres(j, i) << " " << h_internal_e(j, i) << " " + << h_de_dt(j, i) << " " << std::endl; + } + } +} + +} // namespace cooling diff --git a/src/pgen/blast.cpp b/src/pgen/blast.cpp index 659af95f..c4aa0919 100644 --- a/src/pgen/blast.cpp +++ b/src/pgen/blast.cpp @@ -131,16 +131,16 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { Real drat = pin->GetOrAddReal("problem/blast", "density_ratio", 1.0); Real gamma = pin->GetOrAddReal("hydro", "gamma", 5 / 3); Real gm1 = gamma - 1.0; - + // get coordinates of center of blast, and convert to Cartesian if necessary Real x0 = pin->GetOrAddReal("problem/blast", "x1_0", 0.0); Real y0 = pin->GetOrAddReal("problem/blast", "x2_0", 0.0); Real z0 = pin->GetOrAddReal("problem/blast", "x3_0", 0.0); - + IndexRange ib = pmb->cellbounds.GetBoundsI(IndexDomain::interior); IndexRange jb = pmb->cellbounds.GetBoundsJ(IndexDomain::interior); IndexRange kb = pmb->cellbounds.GetBoundsK(IndexDomain::interior); - + // initialize conserved variables auto &rc = pmb->meshblock_data.Get(); auto &u_dev = rc->Get("cons").data; diff --git a/src/pgen/cluster.cpp b/src/pgen/cluster.cpp index 1e21f534..2882c96d 100644 --- a/src/pgen/cluster.cpp +++ b/src/pgen/cluster.cpp @@ -23,6 +23,8 @@ #include // stringstream #include // runtime_error #include // c_str() +// #include // HDF5 +#include "../../external/HighFive/include/highfive/H5Easy.hpp" // Parthenon headers #include "kokkos_abstraction.hpp" @@ -100,7 +102,12 @@ void ApplyClusterClips(MeshData *md, const parthenon::SimTime &tm, if (r2 < clip_r2) { // Cell falls within clipping radius - eos.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i); + + const int gks = 0; + const int gjs = 0; + const int gis = 0; + + eos.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i, gks, gjs, gis); // Three last parameters are just passive here if (dfloor > 0) { const Real rho = prim(IDN, k, j, i); @@ -252,7 +259,7 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *hyd const Real uniform_b_field_bx = pin->GetReal("problem/cluster/uniform_b_field", "bx"); const Real uniform_b_field_by = pin->GetReal("problem/cluster/uniform_b_field", "by"); const Real uniform_b_field_bz = pin->GetReal("problem/cluster/uniform_b_field", "bz"); - + hydro_pkg->AddParam<>("uniform_b_field_bx", uniform_b_field_bx); hydro_pkg->AddParam<>("uniform_b_field_by", uniform_b_field_by); hydro_pkg->AddParam<>("uniform_b_field_bz", uniform_b_field_bz); @@ -276,6 +283,7 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *hyd hydro_pkg->AddParam<>("dipole_b_field_mz", dipole_b_field_mz); } + /************************************************************ * Read Cluster Gravity Parameters ************************************************************/ @@ -403,11 +411,25 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *hyd // plasma beta hydro_pkg->AddField("plasma_beta", m); } - + + + + + + + /************************************************************ - * Add infrastructure for initial pertubations + * Read Density perturbation ************************************************************/ + + //const bool init_perturb_rho = pin->GetOrAddBoolean("problem/cluster/init_perturb", "init_perturb_rho", false); + //hydro_pkg->AddParam<>("init_perturb_rho", init_perturb_rho); + + /************************************************************ + * Read Velocity perturbation + ************************************************************/ + const auto sigma_v = pin->GetOrAddReal("problem/cluster/init_perturb", "sigma_v", 0.0); if (sigma_v != 0.0) { // peak of init vel perturb @@ -506,12 +528,14 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md) { const Real uy = hydro_pkg->Param("uniform_gas_uy"); const Real uz = hydro_pkg->Param("uniform_gas_uz"); const Real pres = hydro_pkg->Param("uniform_gas_pres"); - + const Real Mx = rho * ux; const Real My = rho * uy; const Real Mz = rho * uz; const Real E = rho * (0.5 * (ux * ux + uy * uy + uz * uz) + pres / (gm1 * rho)); - + + std::cout << "Setting uniform gas"; + parthenon::par_for( DEFAULT_LOOP_PATTERN, "Cluster::ProblemGenerator::UniformGas", parthenon::DevExecSpace(), kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, @@ -522,7 +546,7 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md) { u(IM3, k, j, i) = Mz; u(IEN, k, j, i) = E; }); - + // end if(init_uniform_gas) } else { /************************************************************ @@ -532,9 +556,9 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md) { hydro_pkg ->Param>( "hydrostatic_equilibirum_sphere"); - + const auto P_rho_profile = he_sphere.generate_P_rho_profile(ib, jb, kb, coords); - + // initialize conserved variables parthenon::par_for( DEFAULT_LOOP_PATTERN, "cluster::ProblemGenerator::UniformGas", @@ -664,27 +688,136 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md) { } // END if(hydro_pkg->Param("fluid") == Fluid::glmmhd) } - /************************************************************ - * Set initial velocity perturbations (requires no other velocities for now) - ************************************************************/ - + /************************************************************ + * Initial parameters + ************************************************************/ + auto pmb = md->GetBlockData(0)->GetBlockPointer(); IndexRange ib = pmb->cellbounds.GetBoundsI(IndexDomain::interior); IndexRange jb = pmb->cellbounds.GetBoundsJ(IndexDomain::interior); IndexRange kb = pmb->cellbounds.GetBoundsK(IndexDomain::interior); + auto hydro_pkg = pmb->packages.Get("Hydro"); const auto fluid = hydro_pkg->Param("fluid"); auto const &cons = md->PackVariables(std::vector{"cons"}); - const auto num_blocks = md->NumBlocks(); - + const auto num_blocks = md->NumBlocks(); + + /************************************************************ + * Set initial density perturbations + ************************************************************/ + + const bool init_perturb_rho = pin->GetOrAddBoolean("problem/cluster/init_perturb", "init_perturb_rho", false); + const bool overpressure_ring = pin->GetOrAddBoolean("problem/cluster/init_perturb", "overpressure_ring", false); + const bool spherical_collapse = pin->GetOrAddBoolean("problem/cluster/init_perturb", "spherical_collapse", false); + + hydro_pkg->AddParam<>("init_perturb_rho", init_perturb_rho); + + if (init_perturb_rho == true) { + + auto init_perturb_rho_file = pin->GetString("problem/cluster/init_perturb", "init_perturb_rho_file"); + auto init_perturb_rho_keys = pin->GetString("problem/cluster/init_perturb", "init_perturb_rho_keys"); + + hydro_pkg->AddParam<>("cluster/init_perturb_rho_file", init_perturb_rho_file); + hydro_pkg->AddParam<>("cluster/init_perturb_rho_keys", init_perturb_rho_keys); + + std::cout << "Setting density perturbation"; + + // Read HDF5 file containing the density + std::string filename_rho = "/work/bbd0833/test/rho.h5"; + std::string keys_rho = "data"; + H5Easy::File file(filename_rho, HighFive::File::ReadOnly); + auto rho_init = H5Easy::load, 64>, 64>>(file, keys_rho); + + Real passive_scalar = 0.0; // Useless + + std::cout << "entering initialisation of rho field"; + + pmb->par_reduce( + "Init density field", 0, num_blocks - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int b, const int k, const int j, const int i, Real &lsum) { + + auto pmbb = md->GetBlockData(b)->GetBlockPointer(); // Meshblock b + + const auto gis = pmbb->loc.lx1 * pmb->block_size.nx1; + const auto gjs = pmbb->loc.lx2 * pmb->block_size.nx2; + const auto gks = pmbb->loc.lx3 * pmb->block_size.nx3; + + //std::cout << "gis=" << gis << "\n"; + + const auto &coords = cons.GetCoords(b); + const auto &u = cons(b); + + // Adding the value of the density field + + u(IDN, k, j, i) = rho_init[gks + k - 2][gjs + j - 2][gis + i - 2]; + + // Ring of overpressure + // compute radius + if (overpressure_ring){ + + const Real x = coords.Xc<1>(i); + const Real y = coords.Xc<2>(j); + const Real z = coords.Xc<3>(k); + const Real r = std::sqrt(SQR(x) + SQR(y) + SQR(z)); // Computing radius + const Real width = 0.002; + const Real overpressure_radius = pin->GetOrAddReal("problem/cluster/init_perturb", "overpressure_radius", 0.01); + + if (r > 0 && r < overpressure_radius){ + + //u(IEN, k, j, i) *= 100; + + // Pushing a ring of gas outward (kinetic boost) + const Real u_x = x / r; + const Real u_y = y / r; + const Real u_z = z / r; + + const Real velocity_blast = pin->GetOrAddReal("problem/cluster/init_perturb", "velocity_blast", 0.0); + + //u(IEN, k, j, i) *= 100; + + u(IM1, k, j, i) = u_x * velocity_blast * u(IDN, k, j, i); // p = (1 Mpc / Gyr) * rho + u(IM2, k, j, i) = u_y * velocity_blast * u(IDN, k, j, i); // p = (1 Mpc / Gyr) * rho + u(IM3, k, j, i) = u_z * velocity_blast * u(IDN, k, j, i); // p = (1 Mpc / Gyr) * rho + u(IEN, k, j, i) += 0.5 * u(IDN, k, j, i) * (pow(velocity_blast,2)); + } + } + + if (spherical_collapse){ + + const Real x = coords.Xc<1>(i); + const Real y = coords.Xc<2>(j); + const Real z = coords.Xc<3>(k); + const Real r = std::sqrt(SQR(x) + SQR(y) + SQR(z)); // Computing radius + const Real overpressure_radius = pin->GetOrAddReal("problem/cluster/init_perturb", "overpressure_radius", 0.01); + + u(IDN, k, j, i) = 3; + + if (r > 0 && r < overpressure_radius){ + + u(IDN, k, j, i) *= 50; + + } + + } + + }, + passive_scalar); + + } + + /************************************************************ + * Set initial velocity perturbations (requires no other velocities for now) + ************************************************************/ + const auto sigma_v = pin->GetOrAddReal("problem/cluster/init_perturb", "sigma_v", 0.0); - + if (sigma_v != 0.0) { auto few_modes_ft = hydro_pkg->Param("cluster/few_modes_ft_v"); // Init phases on all blocks for (int b = 0; b < md->NumBlocks(); b++) { auto pmb = md->GetBlockData(b)->GetBlockPointer(); few_modes_ft.SetPhases(pmb.get(), pin); + } // As for t_corr in few_modes_ft, the choice for dt is // in principle arbitrary because the inital v_hat is 0 and the v_hat_new will contain @@ -877,14 +1010,15 @@ void UserWorkBeforeOutput(MeshBlock *pmb, ParameterInput *pin) { // compute temperature temperature(k, j, i) = mbar_over_kb * P / rho; }); + if (pkg->Param("enable_cooling") == Cooling::tabular) { auto &cooling_time = data->Get("cooling_time").data; - + // get cooling function const cooling::TabularCooling &tabular_cooling = pkg->Param("tabular_cooling"); const auto cooling_table_obj = tabular_cooling.GetCoolingTableObj(); - + pmb->par_for( "Cluster::UserWorkBeforeOutput::CoolingTime", kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, KOKKOS_LAMBDA(const int k, const int j, const int i) { diff --git a/src/pgen/cluster/agn_feedback.cpp b/src/pgen/cluster/agn_feedback.cpp index 1c76e913..f12af3e9 100644 --- a/src/pgen/cluster/agn_feedback.cpp +++ b/src/pgen/cluster/agn_feedback.cpp @@ -338,7 +338,7 @@ void AGNFeedback::FeedbackSrcTerm(parthenon::MeshData *md, // momentum, and total energy added depend on the triggered power. /////////////////////////////////////////////////////////////////// - eos.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i); + eos.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i,0,0,0); const Real old_specific_internal_e = prim(IPR, k, j, i) / (prim(IDN, k, j, i) * (eos.GetGamma() - 1.)); @@ -348,7 +348,7 @@ void AGNFeedback::FeedbackSrcTerm(parthenon::MeshData *md, cons(IM3, k, j, i) += jet_momentum * sign_jet * jet_axis_z; cons(IEN, k, j, i) += jet_feedback; - eos.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i); + eos.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i,0,0,0); const Real new_specific_internal_e = prim(IPR, k, j, i) / (prim(IDN, k, j, i) * (eos.GetGamma() - 1.)); PARTHENON_REQUIRE( @@ -382,7 +382,7 @@ void AGNFeedback::FeedbackSrcTerm(parthenon::MeshData *md, } } - eos.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i); + eos.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i,0,0,0); PARTHENON_REQUIRE(prim(IPR, k, j, i) > 0, "Kinetic injection leads to negative pressure"); }); diff --git a/src/pgen/cluster/agn_triggering.cpp b/src/pgen/cluster/agn_triggering.cpp index 242de4fa..5f39957e 100644 --- a/src/pgen/cluster/agn_triggering.cpp +++ b/src/pgen/cluster/agn_triggering.cpp @@ -188,7 +188,7 @@ void AGNTriggering::ReduceColdMass(parthenon::Real &cold_mass, AddDensityToConsAtFixedVelTemp(cell_delta_rho, cons, prim, eos.GetGamma(), k, j, i); // Update the Primitives - eos.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i); + eos.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i,0,0,0); } } } @@ -309,7 +309,7 @@ void AGNTriggering::RemoveBondiAccretedGas(parthenon::MeshData i); // Update the Primitives - eos.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i); + eos.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i,0,0,0); } }); } diff --git a/src/pgen/cluster/snia_feedback.cpp b/src/pgen/cluster/snia_feedback.cpp index ed41657f..25087235 100644 --- a/src/pgen/cluster/snia_feedback.cpp +++ b/src/pgen/cluster/snia_feedback.cpp @@ -114,7 +114,7 @@ void SNIAFeedback::FeedbackSrcTerm(parthenon::MeshData *md, cons(IEN, k, j, i) += snia_energy_density; AddDensityToConsAtFixedVel(snia_mass_density, cons, prim, k, j, i); - eos.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i); + eos.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i,0,0,0); }); } diff --git a/src/pgen/cluster_modified.cpp b/src/pgen/cluster_modified.cpp new file mode 100644 index 00000000..623a7e88 --- /dev/null +++ b/src/pgen/cluster_modified.cpp @@ -0,0 +1,1023 @@ +//======================================================================================== +// AthenaPK - a performance portable block structured AMR astrophysical MHD code. +// Copyright (c) 2021-2023, Athena-Parthenon Collaboration. All rights reserved. +// Licensed under the 3-clause BSD License, see LICENSE file for details +//======================================================================================== +//! \file cluster.cpp +// \brief Idealized galaxy cluster problem generator +// +// Setups up an idealized galaxy cluster with an ACCEPT-like entropy profile in +// hydrostatic equilbrium with an NFW+BCG+SMBH gravitational profile, +// optionally with an initial magnetic tower field. Includes AGN feedback, AGN +// triggering via cold gas, simple SNIA Feedback(TODO) +//======================================================================================== + +// C headers + +// C++ headers +#include // min, max +#include // sqrt() +#include // fopen(), fprintf(), freopen() +#include // endl +#include +#include // stringstream +#include // runtime_error +#include // c_str() +// #include // HDF5 +#include "../../external/HighFive/include/highfive/H5Easy.hpp" + + +// Parthenon headers +#include "kokkos_abstraction.hpp" +#include "mesh/domain.hpp" +#include "mesh/mesh.hpp" +#include +#include + +// AthenaPK headers +#include "../eos/adiabatic_glmmhd.hpp" +#include "../eos/adiabatic_hydro.hpp" +#include "../hydro/hydro.hpp" +#include "../hydro/srcterms/gravitational_field.hpp" +#include "../hydro/srcterms/tabular_cooling.hpp" +#include "../main.hpp" +#include "../utils/few_modes_ft.hpp" + +// Cluster headers +#include "cluster/agn_feedback.hpp" +#include "cluster/agn_triggering.hpp" +#include "cluster/cluster_gravity.hpp" +#include "cluster/entropy_profiles.hpp" +#include "cluster/hydrostatic_equilibrium_sphere.hpp" +#include "cluster/magnetic_tower.hpp" +#include "cluster/snia_feedback.hpp" +#include "parthenon_array_generic.hpp" +#include "utils/error_checking.hpp" + +namespace cluster { +using namespace parthenon::driver::prelude; +using namespace parthenon::package::prelude; +using utils::few_modes_ft::FewModesFT; + +template +void ApplyClusterClips(MeshData *md, const parthenon::SimTime &tm, + const Real beta_dt, const EOS eos) { + + auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro"); + + // Apply clips -- ceilings on temperature, velocity, alfven velocity, and + // density floor -- within a radius of the AGN + const auto &dfloor = hydro_pkg->Param("cluster_dfloor"); + const auto &eceil = hydro_pkg->Param("cluster_eceil"); + const auto &vceil = hydro_pkg->Param("cluster_vceil"); + const auto &vAceil = hydro_pkg->Param("cluster_vAceil"); + const auto &clip_r = hydro_pkg->Param("cluster_clip_r"); + + if (clip_r > 0 && (dfloor > 0 || eceil < std::numeric_limits::infinity() || + vceil < std::numeric_limits::infinity() || + vAceil < std::numeric_limits::infinity())) { + // Grab some necessary variables + const auto &prim_pack = md->PackVariables(std::vector{"prim"}); + const auto &cons_pack = md->PackVariables(std::vector{"cons"}); + IndexRange ib = md->GetBlockData(0)->GetBoundsI(IndexDomain::interior); + IndexRange jb = md->GetBlockData(0)->GetBoundsJ(IndexDomain::interior); + IndexRange kb = md->GetBlockData(0)->GetBoundsK(IndexDomain::interior); + const auto nhydro = hydro_pkg->Param("nhydro"); + const auto nscalars = hydro_pkg->Param("nscalars"); + + const Real clip_r2 = SQR(clip_r); + const Real vceil2 = SQR(vceil); + const Real vAceil2 = SQR(vAceil); + const Real gm1 = (hydro_pkg->Param("AdiabaticIndex") - 1.0); + + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "Cluster::ApplyClusterClips", parthenon::DevExecSpace(), 0, + cons_pack.GetDim(5) - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i) { + auto &cons = cons_pack(b); + auto &prim = prim_pack(b); + const auto &coords = cons_pack.GetCoords(b); + + const Real r2 = + SQR(coords.Xc<1>(i)) + SQR(coords.Xc<2>(j)) + SQR(coords.Xc<3>(k)); + + if (r2 < clip_r2) { + // Cell falls within clipping radius + eos.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i); + + if (dfloor > 0) { + const Real rho = prim(IDN, k, j, i); + if (rho < dfloor) { + cons(IDN, k, j, i) = dfloor; + prim(IDN, k, j, i) = dfloor; + } + } + + if (vceil < std::numeric_limits::infinity()) { + // Apply velocity ceiling + const Real v2 = SQR(prim(IV1, k, j, i)) + SQR(prim(IV2, k, j, i)) + + SQR(prim(IV3, k, j, i)); + if (v2 > vceil2) { + // Fix the velocity to the velocity ceiling + const Real v = sqrt(v2); + cons(IM1, k, j, i) *= vceil / v; + cons(IM2, k, j, i) *= vceil / v; + cons(IM3, k, j, i) *= vceil / v; + prim(IV1, k, j, i) *= vceil / v; + prim(IV2, k, j, i) *= vceil / v; + prim(IV3, k, j, i) *= vceil / v; + + // Remove kinetic energy + cons(IEN, k, j, i) -= 0.5 * prim(IDN, k, j, i) * (v2 - vceil2); + } + } + + if (vAceil2 < std::numeric_limits::infinity()) { + // Apply Alfven velocity ceiling by raising density + const Real rho = prim(IDN, k, j, i); + const Real B2 = (SQR(prim(IB1, k, j, i)) + SQR(prim(IB2, k, j, i)) + + SQR(prim(IB3, k, j, i))); + + // compute Alfven mach number + const Real va2 = (B2 / rho); + + if (va2 > vAceil2) { + // Increase the density to match the alfven velocity ceiling + const Real rho_new = std::sqrt(B2 / vAceil2); + cons(IDN, k, j, i) = rho_new; + prim(IDN, k, j, i) = rho_new; + } + } + + if (eceil < std::numeric_limits::infinity()) { + // Apply internal energy ceiling as a pressure ceiling + const Real internal_e = prim(IPR, k, j, i) / (gm1 * prim(IDN, k, j, i)); + if (internal_e > eceil) { + cons(IEN, k, j, i) -= prim(IDN, k, j, i) * (internal_e - eceil); + prim(IPR, k, j, i) = gm1 * prim(IDN, k, j, i) * eceil; + } + } + } + }); + } +} + +void ApplyClusterClips(MeshData *md, const parthenon::SimTime &tm, + const Real beta_dt) { + auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro"); + auto fluid = hydro_pkg->Param("fluid"); + if (fluid == Fluid::euler) { + ApplyClusterClips(md, tm, beta_dt, hydro_pkg->Param("eos")); + } else if (fluid == Fluid::glmmhd) { + ApplyClusterClips(md, tm, beta_dt, hydro_pkg->Param("eos")); + } else { + PARTHENON_FAIL("Cluster::ApplyClusterClips: Unknown EOS"); + } +} + +void ClusterSrcTerm(MeshData *md, const parthenon::SimTime &tm, + const Real beta_dt) { + auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro"); + + const bool &gravity_srcterm = hydro_pkg->Param("gravity_srcterm"); + + if (gravity_srcterm) { + const ClusterGravity &cluster_gravity = + hydro_pkg->Param("cluster_gravity"); + + GravitationalFieldSrcTerm(md, beta_dt, cluster_gravity); + } + + const auto &agn_feedback = hydro_pkg->Param("agn_feedback"); + agn_feedback.FeedbackSrcTerm(md, beta_dt, tm); + + const auto &magnetic_tower = hydro_pkg->Param("magnetic_tower"); + magnetic_tower.FixedFieldSrcTerm(md, beta_dt, tm); + + const auto &snia_feedback = hydro_pkg->Param("snia_feedback"); + snia_feedback.FeedbackSrcTerm(md, beta_dt, tm); + + ApplyClusterClips(md, tm, beta_dt); +}; + +Real ClusterEstimateTimestep(MeshData *md) { + Real min_dt = std::numeric_limits::max(); + + auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro"); + + // TODO time constraints imposed by thermal AGN feedback, jet velocity, + // magnetic tower + const auto &agn_triggering = hydro_pkg->Param("agn_triggering"); + const Real agn_triggering_min_dt = agn_triggering.EstimateTimeStep(md); + min_dt = std::min(min_dt, agn_triggering_min_dt); + + return min_dt; +} + +//======================================================================================== +//! \fn void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor +//! *hydro_pkg) \brief Init package data from parameter input +//======================================================================================== + +void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *hydro_pkg) { + + /************************************************************ + * Read Uniform Gas + ************************************************************/ + + const bool init_uniform_gas = + pin->GetOrAddBoolean("problem/cluster/uniform_gas", "init_uniform_gas", false); + hydro_pkg->AddParam<>("init_uniform_gas", init_uniform_gas); + + if (init_uniform_gas) { + const Real uniform_gas_rho = pin->GetReal("problem/cluster/uniform_gas", "rho"); + const Real uniform_gas_ux = pin->GetReal("problem/cluster/uniform_gas", "ux"); + const Real uniform_gas_uy = pin->GetReal("problem/cluster/uniform_gas", "uy"); + const Real uniform_gas_uz = pin->GetReal("problem/cluster/uniform_gas", "uz"); + const Real uniform_gas_pres = pin->GetReal("problem/cluster/uniform_gas", "pres"); + + hydro_pkg->AddParam<>("uniform_gas_rho", uniform_gas_rho); + hydro_pkg->AddParam<>("uniform_gas_ux", uniform_gas_ux); + hydro_pkg->AddParam<>("uniform_gas_uy", uniform_gas_uy); + hydro_pkg->AddParam<>("uniform_gas_uz", uniform_gas_uz); + hydro_pkg->AddParam<>("uniform_gas_pres", uniform_gas_pres); + } + + /************************************************************ + * Read Uniform Magnetic Field + ************************************************************/ + + const bool init_uniform_b_field = pin->GetOrAddBoolean( + "problem/cluster/uniform_b_field", "init_uniform_b_field", false); + hydro_pkg->AddParam<>("init_uniform_b_field", init_uniform_b_field); + + if (init_uniform_b_field) { + const Real uniform_b_field_bx = pin->GetReal("problem/cluster/uniform_b_field", "bx"); + const Real uniform_b_field_by = pin->GetReal("problem/cluster/uniform_b_field", "by"); + const Real uniform_b_field_bz = pin->GetReal("problem/cluster/uniform_b_field", "bz"); + + hydro_pkg->AddParam<>("uniform_b_field_bx", uniform_b_field_bx); + hydro_pkg->AddParam<>("uniform_b_field_by", uniform_b_field_by); + hydro_pkg->AddParam<>("uniform_b_field_bz", uniform_b_field_bz); + } + + /************************************************************ + * Read Uniform Magnetic Field + ************************************************************/ + + const bool init_dipole_b_field = pin->GetOrAddBoolean("problem/cluster/dipole_b_field", + "init_dipole_b_field", false); + hydro_pkg->AddParam<>("init_dipole_b_field", init_dipole_b_field); + + if (init_dipole_b_field) { + const Real dipole_b_field_mx = pin->GetReal("problem/cluster/dipole_b_field", "mx"); + const Real dipole_b_field_my = pin->GetReal("problem/cluster/dipole_b_field", "my"); + const Real dipole_b_field_mz = pin->GetReal("problem/cluster/dipole_b_field", "mz"); + + hydro_pkg->AddParam<>("dipole_b_field_mx", dipole_b_field_mx); + hydro_pkg->AddParam<>("dipole_b_field_my", dipole_b_field_my); + hydro_pkg->AddParam<>("dipole_b_field_mz", dipole_b_field_mz); + } + + + /************************************************************ + * Read Cluster Gravity Parameters + ************************************************************/ + + // Build cluster_gravity object + ClusterGravity cluster_gravity(pin, hydro_pkg); + // hydro_pkg->AddParam<>("cluster_gravity", cluster_gravity); + + // Include gravity as a source term during evolution + const bool gravity_srcterm = + pin->GetBoolean("problem/cluster/gravity", "gravity_srcterm"); + hydro_pkg->AddParam<>("gravity_srcterm", gravity_srcterm); + + /************************************************************ + * Read Initial Entropy Profile + ************************************************************/ + + // Build entropy_profile object + ACCEPTEntropyProfile entropy_profile(pin); + + /************************************************************ + * Build Hydrostatic Equilibrium Sphere + ************************************************************/ + + HydrostaticEquilibriumSphere hse_sphere(pin, hydro_pkg, cluster_gravity, + entropy_profile); + + /************************************************************ + * Read Precessing Jet Coordinate system + ************************************************************/ + + JetCoordsFactory jet_coords_factory(pin, hydro_pkg); + + /************************************************************ + * Read AGN Feedback + ************************************************************/ + + AGNFeedback agn_feedback(pin, hydro_pkg); + + /************************************************************ + * Read AGN Triggering + ************************************************************/ + AGNTriggering agn_triggering(pin, hydro_pkg); + + /************************************************************ + * Read Magnetic Tower + ************************************************************/ + + // Build Magnetic Tower + MagneticTower magnetic_tower(pin, hydro_pkg); + + // Determine if magnetic_tower_power_scaling is needed + // Is AGN Power and Magnetic fraction non-zero? + bool magnetic_tower_power_scaling = + (agn_feedback.magnetic_fraction_ != 0 && + (agn_feedback.fixed_power_ != 0 || + agn_triggering.triggering_mode_ != AGNTriggeringMode::NONE)); + hydro_pkg->AddParam("magnetic_tower_power_scaling", magnetic_tower_power_scaling); + + /************************************************************ + * Read SNIA Feedback + ************************************************************/ + + SNIAFeedback snia_feedback(pin, hydro_pkg); + + /************************************************************ + * Read Clips (ceilings and floors) + ************************************************************/ + + // Disable all clips by default with a negative radius clip + Real clip_r = pin->GetOrAddReal("problem/cluster/clips", "clip_r", -1.0); + + // By default disable floors by setting a negative value + Real dfloor = pin->GetOrAddReal("problem/cluster/clips", "dfloor", -1.0); + + // By default disable ceilings by setting to infinity + Real vceil = pin->GetOrAddReal("problem/cluster/clips", "vceil", + std::numeric_limits::infinity()); + Real vAceil = pin->GetOrAddReal("problem/cluster/clips", "vAceil", + std::numeric_limits::infinity()); + Real Tceil = pin->GetOrAddReal("problem/cluster/clips", "Tceil", + std::numeric_limits::infinity()); + Real eceil = Tceil; + if (eceil < std::numeric_limits::infinity()) { + if (!hydro_pkg->AllParams().hasKey("mbar_over_kb")) { + PARTHENON_FAIL("Temperature ceiling requires units and gas composition. " + "Either set a 'units' block and the 'hydro/He_mass_fraction' in " + "input file or use a pressure floor " + "(defined code units) instead."); + } + auto mbar_over_kb = hydro_pkg->Param("mbar_over_kb"); + eceil = Tceil / mbar_over_kb / (hydro_pkg->Param("AdiabaticIndex") - 1.0); + } + hydro_pkg->AddParam("cluster_dfloor", dfloor); + hydro_pkg->AddParam("cluster_eceil", eceil); + hydro_pkg->AddParam("cluster_vceil", vceil); + hydro_pkg->AddParam("cluster_vAceil", vAceil); + hydro_pkg->AddParam("cluster_clip_r", clip_r); + + /************************************************************ + * Add derived fields + * NOTE: these must be filled in UserWorkBeforeOutput + ************************************************************/ + + auto m = Metadata({Metadata::Cell, Metadata::OneCopy}, std::vector({1})); + + // log10 of cell-centered radius + hydro_pkg->AddField("log10_cell_radius", m); + // entropy + hydro_pkg->AddField("entropy", m); + // sonic Mach number v/c_s + hydro_pkg->AddField("mach_sonic", m); + // temperature + hydro_pkg->AddField("temperature", m); + + if (hydro_pkg->Param("enable_cooling") == Cooling::tabular) { + // cooling time + hydro_pkg->AddField("cooling_time", m); + } + + if (hydro_pkg->Param("fluid") == Fluid::glmmhd) { + // alfven Mach number v_A/c_s + hydro_pkg->AddField("mach_alfven", m); + + // plasma beta + hydro_pkg->AddField("plasma_beta", m); + } + + + + + + + + /************************************************************ + * Read Density perturbation + ************************************************************/ + + //const bool init_perturb_rho = pin->GetOrAddBoolean("problem/cluster/init_perturb", "init_perturb_rho", false); + //hydro_pkg->AddParam<>("init_perturb_rho", init_perturb_rho); + + + /************************************************************ + * Read Velocity perturbation + ************************************************************/ + + const auto sigma_v = pin->GetOrAddReal("problem/cluster/init_perturb", "sigma_v", 0.0); + if (sigma_v != 0.0) { + // peak of init vel perturb + auto k_peak_v = pin->GetReal("problem/cluster/init_perturb", "k_peak_v"); + auto num_modes_v = + pin->GetOrAddInteger("problem/cluster/init_perturb", "num_modes_v", 40); + auto sol_weight_v = + pin->GetOrAddReal("problem/cluster/init_perturb", "sol_weight_v", 1.0); + uint32_t rseed_v = pin->GetOrAddInteger("problem/cluster/init_perturb", "rseed_v", 1); + // In principle arbitrary because the inital v_hat is 0 and the v_hat_new will contain + // the perturbation (and is normalized in the following to get the desired sigma_v) + const auto t_corr = 1e-10; + + auto k_vec_v = utils::few_modes_ft::MakeRandomModes(num_modes_v, k_peak_v, rseed_v); + + auto few_modes_ft = FewModesFT(pin, hydro_pkg, "cluster_perturb_v", num_modes_v, + k_vec_v, k_peak_v, sol_weight_v, t_corr, rseed_v); + hydro_pkg->AddParam<>("cluster/few_modes_ft_v", few_modes_ft); + + // Add field for initial perturation (must not need to be consistent but defining it + // this way is easier for now) + Metadata m({Metadata::Cell, Metadata::Derived, Metadata::OneCopy}, + std::vector({3})); + hydro_pkg->AddField("tmp_perturb", m); + } + const auto sigma_b = pin->GetOrAddReal("problem/cluster/init_perturb", "sigma_b", 0.0); + if (sigma_b != 0.0) { + PARTHENON_REQUIRE_THROWS(hydro_pkg->Param("fluid") == Fluid::glmmhd, + "Requested initial magnetic field perturbations but not " + "solving the MHD equations.") + // peak of init vel perturb + auto k_peak_b = pin->GetReal("problem/cluster/init_perturb", "k_peak_b"); + auto num_modes_b = + pin->GetOrAddInteger("problem/cluster/init_perturb", "num_modes_b", 40); + uint32_t rseed_b = pin->GetOrAddInteger("problem/cluster/init_perturb", "rseed_b", 2); + // In principle arbitrary because the inital A_hat is 0 and the A_hat_new will contain + // the perturbation (and is normalized in the following to get the desired sigma_b) + const auto t_corr = 1e-10; + // This field should by construction have no compressive modes, so we fix the number. + const auto sol_weight_b = 1.0; + + auto k_vec_b = utils::few_modes_ft::MakeRandomModes(num_modes_b, k_peak_b, rseed_b); + + const bool fill_ghosts = true; // as we fill a vector potential to calc B + auto few_modes_ft = + FewModesFT(pin, hydro_pkg, "cluster_perturb_b", num_modes_b, k_vec_b, k_peak_b, + sol_weight_b, t_corr, rseed_b, fill_ghosts); + hydro_pkg->AddParam<>("cluster/few_modes_ft_b", few_modes_ft); + + // Add field for initial perturation (must not need to be consistent but defining it + // this way is easier for now). Only add if not already done for the velocity. + if (sigma_v == 0.0) { + Metadata m({Metadata::Cell, Metadata::Derived, Metadata::OneCopy}, + std::vector({3})); + hydro_pkg->AddField("tmp_perturb", m); + } + } +} + +//======================================================================================== +//! \fn void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md) +//! \brief Generate problem data for all blocks on rank +// +// Note, this requires that parthenon/mesh/pack_size=-1 during initialization so that +// reductions work +//======================================================================================== + +void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md) { + // This could be more optimized, but require a refactor of init routines being called. + // However, given that it's just called during initial setup, this should not be a + // performance concern. + for (int b = 0; b < md->NumBlocks(); b++) { + auto pmb = md->GetBlockData(b)->GetBlockPointer(); + auto hydro_pkg = pmb->packages.Get("Hydro"); + + IndexRange ib = pmb->cellbounds.GetBoundsI(IndexDomain::interior); + IndexRange jb = pmb->cellbounds.GetBoundsJ(IndexDomain::interior); + IndexRange kb = pmb->cellbounds.GetBoundsK(IndexDomain::interior); + + // Initialize the conserved variables + auto &u = pmb->meshblock_data.Get()->Get("cons").data; + + auto &coords = pmb->coords; + + // Get Adiabatic Index + const Real gam = pin->GetReal("hydro", "gamma"); + const Real gm1 = (gam - 1.0); + + /************************************************************ + * Initialize the initial hydro state + ************************************************************/ + const auto &init_uniform_gas = hydro_pkg->Param("init_uniform_gas"); + if (init_uniform_gas) { + const Real rho = hydro_pkg->Param("uniform_gas_rho"); + const Real ux = hydro_pkg->Param("uniform_gas_ux"); + const Real uy = hydro_pkg->Param("uniform_gas_uy"); + const Real uz = hydro_pkg->Param("uniform_gas_uz"); + const Real pres = hydro_pkg->Param("uniform_gas_pres"); + + const Real Mx = rho * ux; + const Real My = rho * uy; + const Real Mz = rho * uz; + const Real E = rho * (0.5 * (ux * ux + uy * uy + uz * uz) + pres / (gm1 * rho)); + + std::cout << "Setting uniform gas"; + + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "Cluster::ProblemGenerator::UniformGas", + parthenon::DevExecSpace(), kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int &k, const int &j, const int &i) { + u(IDN, k, j, i) = rho; + u(IM1, k, j, i) = Mx; + u(IM2, k, j, i) = My; + u(IM3, k, j, i) = Mz; + u(IEN, k, j, i) = E; + }); + + // end if(init_uniform_gas) + } else { + /************************************************************ + * Initialize a HydrostaticEquilibriumSphere + ************************************************************/ + const auto &he_sphere = + hydro_pkg + ->Param>( + "hydrostatic_equilibirum_sphere"); + + const auto P_rho_profile = he_sphere.generate_P_rho_profile(ib, jb, kb, coords); + + // initialize conserved variables + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "cluster::ProblemGenerator::UniformGas", + parthenon::DevExecSpace(), kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int &k, const int &j, const int &i) { + // Calculate radius + const Real r = sqrt(coords.Xc<1>(i) * coords.Xc<1>(i) + + coords.Xc<2>(j) * coords.Xc<2>(j) + + coords.Xc<3>(k) * coords.Xc<3>(k)); + + // Get pressure and density from generated profile + const Real P_r = P_rho_profile.P_from_r(r); + const Real rho_r = P_rho_profile.rho_from_r(r); + + // Fill conserved states, 0 initial velocity + u(IDN, k, j, i) = rho_r; + u(IM1, k, j, i) = 0.0; + u(IM2, k, j, i) = 0.0; + u(IM3, k, j, i) = 0.0; + u(IEN, k, j, i) = P_r / gm1; + }); + } + + if (hydro_pkg->Param("fluid") == Fluid::glmmhd) { + /************************************************************ + * Initialize the initial magnetic field state via a vector potential + ************************************************************/ + parthenon::ParArray4D A("A", 3, pmb->cellbounds.ncellsk(IndexDomain::entire), + pmb->cellbounds.ncellsj(IndexDomain::entire), + pmb->cellbounds.ncellsi(IndexDomain::entire)); + + IndexRange a_ib = ib; + a_ib.s -= 1; + a_ib.e += 1; + IndexRange a_jb = jb; + a_jb.s -= 1; + a_jb.e += 1; + IndexRange a_kb = kb; + a_kb.s -= 1; + a_kb.e += 1; + + /************************************************************ + * Initialize an initial magnetic tower + ************************************************************/ + const auto &magnetic_tower = hydro_pkg->Param("magnetic_tower"); + + magnetic_tower.AddInitialFieldToPotential(pmb.get(), a_kb, a_jb, a_ib, A); + + /************************************************************ + * Add dipole magnetic field to the magnetic potential + ************************************************************/ + const auto &init_dipole_b_field = hydro_pkg->Param("init_dipole_b_field"); + if (init_dipole_b_field) { + const Real mx = hydro_pkg->Param("dipole_b_field_mx"); + const Real my = hydro_pkg->Param("dipole_b_field_my"); + const Real mz = hydro_pkg->Param("dipole_b_field_mz"); + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "MagneticTower::AddInitialFieldToPotential", + parthenon::DevExecSpace(), a_kb.s, a_kb.e, a_jb.s, a_jb.e, a_ib.s, a_ib.e, + KOKKOS_LAMBDA(const int &k, const int &j, const int &i) { + // Compute and apply potential + const Real x = coords.Xc<1>(i); + const Real y = coords.Xc<2>(j); + const Real z = coords.Xc<3>(k); + + const Real r3 = pow(SQR(x) + SQR(y) + SQR(z), 3. / 2); + + const Real m_cross_r_x = my * z - mz * y; + const Real m_cross_r_y = mz * x - mx * z; + const Real m_cross_r_z = mx * y - mx * y; + + A(0, k, j, i) += m_cross_r_x / (4 * M_PI * r3); + A(1, k, j, i) += m_cross_r_y / (4 * M_PI * r3); + A(2, k, j, i) += m_cross_r_z / (4 * M_PI * r3); + }); + } + + /************************************************************ + * Apply the potential to the conserved variables + ************************************************************/ + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "cluster::ProblemGenerator::ApplyMagneticPotential", + parthenon::DevExecSpace(), kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int &k, const int &j, const int &i) { + u(IB1, k, j, i) = + (A(2, k, j + 1, i) - A(2, k, j - 1, i)) / coords.Dxc<2>(j) / 2.0 - + (A(1, k + 1, j, i) - A(1, k - 1, j, i)) / coords.Dxc<3>(k) / 2.0; + u(IB2, k, j, i) = + (A(0, k + 1, j, i) - A(0, k - 1, j, i)) / coords.Dxc<3>(k) / 2.0 - + (A(2, k, j, i + 1) - A(2, k, j, i - 1)) / coords.Dxc<1>(i) / 2.0; + u(IB3, k, j, i) = + (A(1, k, j, i + 1) - A(1, k, j, i - 1)) / coords.Dxc<1>(i) / 2.0 - + (A(0, k, j + 1, i) - A(0, k, j - 1, i)) / coords.Dxc<2>(j) / 2.0; + + u(IEN, k, j, i) += 0.5 * (SQR(u(IB1, k, j, i)) + SQR(u(IB2, k, j, i)) + + SQR(u(IB3, k, j, i))); + }); + + /************************************************************ + * Add uniform magnetic field to the conserved variables + ************************************************************/ + const auto &init_uniform_b_field = hydro_pkg->Param("init_uniform_b_field"); + if (init_uniform_b_field) { + const Real bx = hydro_pkg->Param("uniform_b_field_bx"); + const Real by = hydro_pkg->Param("uniform_b_field_by"); + const Real bz = hydro_pkg->Param("uniform_b_field_bz"); + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "cluster::ProblemGenerator::ApplyUniformBField", + parthenon::DevExecSpace(), kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int &k, const int &j, const int &i) { + const Real bx_i = u(IB1, k, j, i); + const Real by_i = u(IB2, k, j, i); + const Real bz_i = u(IB3, k, j, i); + + u(IB1, k, j, i) += bx; + u(IB2, k, j, i) += by; + u(IB3, k, j, i) += bz; + + // Old magnetic energy is b_i^2, new Magnetic energy should be 0.5*(b_i + + // b)^2, add b_i*b + 0.5b^2 to old energy to accomplish that + u(IEN, k, j, i) += + bx_i * bx + by_i * by + bz_i * bz + 0.5 * (SQR(bx) + SQR(by) + SQR(bz)); + }); + // end if(init_uniform_b_field) + } + + } // END if(hydro_pkg->Param("fluid") == Fluid::glmmhd) + } + + /************************************************************ + * Initial parameters + ************************************************************/ + + auto pmb = md->GetBlockData(0)->GetBlockPointer(); + IndexRange ib = pmb->cellbounds.GetBoundsI(IndexDomain::interior); + IndexRange jb = pmb->cellbounds.GetBoundsJ(IndexDomain::interior); + IndexRange kb = pmb->cellbounds.GetBoundsK(IndexDomain::interior); + auto hydro_pkg = pmb->packages.Get("Hydro"); + const auto fluid = hydro_pkg->Param("fluid"); + auto const &cons = md->PackVariables(std::vector{"cons"}); + const auto num_blocks = md->NumBlocks(); + + /************************************************************ + * Set initial density perturbations + ************************************************************/ + + const bool init_perturb_rho = pin->GetOrAddBoolean("problem/cluster/init_perturb", "init_perturb_rho", false); + hydro_pkg->AddParam<>("init_perturb_rho", init_perturb_rho); + + if (init_perturb_rho == true) { + + // File load + + auto init_perturb_rho_file = pin->GetString("problem/cluster/init_perturb", "init_perturb_rho_file"); + auto init_perturb_rho_keys = pin->GetString("problem/cluster/init_perturb", "init_perturb_rho_keys"); + + hydro_pkg->AddParam<>("cluster/init_perturb_rho_file", init_perturb_rho_file); + hydro_pkg->AddParam<>("cluster/init_perturb_rho_keys", init_perturb_rho_keys); + + std::cout << "Setting density perturbation"; + + // Read HDF5 file containing the density + std::string filename_rho = "/work/bbd0833/test/rho.h5"; + std::string keys_rho = "data"; + H5Easy::File file(filename_rho, HighFive::File::ReadOnly); + auto rho_init = H5Easy::load, 64>, 64>>(file, keys_rho); + + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "Cluster::ProblemGenerator::UniformGas", + parthenon::DevExecSpace(), kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int &k, const int &j, const int &i) { + u(IDN, k, j, i) = rho_init[k-2][j-2][i-2]; + }); + + /* + auto init_perturb_rho_file = pin->GetString("problem/cluster/init_perturb", "init_perturb_rho_file"); + auto init_perturb_rho_keys = pin->GetString("problem/cluster/init_perturb", "init_perturb_rho_keys"); + + hydro_pkg->AddParam<>("cluster/init_perturb_rho_file", init_perturb_rho_file); + hydro_pkg->AddParam<>("cluster/init_perturb_rho_keys", init_perturb_rho_keys); + + std::cout << "Setting density perturbation"; + + // Read HDF5 file containing the density + std::string filename_rho = "/work/bbd0833/test/rho.h5"; + std::string keys_rho = "data"; + H5Easy::File file(filename_rho, HighFive::File::ReadOnly); + auto rho_init = H5Easy::load, 64>, 64>>(file, keys_rho); + + Real passive_scalar = 0.0; // Useless + + pmb->par_reduce( + "Init density field", 0, num_blocks - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int b, const int k, const int j, const int i, Real &lsum) { + const auto &coords = cons.GetCoords(b); + const auto &u = cons(b); + + // Adding the value of the density field + + u(IDN, k, j, i) = rho_init[k-2][j-2][i-2]; + + if (rho_init[k-2][j-2][i-2] < 0.0) { + std::cout << "Negative density for rho_init[" << k-2 << "][" << j-2 << "][" << i-2 << "]: " << rho_init[k-2][j-2][i-2] << "\n";} + + }, + passive_scalar); + + */ + } + + + /************************************************************ + * Set initial velocity perturbations (requires no other velocities for now) + ************************************************************/ + + const auto sigma_v = pin->GetOrAddReal("problem/cluster/init_perturb", "sigma_v", 0.0); + + if (sigma_v != 0.0) { + auto few_modes_ft = hydro_pkg->Param("cluster/few_modes_ft_v"); + // Init phases on all blocks + for (int b = 0; b < md->NumBlocks(); b++) { + auto pmb = md->GetBlockData(b)->GetBlockPointer(); + few_modes_ft.SetPhases(pmb.get(), pin); + + } + // As for t_corr in few_modes_ft, the choice for dt is + // in principle arbitrary because the inital v_hat is 0 and the v_hat_new will contain + // the perturbation (and is normalized in the following to get the desired sigma_v) + const Real dt = 1.0; + few_modes_ft.Generate(md, dt, "tmp_perturb"); + + Real v2_sum = 0.0; // used for normalization + + auto perturb_pack = md->PackVariables(std::vector{"tmp_perturb"}); + + pmb->par_reduce( + "Init sigma_v", 0, num_blocks - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int b, const int k, const int j, const int i, Real &lsum) { + const auto &coords = cons.GetCoords(b); + const auto &u = cons(b); + auto rho = u(IDN, k, j, i); + // The following restriction could be lifted, but requires refactoring of the + // logic for the normalization/reduction below + PARTHENON_REQUIRE( + u(IM1, k, j, i) == 0.0 && u(IM2, k, j, i) == 0.0 && u(IM3, k, j, i) == 0.0, + "Found existing non-zero velocity when setting velocity perturbations."); + + u(IM1, k, j, i) = rho * perturb_pack(b, 0, k, j, i); + u(IM2, k, j, i) = rho * perturb_pack(b, 1, k, j, i); + u(IM3, k, j, i) = rho * perturb_pack(b, 2, k, j, i); + // No need to touch the energy yet as we'll normalize later + + lsum += (SQR(u(IM1, k, j, i)) + SQR(u(IM2, k, j, i)) + SQR(u(IM3, k, j, i))) * + coords.CellVolume(k, j, i) / SQR(rho); + }, + v2_sum); + +#ifdef MPI_PARALLEL + PARTHENON_MPI_CHECK(MPI_Allreduce(MPI_IN_PLACE, &v2_sum, 1, MPI_PARTHENON_REAL, + MPI_SUM, MPI_COMM_WORLD)); +#endif // MPI_PARALLEL + + const auto Lx = pmesh->mesh_size.x1max - pmesh->mesh_size.x1min; + const auto Ly = pmesh->mesh_size.x2max - pmesh->mesh_size.x2min; + const auto Lz = pmesh->mesh_size.x3max - pmesh->mesh_size.x3min; + auto v_norm = std::sqrt(v2_sum / (Lx * Ly * Lz) / (SQR(sigma_v))); + + pmb->par_for( + "Norm sigma_v", 0, num_blocks - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) { + const auto &u = cons(b); + + u(IM1, k, j, i) /= v_norm; + u(IM2, k, j, i) /= v_norm; + u(IM3, k, j, i) /= v_norm; + + u(IEN, k, j, i) += + 0.5 * (SQR(u(IM1, k, j, i)) + SQR(u(IM2, k, j, i)) + SQR(u(IM3, k, j, i))) / + u(IDN, k, j, i); + }); + } + + /************************************************************ + * Set initial magnetic field perturbations (resets magnetic field field) + ************************************************************/ + const auto sigma_b = pin->GetOrAddReal("problem/cluster/init_perturb", "sigma_b", 0.0); + if (sigma_b != 0.0) { + auto few_modes_ft = hydro_pkg->Param("cluster/few_modes_ft_b"); + // Init phases on all blocks + for (int b = 0; b < md->NumBlocks(); b++) { + auto pmb = md->GetBlockData(b)->GetBlockPointer(); + few_modes_ft.SetPhases(pmb.get(), pin); + } + // As for t_corr in few_modes_ft, the choice for dt is + // in principle arbitrary because the inital b_hat is 0 and the b_hat_new will contain + // the perturbation (and is normalized in the following to get the desired sigma_b) + const Real dt = 1.0; + few_modes_ft.Generate(md, dt, "tmp_perturb"); + + Real b2_sum = 0.0; // used for normalization + + auto perturb_pack = md->PackVariables(std::vector{"tmp_perturb"}); + + pmb->par_reduce( + "Init sigma_b", 0, num_blocks - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int b, const int k, const int j, const int i, Real &lsum) { + const auto &coords = cons.GetCoords(b); + const auto &u = cons(b); + // The following restriction could be lifted, but requires refactoring of the + // logic for the normalization/reduction below + PARTHENON_REQUIRE( + u(IB1, k, j, i) == 0.0 && u(IB2, k, j, i) == 0.0 && u(IB3, k, j, i) == 0.0, + "Found existing non-zero B when setting magnetic field perturbations."); + u(IB1, k, j, i) = + (perturb_pack(b, 2, k, j + 1, i) - perturb_pack(b, 2, k, j - 1, i)) / + coords.Dxc<2>(j) / 2.0 - + (perturb_pack(b, 1, k + 1, j, i) - perturb_pack(b, 1, k - 1, j, i)) / + coords.Dxc<3>(k) / 2.0; + u(IB2, k, j, i) = + (perturb_pack(b, 0, k + 1, j, i) - perturb_pack(b, 0, k - 1, j, i)) / + coords.Dxc<3>(k) / 2.0 - + (perturb_pack(b, 2, k, j, i + 1) - perturb_pack(b, 2, k, j, i - 1)) / + coords.Dxc<1>(i) / 2.0; + u(IB3, k, j, i) = + (perturb_pack(b, 1, k, j, i + 1) - perturb_pack(b, 1, k, j, i - 1)) / + coords.Dxc<1>(i) / 2.0 - + (perturb_pack(b, 0, k, j + 1, i) - perturb_pack(b, 0, k, j - 1, i)) / + coords.Dxc<2>(j) / 2.0; + + // No need to touch the energy yet as we'll normalize later + lsum += (SQR(u(IB1, k, j, i)) + SQR(u(IB2, k, j, i)) + SQR(u(IB3, k, j, i))) * + coords.CellVolume(k, j, i); + }, + b2_sum); + +#ifdef MPI_PARALLEL + PARTHENON_MPI_CHECK(MPI_Allreduce(MPI_IN_PLACE, &b2_sum, 1, MPI_PARTHENON_REAL, + MPI_SUM, MPI_COMM_WORLD)); +#endif // MPI_PARALLEL + + const auto Lx = pmesh->mesh_size.x1max - pmesh->mesh_size.x1min; + const auto Ly = pmesh->mesh_size.x2max - pmesh->mesh_size.x2min; + const auto Lz = pmesh->mesh_size.x3max - pmesh->mesh_size.x3min; + auto b_norm = std::sqrt(b2_sum / (Lx * Ly * Lz) / (SQR(sigma_b))); + + pmb->par_for( + "Norm sigma_b", 0, num_blocks - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) { + const auto &u = cons(b); + + u(IB1, k, j, i) /= b_norm; + u(IB2, k, j, i) /= b_norm; + u(IB3, k, j, i) /= b_norm; + + u(IEN, k, j, i) += + 0.5 * (SQR(u(IB1, k, j, i)) + SQR(u(IB2, k, j, i)) + SQR(u(IB3, k, j, i))); + }); + } +} + +void UserWorkBeforeOutput(MeshBlock *pmb, ParameterInput *pin) { + // get hydro + auto pkg = pmb->packages.Get("Hydro"); + const Real gam = pin->GetReal("hydro", "gamma"); + const Real gm1 = (gam - 1.0); + + // get prim vars + auto &data = pmb->meshblock_data.Get(); + auto const &prim = data->Get("prim").data; + + // get derived fields + auto &log10_radius = data->Get("log10_cell_radius").data; + auto &entropy = data->Get("entropy").data; + auto &mach_sonic = data->Get("mach_sonic").data; + auto &temperature = data->Get("temperature").data; + + // for computing temperature from primitives + auto units = pkg->Param("units"); + auto mbar_over_kb = pkg->Param("mbar_over_kb"); + auto mbar = mbar_over_kb * units.k_boltzmann(); + + // fill derived vars (*including ghost cells*) + auto &coords = pmb->coords; + IndexRange ib = pmb->cellbounds.GetBoundsI(IndexDomain::entire); + IndexRange jb = pmb->cellbounds.GetBoundsJ(IndexDomain::entire); + IndexRange kb = pmb->cellbounds.GetBoundsK(IndexDomain::entire); + + pmb->par_for( + "Cluster::UserWorkBeforeOutput", kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int k, const int j, const int i) { + // get gas properties + const Real rho = prim(IDN, k, j, i); + const Real v1 = prim(IV1, k, j, i); + const Real v2 = prim(IV2, k, j, i); + const Real v3 = prim(IV3, k, j, i); + const Real P = prim(IPR, k, j, i); + + // compute radius + const Real x = coords.Xc<1>(i); + const Real y = coords.Xc<2>(j); + const Real z = coords.Xc<3>(k); + const Real r2 = SQR(x) + SQR(y) + SQR(z); + log10_radius(k, j, i) = 0.5 * std::log10(r2); + + // compute entropy + const Real K = P / std::pow(rho / mbar, gam); + entropy(k, j, i) = K; + + const Real v_mag = std::sqrt(SQR(v1) + SQR(v2) + SQR(v3)); + const Real c_s = std::sqrt(gam * P / rho); // ideal gas EOS + const Real M_s = v_mag / c_s; + mach_sonic(k, j, i) = M_s; + + // compute temperature + temperature(k, j, i) = mbar_over_kb * P / rho; + }); + + if (pkg->Param("enable_cooling") == Cooling::tabular) { + auto &cooling_time = data->Get("cooling_time").data; + + // get cooling function + const cooling::TabularCooling &tabular_cooling = + pkg->Param("tabular_cooling"); + const auto cooling_table_obj = tabular_cooling.GetCoolingTableObj(); + + pmb->par_for( + "Cluster::UserWorkBeforeOutput::CoolingTime", kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int k, const int j, const int i) { + // get gas properties + const Real rho = prim(IDN, k, j, i); + const Real P = prim(IPR, k, j, i); + + // compute cooling time + const Real eint = P / (rho * gm1); + const Real edot = cooling_table_obj.DeDt(eint, rho); + cooling_time(k, j, i) = (edot != 0) ? -eint / edot : NAN; + }); + } + + if (pkg->Param("fluid") == Fluid::glmmhd) { + auto &plasma_beta = data->Get("plasma_beta").data; + auto &mach_alfven = data->Get("mach_alfven").data; + + pmb->par_for( + "Cluster::UserWorkBeforeOutput::MHD", kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int k, const int j, const int i) { + // get gas properties + const Real rho = prim(IDN, k, j, i); + const Real P = prim(IPR, k, j, i); + const Real Bx = prim(IB1, k, j, i); + const Real By = prim(IB2, k, j, i); + const Real Bz = prim(IB3, k, j, i); + const Real B2 = (SQR(Bx) + SQR(By) + SQR(Bz)); + + // compute Alfven mach number + const Real v_A = std::sqrt(B2 / rho); + const Real c_s = std::sqrt(gam * P / rho); // ideal gas EOS + mach_alfven(k, j, i) = mach_sonic(k, j, i) * c_s / v_A; + + // compute plasma beta + plasma_beta(k, j, i) = (B2 != 0) ? P / (0.5 * B2) : NAN; + }); + } +} + +} // namespace cluster