From 35dc92388709394aa4ff0074722115c021fb4e91 Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Wed, 24 Jun 2026 10:23:49 +0200 Subject: [PATCH 1/8] restore ctrl-c with astrapy --- src/astrapy.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/astrapy.cpp b/src/astrapy.cpp index 5edc9e2..a785811 100644 --- a/src/astrapy.cpp +++ b/src/astrapy.cpp @@ -12,6 +12,7 @@ #include // everything needed for embedding #include // for numpy arrays #include // For STL vectors and containers +#include #include #include #include @@ -108,8 +109,9 @@ AstraPy::AstraPy(Input &input) { // Check whether we need to start an interpreter if(ninstance==1) { astra::cout << "AstraPy: start Python interpreter." << std::endl; - py::initialize_interpreter(); + std::signal(SIGINT, SIG_DFL); // restore "terminate on Ctrl-C" signal handler + py::exec("import sys; print(f'AstraPy: Python Version: {sys.version}')"); py::exec("print(f'AstraPy: Executable Path: {sys.executable}')"); py::exec("print(f'AstraPy: Sys Path: {sys.path}')"); From c6d3fce25a88d451057f812f23eac6547981a1b5 Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Wed, 24 Jun 2026 10:24:02 +0200 Subject: [PATCH 2/8] add sod solution for tests --- pytools/sod.py | 279 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 279 insertions(+) create mode 100644 pytools/sod.py diff --git a/pytools/sod.py b/pytools/sod.py new file mode 100644 index 0000000..07a60d2 --- /dev/null +++ b/pytools/sod.py @@ -0,0 +1,279 @@ +""" +Created on Thu Mar 5 11:27:16 2020 + +@author: glesur +Source from +https://github.com/ibackus/sod-shocktube/blob/master/sod.py +""" + +import numpy as np +import scipy +import scipy.optimize + + +def sound_speed(gamma, pressure, density, dustFrac=0.): + """ + Calculate sound speed, scaled by the dust fraction according to: + + .. math:: + \widetilde{c}_s = c_s \sqrt{1 - \epsilon} + + Where :math:`\epsilon` is the dustFrac + """ + scale = np.sqrt(1 - dustFrac) + return np.sqrt(gamma * pressure/ density) * scale + +def shock_tube_function(p4, p1, p5, rho1, rho5, gamma, dustFrac=0.): + """ + Shock tube equation + """ + z = (p4 / p5 - 1.) + c1 = sound_speed(gamma, p1, rho1, dustFrac) + c5 = sound_speed(gamma, p5, rho5, dustFrac) + + gm1 = gamma - 1. + gp1 = gamma + 1. + g2 = 2. * gamma + + fact = gm1 / g2 * (c5 / c1) * z / np.sqrt(1. + gp1 / g2 * z) + fact = (1. - fact) ** (g2 / gm1) + + return p1 * fact - p4 + + +def calculate_regions(pl, ul, rhol, pr, ur, rhor, gamma=1.4, dustFrac=0.): + """ + Compute regions + :rtype : tuple + :return: returns p, rho and u for regions 1,3,4,5 as well as the shock speed + """ + # if pl > pr... + rho1 = rhol + p1 = pl + u1 = ul + rho5 = rhor + p5 = pr + u5 = ur + + # unless... + if pl < pr: + rho1 = rhor + p1 = pr + u1 = ur + rho5 = rhol + p5 = pl + u5 = ul + + # solve for post-shock pressure + p4 = scipy.optimize.fsolve(shock_tube_function, p1, (p1, p5, rho1, rho5, gamma))[0] + + # compute post-shock density and velocity + z = (p4 / p5 - 1.) + c5 = sound_speed(gamma, p5, rho5, dustFrac) + + gm1 = gamma - 1. + gp1 = gamma + 1. + gmfac1 = 0.5 * gm1 / gamma + gmfac2 = 0.5 * gp1 / gamma + + fact = np.sqrt(1. + gmfac2 * z) + + u4 = c5 * z / (gamma * fact) + rho4 = rho5 * (1. + gmfac2 * z) / (1. + gmfac1 * z) + + # shock speed + w = c5 * fact + + # compute values at foot of rarefaction + p3 = p4 + u3 = u4 + rho3 = rho1 * (p3 / p1)**(1. / gamma) + return (p1, rho1, u1), (p3, rho3, u3), (p4, rho4, u4), (p5, rho5, u5), w + + +def calc_positions(pl, pr, region1, region3, w, xi, t, gamma, dustFrac=0.): + """ + :return: tuple of positions in the following order -> + Head of Rarefaction: xhd, Foot of Rarefaction: xft, + Contact Discontinuity: xcd, Shock: xsh + """ + p1, rho1 = region1[:2] # don't need velocity + p3, rho3, u3 = region3 + c1 = sound_speed(gamma, p1, rho1, dustFrac) + c3 = sound_speed(gamma, p3, rho3, dustFrac) + + if pl > pr: + xsh = xi + w * t + xcd = xi + u3 * t + xft = xi + (u3 - c3) * t + xhd = xi - c1 * t + else: + # pr > pl + xsh = xi - w * t + xcd = xi - u3 * t + xft = xi - (u3 - c3) * t + xhd = xi + c1 * t + + return xhd, xft, xcd, xsh + + +def region_states(pl, pr, region1, region3, region4, region5): + """ + :return: dictionary (region no.: p, rho, u), except for rarefaction region + where the value is a string, obviously + """ + if pl > pr: + return {'Region 1': region1, + 'Region 2': 'RAREFACTION', + 'Region 3': region3, + 'Region 4': region4, + 'Region 5': region5} + else: + return {'Region 1': region5, + 'Region 2': region4, + 'Region 3': region3, + 'Region 4': 'RAREFACTION', + 'Region 5': region1} + + +def create_arrays(pl, pr, xl, xr, positions, state1, state3, state4, state5, + npts, gamma, t, xi, dustFrac=0.): + """ + :return: tuple of x, p, rho and u values across the domain of interest + """ + xhd, xft, xcd, xsh = positions + p1, rho1, u1 = state1 + p3, rho3, u3 = state3 + p4, rho4, u4 = state4 + p5, rho5, u5 = state5 + gm1 = gamma - 1. + gp1 = gamma + 1. + + x_arr = np.linspace(xl, xr, npts) + rho = np.zeros(npts, dtype=float) + p = np.zeros(npts, dtype=float) + u = np.zeros(npts, dtype=float) + c1 = sound_speed(gamma, p1, rho1, dustFrac) + if pl > pr: + for i, x in enumerate(x_arr): + if x < xhd: + rho[i] = rho1 + p[i] = p1 + u[i] = u1 + elif x < xft: + u[i] = 2. / gp1 * (c1 + (x - xi) / t) + fact = 1. - 0.5 * gm1 * u[i] / c1 + rho[i] = rho1 * fact ** (2. / gm1) + p[i] = p1 * fact ** (2. * gamma / gm1) + elif x < xcd: + rho[i] = rho3 + p[i] = p3 + u[i] = u3 + elif x < xsh: + rho[i] = rho4 + p[i] = p4 + u[i] = u4 + else: + rho[i] = rho5 + p[i] = p5 + u[i] = u5 + else: + for i, x in enumerate(x_arr): + if x < xsh: + rho[i] = rho5 + p[i] = p5 + u[i] = -u1 + elif x < xcd: + rho[i] = rho4 + p[i] = p4 + u[i] = -u4 + elif x < xft: + rho[i] = rho3 + p[i] = p3 + u[i] = -u3 + elif x < xhd: + u[i] = -2. / gp1 * (c1 + (xi - x) / t) + fact = 1. + 0.5 * gm1 * u[i] / c1 + rho[i] = rho1 * fact ** (2. / gm1) + p[i] = p1 * fact ** (2. * gamma / gm1) + else: + rho[i] = rho1 + p[i] = p1 + u[i] = -u1 + + return x_arr, p, rho, u + + +def solve(left_state, right_state, geometry, t, gamma=1.4, npts=500, + dustFrac=0.): + """ + Solves the Sod shock tube problem (i.e. riemann problem) of discontinuity + across an interface. + + Parameters + ---------- + left_state, right_state: tuple + A tuple of the state (pressure, density, velocity) on each side of the + shocktube barrier for the ICs. In the case of a dusty-gas, the density + should be the gas density. + geometry: tuple + A tuple of positions for (left boundary, right boundary, barrier) + t: float + Time to calculate the solution at + gamma: float + Adiabatic index for the gas. + npts: int + number of points for array of pressure, density and velocity + dustFrac: float + Uniform fraction for the gas, between 0 and 1. + + Returns + ------- + positions: dict + Locations of the important places (rarefaction wave, shock, etc...) + regions: dict + constant pressure, density and velocity states in distinct regions + values: dict + Arrays of pressure, density, and velocity as a function of position. + The density ('rho') is the gas density, which may differ from the + total density in a dusty-gas. + Also calculates the specific internal energy + """ + + pl, rhol, ul = left_state + pr, rhor, ur = right_state + xl, xr, xi = geometry + + # basic checking + if xl >= xr: + print('xl has to be less than xr!') + exit() + if xi >= xr or xi <= xl: + print('xi has in between xl and xr!') + exit() + + # calculate regions + region1, region3, region4, region5, w = \ + calculate_regions(pl, ul, rhol, pr, ur, rhor, gamma, dustFrac) + + regions = region_states(pl, pr, region1, region3, region4, region5) + + # calculate positions + x_positions = calc_positions(pl, pr, region1, region3, w, xi, t, gamma, + dustFrac) + + pos_description = ('Head of Rarefaction', 'Foot of Rarefaction', + 'Contact Discontinuity', 'Shock') + positions = dict(zip(pos_description, x_positions)) + + # create arrays + x, p, rho, u = create_arrays(pl, pr, xl, xr, x_positions, + region1, region3, region4, region5, + npts, gamma, t, xi, dustFrac) + + energy = p/(rho * (gamma - 1.0)) + rho_total = rho/(1.0 - dustFrac) + val_dict = {'x':x, 'p':p, 'rho':rho, 'u':u, 'energy':energy, + 'rho_total':rho_total} + + return positions, regions, val_dict From 8ee2306ddead1454b227da0d7a83f4bb1d52ce20 Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Wed, 24 Jun 2026 10:24:19 +0200 Subject: [PATCH 3/8] add mass diffusion to compressible hydro --- src/rightHandSide/compressible_hydro.hpp | 33 +++++++++++++++++++----- 1 file changed, 27 insertions(+), 6 deletions(-) diff --git a/src/rightHandSide/compressible_hydro.hpp b/src/rightHandSide/compressible_hydro.hpp index dfbc274..70dcf1f 100644 --- a/src/rightHandSide/compressible_hydro.hpp +++ b/src/rightHandSide/compressible_hydro.hpp @@ -33,7 +33,10 @@ class CompressibleHydro : public RightHandSide, Shear> { private: real nu; + real etaRho; int viscosityOrder{1}; + int etaRhoOrder{1}; + real rhoFloor; bool haveSourceTerm{false}; real Omega; real cs{1.0}; // Sound speed (for isothermal equation of state) @@ -80,6 +83,12 @@ CompressibleHydro::CompressibleHydro(Input &input, Grid *grid) : RightHan this->nu = input.GetOrSet("Physics","viscosity",0,1e-3); this->viscosityOrder = input.GetOrSet("Physics","viscosity",1,1); + + this->etaRho = input.GetOrSet("Physics","eta_rho",0,0.0); + this->etaRhoOrder = input.GetOrSet("Physics","eta_rho",1,1); + + this->rhoFloor = input.GetOrSet("Physics","rho_floor",0,1e-6); + this->Omega = input.GetOrSet("Physics","omega",0,0.0); this->haveSourceTerm = (input.CheckEntry("Physics","omega") >0 || Shear::isEnabled); this->cs = input.GetOrSet("Physics","cs",0,1.0); @@ -124,18 +133,23 @@ void CompressibleHydro::ExplicitStep(Field>& fldin, Fiel Shear &shear = this->shear; shear.Refresh(t); + + real rhoFloor = this->rhoFloor; astra_for("hydro_windup", 0,npr[IDIR],0,npr[JDIR],0,npr[KDIR], KOKKOS_LAMBDA(int64_t i, int64_t j, int64_t k) { real p1 = pr1(i, j, k); real p2 = pr2(i, j, k); real p3 = pr3(i, j, k); + + real rhor_val = std::fmax(rhor(i, j, k), rhoFloor); // Ensure density is above the floor value + rhor(i, j, k) = rhor_val; // Update the density field with the floored value // Possibly a check on the minimum allowable density to avoid division by zero? - wr11(i, j, k) = p1 * p1 / rhor(i, j, k); - wr12(i, j, k) = p1 * p2 / rhor(i, j, k); - wr13(i, j, k) = p1 * p3 / rhor(i, j, k); - wr22(i, j, k) = p2 * p2 / rhor(i, j, k); - wr23(i, j, k) = p2 * p3 / rhor(i, j, k); - wr33(i, j, k) = p3 * p3 / rhor(i, j, k); + wr11(i, j, k) = p1 * p1 / rhor_val; + wr12(i, j, k) = p1 * p2 / rhor_val; + wr13(i, j, k) = p1 * p3 / rhor_val; + wr22(i, j, k) = p2 * p2 / rhor_val; + wr23(i, j, k) = p2 * p3 / rhor_val; + wr33(i, j, k) = p3 * p3 / rhor_val; }); // Fourier transform back to spectral space @@ -218,9 +232,13 @@ void CompressibleHydro::ImplicitStep(Field>& fldin, real auto px1 = fldin["px1"]; auto px2 = fldin["px2"]; auto px3 = fldin["px3"]; + auto rho = fldin["rho"]; real nu= this->nu; int n = this->viscosityOrder; + real etaRho = this->etaRho; + int nRho = this->etaRhoOrder; + Shear shear = this->shear; shear.Refresh(t); astra_for("CompressibleHydro_viscosity", fldin, @@ -235,6 +253,9 @@ void CompressibleHydro::ImplicitStep(Field>& fldin, real px1(i,j,k) *= factor; px2(i,j,k) *= factor; px3(i,j,k) *= factor; + + real factorRho = std::exp(-dt * pow(etaRho*k2t, nRho)); // Exact integration + rho(i,j,k) *= factorRho; }); astra::popRegion(); } From cf54d88bccb573a52f8a869ac7484c73ff825644 Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Wed, 24 Jun 2026 10:24:57 +0200 Subject: [PATCH 4/8] add sod test for compressible hydro --- test/sod-iso/astra.ini | 28 ++++++++++++++ test/sod-iso/sod_test.py | 11 ++++++ test/sod-iso/validate.py | 79 ++++++++++++++++++++++++++++++++++++++++ 3 files changed, 118 insertions(+) create mode 100644 test/sod-iso/astra.ini create mode 100644 test/sod-iso/sod_test.py create mode 100755 test/sod-iso/validate.py diff --git a/test/sod-iso/astra.ini b/test/sod-iso/astra.ini new file mode 100644 index 0000000..5a8358d --- /dev/null +++ b/test/sod-iso/astra.ini @@ -0,0 +1,28 @@ +[Grid] +X1-grid -1 1536 2.0 +X2-grid -0.5 2 0.5 +X3-grid -0.5 2 0.5 + +[Python] +script sod_test + +[TimeIntegrator] +method rk3 +cfl 0.9 +tstop 0.2 + +[Physics] +rhs compressible_hydro +viscosity 4e-5 2 +eta_rho 4e-5 2 +rho_floor 1e-2 +cs 1.0 + +[InitFlow] +python init_flow real +# large_scale_3d_noise 1.0e-2 0.2 # Amplitude and cut-off length +# mean_field 0.0 0.0 3.0e-2 # Mean velocity in each direction + +[Output] +log 10 +vtk 0.1 diff --git a/test/sod-iso/sod_test.py b/test/sod-iso/sod_test.py new file mode 100644 index 0000000..c07ff86 --- /dev/null +++ b/test/sod-iso/sod_test.py @@ -0,0 +1,11 @@ +from astrapy import * +import numpy as np +#import matplotlib.pyplot as plt + + +def init_flow(grid, field): + x,y,z = np.meshgrid(grid.x[IDIR], grid.x[JDIR], grid.x[KDIR], indexing='ij') + field["rho"][:,:,:] = np.where(x<0.5, 1.0, 0.125) + field["px1"][:,:,:] = 0 + field["px2"][:,:,:] = 0 + field["px3"][:,:,:] = 0 diff --git a/test/sod-iso/validate.py b/test/sod-iso/validate.py new file mode 100755 index 0000000..611ebdc --- /dev/null +++ b/test/sod-iso/validate.py @@ -0,0 +1,79 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Mon Jun 21 15:42:19 2021 + +@author: lesurg +""" +import sys +import numpy as np +import matplotlib.pyplot as plt +from scipy.interpolate import interp1d +import argparse +import os +sys.path.append(os.path.join(os.path.dirname(os.path.abspath(__file__)), "../../pytools/")) +from vtk_io import readVTK +import sod + + +parser = argparse.ArgumentParser() +parser.add_argument("-noplot", + default=False, + help="disable plotting", + action="store_true") + + +args, unknown=parser.parse_known_args() + +gamma = 1.00000000001 +npts = 5000 + +# left_state and right_state set p, rho and u +# geometry sets left boundary on 0., right boundary on 1 and initial +# position of the shock xi on 0.5 +# t is the time evolution for which positions and states in tube should be calculated +# gamma denotes specific heat +# note that gamma and npts are default parameters (1.4 and 500) in solve function +positions, regions, values = sod.solve(left_state=(1, 1, 0), right_state=(0.125, 0.125, 0.), + geometry=(0., 1., 0.5), t=0.2, gamma=gamma, npts=npts) + + +# Finally, let's plot solutions +p = values['p'] +rho = values['rho'] +u = values['u'] +x= values['x'] + + +solinterp=interp1d(x,u) + +V=readVTK('data.0002.vtk', geometry='cartesian') +istart=np.argwhere(V.x>=0.0)[0][0] +iend=np.argwhere(V.x>1.0)[0][0] +x_simu=V.x[istart:iend] +v_simu=V.data['px1'][istart:iend,0,0]/V.data['rho'][istart:iend,0,0] +rho_simu=V.data['rho'][istart:iend,0,0] +error=np.mean(np.fabs((v_simu-solinterp(x_simu)))) + + +if(not args.noplot): + plt.figure(1) + plt.plot(x,rho) + plt.plot(x_simu,rho_simu,'+',markersize=2) + plt.title('Density') + + plt.figure(2) + plt.plot(x,u) + plt.plot(x_simu,v_simu,'+',markersize=2) + plt.title('Velocity') + + plt.ioff() + plt.show() + +print("Error=%e"%error) +if error<5e-3: + print("SUCCESS!") + sys.exit(0) +else: + print("FAILURE!") + sys.exit(1) From b8121da875336704c781fe66564248af6aa0c3cd Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Wed, 24 Jun 2026 10:27:28 +0200 Subject: [PATCH 5/8] add sod-iso to CI --- .github/workflows/astra-ci-jobs.yml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.github/workflows/astra-ci-jobs.yml b/.github/workflows/astra-ci-jobs.yml index 9f3dfd8..dfe0fee 100644 --- a/.github/workflows/astra-ci-jobs.yml +++ b/.github/workflows/astra-ci-jobs.yml @@ -56,6 +56,7 @@ jobs: source .venv/bin/activate cd job-serial python3 ../test/test.py run shearing_wave_compressible $TEST_OPTIONS + python3 ../test/test.py run sod-iso $TEST_OPTIONS Parallel: runs-on: self-hosted @@ -98,4 +99,5 @@ jobs: run: | source .venv/bin/activate cd job-parallel + python3 ../test/test.py run sod-iso $TEST_OPTIONS -mpi python3 ../test/test.py run shearing_wave_compressible $TEST_OPTIONS -mpi From 58434bb6cf21644ca920ee494f5b5816dd15d529 Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Wed, 24 Jun 2026 17:35:26 +0200 Subject: [PATCH 6/8] update documentation --- doc/source/input_file.rst | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/doc/source/input_file.rst b/doc/source/input_file.rst index d34846c..1a21beb 100644 --- a/doc/source/input_file.rst +++ b/doc/source/input_file.rst @@ -83,9 +83,7 @@ Depending on the choice of the right-hand side, the following entries may be req +-----------------+--------------------+-----------------------------------------------------------------------------------------------------------+ | omega | float | (optional) rotation rate along the x3 (=z) axis | +-----------------+--------------------+-----------------------------------------------------------------------------------------------------------+ -| viscosity_order | int | (optional) order of the viscosity term *n* in :math:`=\nu \Delta ^n v` | -+-----------------+--------------------+-----------------------------------------------------------------------------------------------------------+ -| shear_type | string | (optional) type of large-scale shear. Value allowed: ``linear`` | +| shear_type | string | (optional) type of large-scale shear. Value allowed: ``linear`` | +-----------------+--------------------+-----------------------------------------------------------------------------------------------------------+ | shear_rate | float | (optional) shear rate when `shear_type` is `linear` | +-----------------+--------------------+-----------------------------------------------------------------------------------------------------------+ @@ -107,9 +105,7 @@ Depending on the choice of the right-hand side, the following entries may be req +-----------------+--------------------+-----------------------------------------------------------------------------------------------------------+ | omega | float | (optional) rotation rate along the x3 (=z) axis | +-----------------+--------------------+-----------------------------------------------------------------------------------------------------------+ -| viscosity_order | int | (optional) order of the viscosity term *n* in :math:`=\nu \Delta ^n v` | -+-----------------+--------------------+-----------------------------------------------------------------------------------------------------------+ -| shear_type | string | (optional) type of large-scale shear. Value allowed: ``linear`` | +| shear_type | string | (optional) type of large-scale shear. Value allowed: ``linear`` | +-----------------+--------------------+-----------------------------------------------------------------------------------------------------------+ | shear_rate | float | (optional) shear rate when `shear_type` is `linear` | +-----------------+--------------------+-----------------------------------------------------------------------------------------------------------+ @@ -124,13 +120,16 @@ Depending on the choice of the right-hand side, the following entries may be req | viscosity | float, (int) | | 1st parameter: kinematic viscosity | | | | | 2nd parameter (optional): order of the viscosity term *n* in :math:`=\nu \Delta^n v (default 1)` | +-----------------+--------------------+-----------------------------------------------------------------------------------------------------------+ +| eta_rho | float, (int) | | 1st parameter: mass diffusion | +| | | | 2nd parameter (optional): order of the diffusion term *n* in :math:`=\eta_\rho \Delta^n v (default 1)` | ++-----------------+--------------------+-----------------------------------------------------------------------------------------------------------+ +| rho_floor | float | (optional) density floor (default 1e-6) | ++-----------------+--------------------+-----------------------------------------------------------------------------------------------------------+ | omega | float | (optional) rotation rate along the x3 (=z) axis | +-----------------+--------------------+-----------------------------------------------------------------------------------------------------------+ | cs | float | (optional) isothermal sound speed (default 1) | +-----------------+--------------------+-----------------------------------------------------------------------------------------------------------+ -| viscosity_order | int | (optional) order of the viscosity term *n* in :math:`=\nu \Delta ^n v` | -+-----------------+--------------------+-----------------------------------------------------------------------------------------------------------+ -| shear_type | string | (optional) type of large-scale shear. Value allowed: ``linear`` | +| shear_type | string | (optional) type of large-scale shear. Value allowed: ``linear`` | +-----------------+--------------------+-----------------------------------------------------------------------------------------------------------+ | shear_rate | float | (optional) shear rate when `shear_type` is `linear` | +-----------------+--------------------+-----------------------------------------------------------------------------------------------------------+ From 80e22ad43e2f5a3e72e1bdc4516b79ab3897a860 Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Wed, 24 Jun 2026 17:48:25 +0200 Subject: [PATCH 7/8] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- doc/source/input_file.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/source/input_file.rst b/doc/source/input_file.rst index 1a21beb..3ee0b96 100644 --- a/doc/source/input_file.rst +++ b/doc/source/input_file.rst @@ -121,7 +121,7 @@ Depending on the choice of the right-hand side, the following entries may be req | | | | 2nd parameter (optional): order of the viscosity term *n* in :math:`=\nu \Delta^n v (default 1)` | +-----------------+--------------------+-----------------------------------------------------------------------------------------------------------+ | eta_rho | float, (int) | | 1st parameter: mass diffusion | -| | | | 2nd parameter (optional): order of the diffusion term *n* in :math:`=\eta_\rho \Delta^n v (default 1)` | +| | | | 2nd parameter (optional): order of the diffusion term *n* in :math:`=\eta_\rho \Delta^n \rho (default 1)` | +-----------------+--------------------+-----------------------------------------------------------------------------------------------------------+ | rho_floor | float | (optional) density floor (default 1e-6) | +-----------------+--------------------+-----------------------------------------------------------------------------------------------------------+ From ab700f15659d80a53d171f52a111ee80b88f0d7a Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Wed, 24 Jun 2026 17:48:36 +0200 Subject: [PATCH 8/8] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- src/rightHandSide/compressible_hydro.hpp | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/rightHandSide/compressible_hydro.hpp b/src/rightHandSide/compressible_hydro.hpp index 70dcf1f..e45b726 100644 --- a/src/rightHandSide/compressible_hydro.hpp +++ b/src/rightHandSide/compressible_hydro.hpp @@ -141,9 +141,7 @@ void CompressibleHydro::ExplicitStep(Field>& fldin, Fiel real p2 = pr2(i, j, k); real p3 = pr3(i, j, k); - real rhor_val = std::fmax(rhor(i, j, k), rhoFloor); // Ensure density is above the floor value - rhor(i, j, k) = rhor_val; // Update the density field with the floored value - // Possibly a check on the minimum allowable density to avoid division by zero? + const real rhor_val = std::fmax(rhor(i, j, k), rhoFloor); // avoid division by ~0 when computing p_i p_j / rho wr11(i, j, k) = p1 * p1 / rhor_val; wr12(i, j, k) = p1 * p2 / rhor_val; wr13(i, j, k) = p1 * p3 / rhor_val;