Add input with proper energ conservation. Modified filtering frequency when time derivative is computed

Author: VARENNES Robin

diagnostics/diag_main.ipynb

@@ -149,6 +149,8 @@
     "\n",
     "from simulation_diag_handler import Simulation, set_plot_defaults\n",
     "\n",
+    "sim = Simulation(\"/ZONE_TRAVAIL/RV255575/tokam2d_conservation/HW_cons_rst\")\n",
+    "\n",
     "set_plot_defaults()\n",
     "\n",
     "plt.close('all')\n",
@@ -290,6 +292,160 @@
     "\n",
     "sim.make_movie(field='density', path=None, filename='mysim', it_slice=None, parallel=True, num_cores=None, for_IA=False, scheme=None, cmap='plasma', vmin=None, vmax=None, fps=30, save_frames=False)"
    ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "869be0f0",
+   "metadata": {},
+   "source": [
+    "## **Conservation**"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "175e983b",
+   "metadata": {},
+   "source": [
+    "### **Compute energy conservation (Hasegawa-Wakatani)**\n",
+    "\n",
+    "\\begin{equation}\n",
+    "\\frac{dE}{dt} = \\Gamma_n - D_\\alpha - D_E\n",
+    "\\end{equation}\n",
+    "\n",
+    "Where \n",
+    "\n",
+    "\\begin{equation}\n",
+    "E = \\frac{1}{2} \\int n^2 + |\\nabla \\phi|^2 \\, d\\mathbf{x}\n",
+    "\\end{equation}\n",
+    "\n",
+    "The time derivative of the energy $E$ can be computed from saved data as\n",
+    "\n",
+    "\\begin{equation}\n",
+    "\\frac{dE}{dt} = \\int \\left( n \\cdot \\frac{dn}{dt} \n",
+    "+ |\\nabla_x \\phi| \\cdot \\frac{d |\\nabla_x \\phi|}{dt}\n",
+    "+ |\\nabla_y \\phi| \\cdot \\frac{d |\\nabla_y \\phi|}{dt} \\right)\n",
+    "\\, d\\mathbf{x}\n",
+    "\\end{equation}\n",
+    "\n",
+    "and \n",
+    "\n",
+    "\\begin{equation}\n",
+    "\\Gamma_n = -\\kappa \\int n \\cdot \\frac{\\partial \\phi}{\\partial y} \\, d\\mathbf{x} \\text{~~(``source'')}\n",
+    "\\end{equation}\n",
+    "\n",
+    "\\begin{equation}\n",
+    "D_C = C \\int (n - \\phi)^2 \\, d\\mathbf{x} \\text{~~(``coupling'')}\n",
+    "\\end{equation}\n",
+    "\n",
+    "\\begin{equation}\n",
+    "D_{E,n} = (-1)^{N_n} D_n \\int n \\nabla^{2N_n} n \\, d\\mathbf{x} \\text{~~(``density dissipation'')}\n",
+    "\\end{equation}\n",
+    "\n",
+    "\\begin{equation}\n",
+    "D_{E,\\Omega} = (-1)^{N_\\phi} D_\\phi \\int n \\nabla^{2N_\\phi+2} \\phi \\, d\\mathbf{x} \\text{~~(``vorticity dissipation'')}\n",
+    "\\end{equation}\n",
+    "\n",
+    "**<span style=\"color:red\">\n",
+    "/!\\ Some important notes\n",
+    "</span>**\n",
+    "\n",
+    "- Accuracy of time derivatives are very important to compute the conservations. It is advised to switch on the 'compute_time_derivatives' in the simulation inputs (which computes the 8th order central finite difference from the rk4 time steps which are not available in the stored data outputs) and derive the time derivative of energy according to the formula above.\n",
+    "\n",
+    "- Semi-spectral code such as Tokam2D usually relies on a 2/3 dealiasing rule for avoid the emergence of spurious (non-physical) modes by truncating a part of high k (small scales) while the simulation is running.\n",
+    "One consequence is that the energy transferred to these modes is lost.\n",
+    "While the difference is (or should be) small, test cases for evaluating conservation properties of the code can be performed without this filter by switching of the 'fft_filter' in input file.\n",
+    "\n",
+    "- For Hasegawa-Wakatani schemes, noise accumulates during Tokam2D simulation and takes the form of non-zero imaginary part for some modes for the real space fields (which should be real signals). A filter is thus applied frequently during the simulation to remove this noise. Consequence on conservation are not clear.\n",
+    "\n",
+    "- Depending on the level of accuracy of conservation one needs, diminishing 'dt_rk4' (or equivalently increasing 'rk4_per_diag') in the input file should improve agreement. Increasing the spatial resolution 'Nx' and 'Ny', however, does not garantee a better conservation"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9f70cf8a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from simulation_diag_handler import Simulation, set_plot_defaults\n",
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "set_plot_defaults()\n",
+    "\n",
+    "## Fetch simulation data\n",
+    "sim_path = '/path/to/sim'\n",
+    "sim = Simulation(sim_path)\n",
+    "\n",
+    "## Define proxy function\n",
+    "def rifft(x): return np.real(np.fft.ifft2(x, axes=(-2, -1)))\n",
+    "\n",
+    "## Choose\n",
+    "it_beg = 0\n",
+    "it_end = len(sim.time) - 1\n",
+    "\n",
+    "## Cast useful quantities\n",
+    "potential_fft    = sim.get_data_slice('potential_fft'   , slice(it_beg, it_end))\n",
+    "density_fft      = sim.get_data_slice('density_fft'     , slice(it_beg, it_end))\n",
+    "density          = sim.get_data_slice('density'         , slice(it_beg, it_end))\n",
+    "potential        = sim.get_data_slice('potential'       , slice(it_beg, it_end))\n",
+    "dt_density       = sim.get_data_slice('dt_density'      , slice(it_beg, it_end))\n",
+    "dt_potential     = sim.get_data_slice('dt_potential'    , slice(it_beg, it_end))\n",
+    "dt_potential_fft = sim.get_data_slice('dt_potential_fft', slice(it_beg, it_end))\n",
+    "kappa = sim[\"kappa\"][()]\n",
+    "C     = sim[\"C\"][()]\n",
+    "D_n   = sim[\"Dn\"][()][0]\n",
+    "D_phi = sim[\"Dphi\"][()][0]\n",
+    "N_n   = 2 # 1: Diffusion, 2: Hyperdiffusion\n",
+    "N_phi = 2 # 1: Diffusion, 2: Hyperdiffusion\n",
+    "kx2D  = sim.kx2D\n",
+    "ky2D  = sim.ky2D\n",
+    "k2D   = sim.k2D\n",
+    "\n",
+    "## Compute energy\n",
+    "gradx_phi    = rifft( (1j * kx2D) * potential_fft)\n",
+    "grady_phi    = rifft( (1j * ky2D) * potential_fft)\n",
+    "energy = 0.5*np.mean(density**2 + (gradx_phi**2 + grady_phi**2), axis=(-2, -1))\n",
+    "\n",
+    "# FIGURE : time evolution of energy\n",
+    "fig = plt.figure(figsize=(10,7))\n",
+    "ax = fig.add_subplot(111)\n",
+    "ax.plot(sim.time[it_beg:it_end], energy)\n",
+    "ax.set_xlabel('Time')\n",
+    "ax.set_ylabel('Energy')\n",
+    "ax.set_title(f'Averaged between {sim.time[it_beg]:.2f} < t < {sim.time[it_end-1]:.2f}')\n",
+    "\n",
+    "## Compute terms of energy conservation equation \n",
+    "# Energy time derivative\n",
+    "dEdt_n       = density * dt_density\n",
+    "gradx_dt_phi = rifft( (1j * kx2D) * dt_potential_fft)\n",
+    "grady_dt_phi = rifft( (1j * ky2D) * dt_potential_fft)\n",
+    "dEdt_phi      = gradx_phi * gradx_dt_phi + grady_phi * grady_dt_phi\n",
+    "dEdt_avg      = np.mean(dEdt_n + dEdt_phi, axis=(-2, -1))\n",
+    "\n",
+    "# Source term\n",
+    "Gamma_n = - kappa * np.mean( density * grady_phi , axis=(-2, -1))\n",
+    "\n",
+    "# Coupling term\n",
+    "D_C = C * np.mean( (density - potential)**2 , axis=(-2, -1))\n",
+    "\n",
+    "# Density dissipation\n",
+    "nabla_2N_dens = rifft( (1j * k2D)**(2*N_n) * density_fft)\n",
+    "DE_n = (-1)**N_n * D_n * np.mean(density * nabla_2N_dens, axis=(-2, -1))\n",
+    "\n",
+    "# Vorticity dissipation\n",
+    "nabla_2N_vort = rifft( (1j * k2D)**(2*N_phi+2) * potential_fft)\n",
+    "DE_vort = (-1)**N_phi * D_phi * np.mean(potential * nabla_2N_vort, axis=(-2, -1))\n",
+    "\n",
+    "# FIGURE : time evolution of energy conservation LHS and RHS\n",
+    "fig = plt.figure(figsize=(10,7))\n",
+    "ax = fig.add_subplot(111)\n",
+    "ax.plot(sim.time[it_beg:it_end], dEdt_avg, ls='-',label=r'LHS: $\\partial_t E$')\n",
+    "ax.plot(sim.time[it_beg:it_end], Gamma_n - D_C - DE_n + DE_vort, ls='--',label=r'RHS: $\\Gamma_n - D_C - D_{E,n} + D_{E,\\Omega}$')\n",
+    "ax.set_xlabel('Time')\n",
+    "ax.set_title(f'Averaged between {sim.time[it_beg]:.2f} < t < {sim.time[it_end-1]:.2f}')\n",
+    "ax.legend()"
+   ]
   }
  ],
  "metadata": {

diagnostics/simulation_diag_handler.py

@@ -242,7 +242,11 @@ def get_data_slice(self, field, it=None, iy=None, ix=None):
             return np.fft.fft2(self["potential"][it_slice], axes=(-2, -1))[..., y_slice, x_slice]
         if field == "density_fft":
             return np.fft.fft2(self["density"][it_slice], axes=(-2, -1))[..., y_slice, x_slice]
-
+        if field == "dt_potential_fft":
+            return np.fft.fft2(self["dt_potential"][it_slice], axes=(-2, -1))[..., y_slice, x_slice]
+        if field == "dt_density_fft":
+            return np.fft.fft2(self["dt_density"][it_slice], axes=(-2, -1))[..., y_slice, x_slice]
+        
         # If the name of the field is not in the mapping, raise an error
         if field not in self.field_mapping:
             raise KeyError(f"Field {field} not found in mapping.")

main.py

@@ -10,7 +10,7 @@
 from src.interfaces.output_saver import OutputSaver
 
 import os
-os.environ["CUDA_VISIBLE_DEVICES"] = "0"
+os.environ["CUDA_VISIBLE_DEVICES"] = "1"
 
 import time
 

src/simulation/run_simulation.py

@@ -114,6 +114,8 @@ def run(self, fields):
         while self.step_rk4_count < self.Nt_rk4 and not self._stop:
 
             if self.user_inline_operations.get('compute_time_derivatives', False):
+                for op in self.inline_operations:
+                    fields = op(fields)
                 self.inline_compute_time_derivative(fields)
 
             if self.step_rk4_count % self.rk4_per_diag == 0:
@@ -128,8 +130,9 @@ def run(self, fields):
                 # Trigger inline operations (e.g. FFT filtering) at diagnostic time step
                 for op in self.inline_operations:
                     fields = op(fields)
-
+                    
                 self.step_diag_count += 1
+
             fields = step(fields, self.pde, self.dt_rk4, t=self.time)
 
             self.time += self.dt_rk4