Computer Physics Communications Research Articles

Analysis of the multi-term fractional Bateman equations in radioactive decay by means of Mikusiński algebraic calculus

Fractional Bateman equations model the concentrations of nuclear species through a radioactive decay chain, by considering non-Markovian behavior. Cruz-López and Espinosa-Paredes (Nuclear Engineering and Technology, 2022) solved the mathematical problem with a single memory coefficient, based on the Laplace-transform technique. Cruz-López, Espinosa-Paredes and François (Computer Physics Communications, 2022) extended the analysis to several fractional orders, with Laplace transforms and series expansions checked by induction. In this paper, we solve such a case of distinct memory effects, via the more general method of Mikusiński operational calculus. We do not distinguish between equal and different decay rates. Multivariate Mittag-Leffler functions appear in the closed-form expression.

Enabling particle transport on CAD-based geometries for radiation simulations with penRed

Geometry construction is a fundamental aspect of any radiation transport simulation, regardless of the Monte Carlo code being used. Typically, this process is tedious, time-consuming, and error-prone. The conventional approach involves defining geometries using mathematical objects or surfaces. However, this method comes with several limitations, especially when dealing with complex models, particularly those with organic shapes. Furthermore, since each code employs its own format and methodology for defining geometries, sharing and reproducing simulations among researchers becomes a challenging task. Consequently, many codes have implemented support for simulating over geometries constructed via Computer-Aided Design (CAD) tools. Unfortunately, this feature is lacking in penRed and other PENELOPE physics-based codes. Therefore, the objective of this work is to implement such support within the penRed framework. New version program summaryProgram Title: Parallel Engine for Radiation Energy Deposition (penRed)CPC Library link to program files:https://doi.org/10.17632/rkw6tvtngy.2Developer's repository link:https://github.com/PenRed/PenRedCode Ocean capsule:https://codeocean.com/capsule/1041417/treeLicensing provisions: GNU Affero General Public License v3Programming language: C++ standard 2011.Journal reference of previous version: V. Giménez-Alventosa, V. Giménez Gómez, S. Oliver, PenRed: An extensible and parallel Monte-Carlo framework for radiation transport based on PENELOPE, Computer Physics Communications 267 (2021) 108065. doi:https://doi.org/10.1016/j.cpc.2021.108065.Does the new version supersede the previous version?: YesReasons for the new version: Implements the capability to simulate on CAD constructed geometries, among many other features and fixes.Summary of revisions: All changes applied through the code versions are summarized in the file CHANGELOG.md in the repository package.Nature of problem: While Monte Carlo codes have proven valuable in simulating complex radiation scenarios, they rely heavily on accurate geometrical representations. In the same way as many other Monte Carlo codes, penRed employs simple geometric quadric surfaces like planes, spheres and cylinders to define geometries. However, since these geometric models offer a certain level of flexibility, these representations have limitations when it comes to simulating highly intricate and irregular shapes. Anatomic structures, for example, require detailed representations of organs, tissues and bones, which are difficult to achieve using basic geometric objects. Similarly, complex devices or intricate mechanical systems may have designs that cannot be accurately represented within the constraints of such geometric models. Moreover, when the complexity of the model increases, geometry construction process becomes more difficult, tedious, time-consuming and error-prone [2]. Also, as each Monte Carlo geometry library uses its own format and construction method, reproducing the same geometry among different codes is a challenging task.Solution method: To face the problems stated above, the objective of this work is to implement the capability to simulate using irregular and adaptable meshed geometries in the penRed framework. This kind of meshes can be constructed using Computer-Aided Design (CAD) tools, the use of which is very widespread and streamline the design process. This feature has been implemented in a new geometry module named “MESH_BODY” specific for this kind of geometries. This one is freely available and usable within the official penRed package1. It can be used since penRed version 1.9.3b and above.

OpenFuelCell2: A new computational tool for fuel cells, electrolyzers, and other electrochemical devices and processes

Fuel cells/electrolyzers are efficient and clean electrochemical devices that convert chemical energy directly into electricity and vice versa. They have attracted sustainable attention over the past decade from multiple experimental and numerical studies. However, detailed experimental investigations are typically expensive and challenging for providing a number of operating conditions and designs. Computational analysis offers an alternative approach for these studies. With the steadily increasing high-performance computing resources available, the limitations of numerical simulations have substantially decreased. This contribution details the design choice and code structure of modern electrochemical devices, which have been implemented as a versatile C++ library named openFuelCell2 within the open-source platform OpenFOAM, allowing for large-scale parallel calculations to be performed. The solver considers the major transport phenomena in a typical electrochemical device, including fluid flow, heat and mass transfer, species and charge transfer, and electrochemical reaction. This enables numerical simulations on popular electrochemical devices, such as fuel cells and electrolyzers, to be conducted. The paper also describes the domain decomposition, and parallel performance issues, as well as future applications. Program summaryProgram Title:openFuelCell2CPC Library link to program files:https://doi.org/10.17632/cvmb5xgy7h.1Developer's repository link:https://github.com/openFuelCell2/openFuelCell2Licensing provisions: GPLv3Programming language: C++Journal reference of previous version: S.B. Beale, H.-W. Choi, J.G. Pharoah, H.K. Roth, H. Jasak, D.H. Jeon, Open-source computational model of a solid oxide fuel cell, Computer Physics Communications 200 (2016) 15–26.Does the new version supersede the previous version?: NoReasons for new version: This new version aims to account for additional electrochemical applications.Summary of revisions: The new version employs new code structure designs and is capable of considering two different electric potential fields and two-phase flows.Nature of problem: This software library provides a set of models for simulating electrochemical devices, such as fuel cells and electrolyzers, as well as other similar devices. These systems consist of multiple components with distinct features but are coupled through properties such as temperature, species concentrations, and electric current/potentials. This software utilizes a computational fluid dynamics (CFD) approach to enable the simulation of multiregion and multiphysics problems, accounting for heat and mass transfer, single and two-phase flow, multicomponent diffusion, electron and ion transfer, and electrochemical reactions. This toolbox provides researchers with a powerful tool for studying electrochemical devices in a comprehensive and efficient manner.Solution method: This software utilizes the finite volume method (FVM) to discretize and solve all conserved variables governed by partial differential equations. It employs a variant of the SIMPLE and PISO algorithms, known as PIMPLE, to solve the pressure-velocity coupling. All of the transport equations are solved sequentially, except for the enthalpy equation which is solved in an implicitly coupled manner within the global domain.Additional comments including restrictions and unusual features:•The two-phase flow consists of continuous and dispersed phases, without resolving the phase interfaces.•The two-phase flow is solved with an Eulerian-Eulerian approach.•A local-time-stepping method was typically used in the two-phase solution.•The electrochemical reactions take place in volumetric regions.•Species diffusion is described by Fick's law.

SATLAS2: An update to the package for analysis of counting data

SATLAS2 is a Python library that enables the user to fit counting data from laser spectroscopy experiments, in particular those that measure atomic hyperfine structures. In this analysis, the user can choose how the uncertainties are treated and can also opt to generate a random walk in order to present a fuller picture of the parameter space. The major upgrade compared to the previous version of SATLAS [1] is the different architecture of the codebase, which enabled a performance boost, with speed-up factors ranging from 20 to 300 times for various use cases. For backward compatibility, a translation layer between the two architectures is available, implementing only the core functionality of SATLAS. New version program summaryProgram Title: SATLAS2CPC Library link to program files:https://doi.org/10.17632/3hr8f5nkhb.2Developer's repository link:https://github.com/IKS-nm/satlas2Licensing provisions: MITProgramming language: PythonJournal reference of previous version: Computer Physics Communications, Volume 222, 2018, Pages 286-294, ISSN 0010-4655, https://doi.org/10.1016/j.cpc.2017.09.012Does the new version supersede the previous version?: YesReasons for the new version: Improved and more stable performanceSummary of revisions: The architecture of the SATLAS package has been completely changed. Instead of operating on Models, using SumModel and LinkedModel for various analysis options, SATLAS2 works with a central Fitter object to which Sources get assigned, which themselves get assigned Models. This structure improves the performance by streamlining the code and benefits from using a pass by reference approach for changing the parameters, rather than a pass by value. Most of the functionality is maintained, with a small interface for most basic things available so SATLAS2 can be used as a drop-in replacement.Nature of problem: The analysis of specifically counting data has some special considerations compared to more common datasets with Gaussian-distributed uncertainties. Application of Bayesian inference through random walk exploration is useful for more accurate exploration of parameter space. The fitting of multiple models with shared parameters can be desirable for either fitting multiple datasets with the same parameters or imposing additional restrictions on the parameters.Solution method: SATLAS2 implements the correct statistical costs for fitting of both counting data and Gaussian distributed uncertainties, and allows an easy implementation of custom cost functions.. Through the implementation of multiple Sources and Models, a natural extension of the fitting is made, so multiple models can be both summed together to fit to the same data and multiple fits can be linked together through parameter linking, using the underlying LMFIT [2] library. The ability to fit simultaneously increases the possibilities for extensive statistical analyses, and the interface through LMFIT with the emcee [3] library enables Bayesian exploration of the parameter space.

Multiscale and Multiphysical Numerical Simulations of Solid Oxide Cell (SOC)

Electrochemical applications play a key role for the topic of “green hydrogen” for the de-carbonization of the energy and mobility sectors. Electrochemical systems and processes, including fuel cells and electrolysers, have witnessed several benefits over conventional combustion-based technologies currently being widely used in power plants and vehicles. The ceramic high-temperature technologies by means of SOC exhibits high efficiencies, with a thermoelectric conversion efficiency as high as 60% and a total efficiency of up to 90% in fuel cell operation and even higher in electrolysis mode. The SOC technology, therefore, sees a promising future in the production of green hydrogen and electricity.Providing high operating temperature, over 600 oC, the SOC system shows the capability to operate with diverse types of gas mixtures, for example, hydrogen, ammonia, and carbon-containing mixtures such as methane (CH4), carbon monoxide (CO) in fuel cell operation (SOFC) and steam and/or carbon dioxide (CO2) in electrolysis operation (SOEC). The design of the SOC stack enables a reversible operation (rSOC) between fuel cell and electrolysis modes. It indicates the SOC system can perform with high efficiencies in both operating modes, which also widens the scope of possible applications.Challenges remain when it comes to commercialization of the SOC technology, in both the investment costs (CAPEX) and operating costs (OPEX) aspects. From the technological and scientific point of view, the physical transport phenomena in SOCs need to be understood which can be done by the help of experimental and numerical investigations. Cheaper and long-lasting material alternatives may be found for the cell, stack and system development afterwards.Detailed experimental investigations usually require a lot of effort in time and data analysis. To promote the scientific and technological studies on the SOC technology, numerical investigations by using multiscale and multiphysical models are carried out in this work. The models include an in-house designed/written phase field model (PFM) 1, and an open-source based computational fluid dynamics (CFD) model, openFuelCell2 2 (based on OpenFOAM). The former accounts for the microstructure evolution of the Ni/YSZ composition. The latter addresses the multiphysical transport processes in different phases, i.e., ionic transfer in YSZ, electronic transfer in Ni, and the gas diffusion in the gas phase.The figure below shows the computational domain and different phases in the numerical simulations. The evolution of Ni/YSZ composition can be predicted by the PFM. It is supposed to reproduce the Ni agglomeration that has been observed in SOFC long-term experiments 3. The simulation result shown at the left-most is obtained by performing the PFM simulation (with 96 x 96 x 96 voxels) representing the microstructure change for a certain time duration. The computational domain (192 x 192 x 192 voxels) consists of three phases, namely, YSZ, Ni, and gas, as shown in the middle and on the right side. It is refined in each direction to better capture the triple phase regions (lower right) shared by the three phases, which refers to the active sites that enable the electrochemical reaction to be conducted. By applying different governing equations on these phases, the CFD model, openFuelCell2, can describe the transport phenomena numerically. Hence, the performance degradation of a SOFC due to Ni agglomeration can be captured by carrying out the simulations for different time durations. The effective properties may be derived as well so that they can be used in numerical simulations with larger scales.AcknowledgementThe authors would like to thank their colleagues at Forschungszentrum Jülich GmbH for their great support and the Helmholtz Society, the German Federal Ministry of Education and Research as well as the Ministry of Culture and Science of the Federal State of North Rhine-Westphalia for financing these activities as part of the Living Lab Energy Campus.References Q. Li, L. Liang, K. Gerdes, and L.-Q. Chen, Appl. Phys. Lett., 101, 033909 (2012). S. Zhang, S. Hess, H. Marschall, U. Reimer, S. B. Beale, and W. Lehnert, Computer Physics Communications, to be submitted (2023). C. E. Frey, Q. Fang, D. Sebold, L. Blum, and N. H. Menzler, J. Electrochem. Soc., 165, F357 (2018). Figure 1

A C++ code for simulations of RHEED intensity oscillations within the kinematical approximation

This paper proposes a new software implementation for the rapid development and deployment of computer simulations of reflection high-energy electron diffraction (RHEED) intensity oscillations during the growth of thin epitaxial films. The calculations are based on the kinematical approximation. In this model the RHEED intensity oscillations from growing layers is considered as an effect of changes of surface roughness. New version program summaryProgram Title: growth22CPC Library link to program files:https://doi.org/10.17632/zkpwhrrfng.1Licensing provisions: GNU General Public License 3Programming language: C++Journal reference of previous version: Computer Physics Communications 170 (2005) 265–286Does the new version supersede the previous version?: YesReasons for the new version:The previous version of the program was written in Object Pascal [1]. This programming language is no longer popular in the scientific community. Responding to user's feedback, the current distribution of the software has been implemented in C++. This version implements improvements for ergonomics, computational performances and code readability. In addition, the output generated by the growth22 (saved in the coverageProfiles.dat file) can be used directly by program [2].Summary of revisions:The growth22 is written in standard C++, therefore one can use the GNU g++ compiler or any ANSI C++ compliant compiler to build the executable program. The code is distributed in the form of one growth22.cpp source file. A Table 1 lists the retired Object Pascal procedures and shows replacements by C++ functions. The present version of the program enables the user to easily manage the input data. Input data are loaded from a manually prepared inputData.dat text file where the user provides information on the parameters of the simulation, i.e. the growth model, the number of growing layers, the upper limit of growth time interval, the values of the phenomenological parameters that measure the net rate of transfer from one layer to the next, and the values of the growth rate for each layer. The user can edit this file as per the requirement of the problem domain. Details of the parameters in the sample inputData.dat file are described in the readme.pdf file provided with the software.The nature of problem:Generally, in order to explain the experimental data two theories—kinematical and dynamical—associated with high-energy electron diffraction in epitaxial layers are used. The dynamical theories take into account multiple scattering within the crystal and are used whenever diffraction from a perfect epitaxial layers is considered [3]. On the other hand, the kinematical theory treats the scattering from each plane in the crystal sample as being independent of others [1,4].In the kinematical theory of RHEED it is assumed that electron wave scatters from atom in the surface layer only once and then interferes with waves scattered on other atoms. This model can be used to describe diffraction effects for small angles of incidence, when the electron beam does not penetrate the deeper layers of the crystal. In the kinematical model the intensity oscillations of specularly reflected electron beam from growing layer is considered as an effect of changes of surface roughness.In this paper we present a time-efficient computer model for kinematical calculations of RHEED intensity oscillations during the growth of thin epitaxial films.The method of solution:According to kinematical approximation the diffracted-beam intensity is given by [4]:(1)I(S‖→=0,Sz)=|∑n=0∞(θn−θn+1)exp(−iSznd)|2 where S→ is the scattering vector, d is the interlayer spacing and n is an integer. θn−θn+1 is the exposed layer coverage which is treated as the atomic distribution function gives the position of surface atoms. For Sz=π/d Eq. (1) simplifies to:(2)I=|∑n=0∞(θn−θn+1)cos⁡(nπ)|2.Variation of the coverage ratio θn of the nth surface layer is modelled by a set of non-linear differential equations [1,4]. The Cash-Karp fifth-order Runge-Kutta method [5] with error estimation and adaptive stepsize control was used for solving initial value problem for non-linear differential equations [6]. The absolute and relative error tolerances and number of integration intervals are fully user-configurable.The proposed code approaches the modelling of RHEED from two different points of view, i.e. the data and functional perspectives. The data perspective (i.e. calculating the layer coverages as a function of their growth time) is used to characterize the data sources of the program whereas the functional perspective includes services provided by the program. It should be noted that the perspective of the data can be expressed by a completely different model introduced by the user.The coverage evolution curves for the diffusive growth model, corresponding kinematical diffracted intensity and total roughness of the surface are shown in Fig. 1. During the numerical calculations showed in Fig. 1 it was assumed that the number of growing layers (numLayers) equals 15, the upper limit of growth time interval (tMax) equals 12, the growth rate equals 1 MLs−1 (for all layers), and the values of the parameter kn (that determines the amount of interlayer diffusion) changing according to the linear relation:(3)kn=k1−(k1−kN)n−1N−1, where N is the total number of growing layers.

GenASiS Mathematics: Object-oriented manifolds, operations, and solvers for large-scale physics simulations (version 2)

GenASiSMathematics provides modern Fortran classes furnishing extensible object-oriented functionality for the solution of fields governed by selected partial differential equations. The initial release included extensible object-oriented implementations of simple meshes and the evolution of generic conserved currents thereon. This revision—Version 2 of Mathematics—includes significant reorganization and streamlining of these classes, higher-order reconstruction by a different method, a Poisson solver, coarsening to avoid Courant time step limitations near coordinate singularities, and the offloading of computational kernels to GPUs. New version program summaryProgram Title:HomogeneousSpheroid, RectangularContraction, RectangularExpansion, SphericalContraction, SphericalExpansion (example problems illustrating GenASiSMathematics)CPC Library link to program files:https://doi.org/10.17632/mzvxngwtw6.2Developer's repository link:https://github.com/GenASiSCode Ocean capsule:https://codeocean.com/capsule/4106782Licensing provisions: GPLv3Programming language: Modern Fortran; OpenMP (tested with recent versions of GNU Compiler Collection (GCC), Cray Compiler Environment (CCE), IBM XL Fortran compiler)Journal reference of previous version: Computer Physics Communications 222 (2018) 384Does the new version supersede the previous version?: YesReasons for the new version: This version includes significant reorganization and streamlining, higher-order reconstruction by a different method, a Poisson solver, coarsening to avoid Courant time step limitations near coordinate singularities, and the offloading of computational kernels to GPUs.Summary of revisions: The left part of Fig. 1 shows the revised structure of GenASiSMathematics. The Algebra and Calculus divisions contain relatively simple stand-alone functionality for such tasks as root finding, interpolation, numerical integration, and the solution of ordinary differential equations. The major functionality currently provided in Mathematics, which focuses on the solution of selected classes of partial differential equations, is now collected in the CauchyProblems division.The structure of CauchyProblems is displayed in the right part of Fig. 1. The solution of partial differential equations presumes the existence of Manifolds (or ‘spaces’) and Fields thereon whose configuration we seek. The two basic types of Cauchy problems are boundary value problems, or Constraints; and initial value problems, or Evolutions.The structures of Manifolds, Fields, Constraints, and Evolutions are displayed in Fig. 2. Individual coordinate patches or Charts are combined into Atlases in order to represent Manifolds. As for Fields, the most basic functionality for a set of fields on a manifold (including I/O) is representedby classes in FieldSets. Specialized sets of fields needed for Cauchy problems are Geometries, including coarsening functionality to avoid Courant time step limitations near coordinate singularities; and CurrentSets, which represent ‘conserved’—or more properly, when source terms are present, ‘balanced’—currents, that is, densities and their corresponding fluxes. CalculusFields provides integrals of fields on manifolds, including volume and surface integrals used to tally conserved and balanced quantities, as well as spherical and azimuthal averaging. New to this version of Mathematics is the solution of PoissonEquations via multipole expansion, the first capability for Constraints or boundary value problems.The solution of balance equations is the type of initial value problem currently handled by Evolutions. As in Version 1, finite-volume discretization results in a large set of ordinary differential equations solved by the Runge-Kutta method. Now however, the right-hand sides, or Slopes, are handled by variations of a recursive class that allows the construction of hierarchies of terms. Thus complicated equations can be computed and included in I/O in pieces representing various terms, with results assembled in an extensible and quasi-automated fashion. For example, the partial derivatives, and geometric source terms from connection coefficients, can be computed and visualized separately, while being automatically combined to give the full divergence of an energy-momentum tensor. Third-order reconstruction (with limiting) to discretize cell interfaces is now included. Multi-stage Runge-Kutta Steps are then taken by Integrators in order to advance the solution of fields from their initial condition.Relative to Version 1 of Mathematics, the Manifolds and Fields classes are significantly reorganized and streamlined. Previously, classes for fields, geometry, and I/O were interleaved with those of charts and atlases. This provided for chart and atlas classes with self-contained geometry and I/O functionality, but in a manner that proved unwieldy and indeed unnecessary. This has now been disentangled: classes for fields, geometry and I/O are built entirely on top of full manifolds (i.e. classes for atlases).All computational kernels (except those specifically needed for checkpointing, that is, I/O) are coded for offloading to GPUs via OpenMP.Nature of problem: By way of illustrating GenASiSMathematics functionality, solve an example potential problem and imposed advection problems.Solution method: Multipole expansion for solution of the Poisson equation; finite-volume discretization; second- and third-order reconstruction with limiting; HLL Riemann Solver; Runge-Kutta integration.Additional comments including restrictions and unusual features: Uses the MPI [1] and Silo [2] libraries. The example problems named above are not ends in themselves, but serve to illustrate the functionality available though GenASiSMathematics. In addition to these more substantial examples, we provide individual unit test programs for the classes comprised by GenASiSMathematics.GenASiSMathematics is available in the CPC Program Library and also at https://github.com/GenASiS.

PorousMultiphaseFoam v2107: An open-source tool for modeling saturated/unsaturated water flows and solute transfers at watershed scale

The first releases of porousMultiphaseFoam proposed an open-source software suite to solve the equations for multiphase flow (generalized Darcy’s law) in porous media or groundwater flows (Richards' equation) by taking advantage of OpenFOAM, a finite volume platform with automatic discretization on three-dimensional unstructured grids and good parallel efficiency. Recently, the porousMultiphaseFoam toolbox has been confronted with complex cases of fast water flows and solute transfers in realistic hydrological configurations with variable forcing conditions (heterogeneous infiltration and local tracer injection). Several developments have been carried out to make it possible to simulate those cases, which extend the toolbox with: (i) a set of solvers dedicated to groundwater flows, including coupled water flow and solute transport and simplified 2D approaches, (ii) improved numerical techniques for problems with strong non-linearities, (iii) libraries/executables for pre-processing of input data (geographical information and time-variable forcing terms) and (iv) passive or coupled scalar transport (tracer) with groundwater solvers that support any number of species. New solvers are validated on several (un-)saturated configurations by a direct comparison with a well validated finite element code. Program summaryProgram title: porousMultiphaseFoam v2107CPC Library link to program files:https://doi.org/10.17632/hphn58ksfy.1Developer's repository link:https://github.com/phorgue/porousMultiphaseFoam.gitLicensing provisions: GNU GPLv3Programming language: C++ / OpenFOAMJournal reference of previous version: Computer Physics Communications 187 (2015) 217-226Does the new version supersede the previous version?: YesReasons for the new version: Multiple new features of the toolbox dedicated to groundwater flows (with pre-processing tools)Summary of revisions:–Solvers:–groundwaterFoam/stationaryGroundwaterFoam: new unsaturated flow (transient/stationary) solvers in porous media (Richards' equation) using Picard or Newton's method for linearization.–groundwater2DFoam/stationaryGroundwater2DFoam: new 2D (transient/stationary) solver for saturated flow in porous media.–porousScalarTransportFoam/porousScalarTransport2DFoam: new passive multi-scalar transport solvers on pre-computed flow fields (3D unsaturated or 2D saturated fields).–groundwaterTransportFoam/groundwaterTransport2DFoam: new solvers coupling flow modeling (2D saturated or 3D unsaturated) with multi-scalar tracers transports.–Improved numerical techniques–Timestep management: time-stepping is now based on the temporal truncation of the time scheme used (backward Euler 1st and 2nd order, Crank-Nicolson).–Linearization techniques: problems of unsaturated flow with strong non-linearities can now be solved using an algorithm based on Newton's method.–Libraries:–numericalMethods: the EulerD3dt3Scheme class allows computation of the 3rd-order time derivative to estimate the truncation error and the JacobianMatrix class computes the Jacobian using finite differencing (expensive method used only for debugging).–toolsGIS: new dynamic library for managing GIS (Geographical Information System) files and event files. Events are variable and localized source terms for fluid or scalar transport.–Ippisch: an improved two-phase model for relative permeability and capillary pressure.–porousBoundaryConditions: several new boundary conditions have been added to handle heterogeneous pressure head boundary conditions and variable fluid or tracer injections (using event files).–Utilities:–setBoundaryHeadPressure: pressure head initializer for boundary conditions using uniform value of GIS files.–setFieldsFromDEM/setFieldsFromXY: heterogeneous initializer for internal fields, similar to the classical OpenFOAM setFields, but using ordered (DEM) or unordered (XY) grid files.–darcyFoam: simple utility to calculate a single-phase flow in porous media.Nature of problem: The porousMultiphaseFoam toolbox was initially released with generic two-phase solvers for flow in porous media. Developments were carried out to model more specifically hydrogeological flows by including a groundwaterFoam solver for Richards' equation. However, many features required to simulate realistic configurations, such as event management for variable and heterogeneous forcing terms (infiltration / scalar source term) or suitable numerical techniques for complex configurations (high rate injection for unsaturated flow for example), were still missing. Moreover, the 2D approach with integration of vertical properties (for flow and tracer) was not present even though this is an indispensable tool for groundwater flow modeling due to its ease of implementation and its low computation cost.Solution method: Improved numerical techniques have been developed to solve realistic flows by handling strong non-linearities with Newton's method, and a wide range of time scales with time-stepping based on time truncation errors. Various solvers for groundwater flows have been developed to handle cases under different assumptions, including the aforementioned 2D approaches. Variants of these solvers were also developed to simultaneously solve the passive transport of any number of scalars (tracers) by the flow. A dynamic library has been developed to pre-process input data: event information for variable conditions, GIS file for spatial information.Additional comments including restrictions and unusual features: porousMultiphaseFoam v2107 requires an installation of OpenFOAM-v9. Dedicated branches are available on the developer's repository for compatibility with other versions of OpenFOAM. The following versions are compatible with the current version of porousMultiphaseFoam: v7, v8, v1906, v2006 and v2106. Older branches may be found in the repository link but do not contain all features.

Medical module for Hyper-Fractal Analysis

In Grossu (2022) [1] it was discussed the migration of Hyper-Fractal Analysis from Visual Basic 6 to C# .Net. The main goal of current work was developing a medical module for fractal analysis of computed tomography multi-channel images. A new tool for comparing images by RGB channels superposition was also considered. This could be of particular interest in medical image analysis (e.g. compare native with contrast CT slices). New version program summaryProgram Title: Hyper-Fractal Analysis v06CPC Library link to program files:https://doi.org/10.17632/z9knmny56p.2Licensing provisions: GPLv2Programming language: C# 7.3 /.Net Framework 4.7.1Journal reference of previous version: Computer Physics Communications 271, (2022) 108189Does the new version supersede the previous version?: YesNature of problem: Estimating the fractal dimension of images, multi-dimensional objects, and medical computed tomography images.Solution method: Optimized version of the fuzzy box-counting algorithm.Reasons for the new version: Develop a medical module for fractal analysis of computed tomography multi-channel images.Summary of revisions:•The migration of Hyper-Fractal Analysis from Visual Basic 6 to C# .Net was discussed in [1]. The new SOLID [2] application design is facilitating the development process. Thus, the main goal of current version was adding a medical module for fractal analysis [3-5] of computed tomography (CT) images [6,7].•Hyper-Fractal Analysis v06 allows opening specific csv files (.ct.csv) containing up to three Hounsfield units matrices (CT slices) separated by empty lines. Most commonly, those three “channels” were used for storing the corresponding slices from native, arterial contrast, and venous contrast series.•Each of the previously described CT slices is represented on a different RGB channel [8], which results in highlighting the images differences (Fig. 1).•A new tab (CT Band-Pass) was added to allow specifying the Hounsfield units input parameters for the CT multi-channel fuzzy band-pass filter.

A Deep Neural Network Based on ResNet for Predicting Solutions of Poisson–Boltzmann Equation

The Poisson–Boltzmann equation (PBE) arises in various disciplines including biophysics, electrochemistry, and colloid chemistry, leading to the need for efficient and accurate simulations of PBE. However, most of the finite difference/element methods developed so far are rather complicated to implement. In this study, we develop a ResNet-based artificial neural network (ANN) to predict solutions of PBE. Our networks are robust with respect to the locations of charges and shapes of solvent–solute interfaces. To generate train and test sets, we have solved PBE using immersed finite element method (IFEM) proposed in (Kwon, I.; Kwak, D. Y. Discontinuous bubble immersed finite element method for Poisson–Boltzmann equation. Communications in Computational Physics 2019, 25, pp. 928–946). Once the proposed ANNs are trained, one can predict solutions of PBE in almost real time by a simple substitution of information of charges/interfaces into the networks. Thus, our algorithms can be used effectively in various biomolecular simulations including ion-channeling simulations and calculations of diffusion-controlled enzyme reaction rate. The performance of the ANN is reported in the result section. The comparison between IFEM-generated solutions and network-generated solutions shows that root mean squared error are below 5·10−7. Additionally, blow-ups of electrostatic potentials near the singular charge region and abrupt decreases near the interfaces are represented in a reasonable way.

Migration of hyper-fractal analysis from visual basic 6 to C# .Net

In Grossu et al. (Comput. Phys. Commun. 184 (2013) 1812–1813) we presented Hyper-Fractal Analysis, a visual tool for estimating the fuzzy fractal dimension of images and 4D objects. As Visual Basic 6 could be considered an outdated language, with limited Object-Oriented Programming capabilities, migrating the application to C# .Net was treated in high priority. Following the goal of creating a highly reusable fractal analysis library, the code was also refactored to SOLID. Together with various improvements, the.Net version is also providing new tools for iso-fractal areas identification. The project success was confirmed by a comparative old/new version study. On the other hand, the most relevant functionalities were covered by unit tests. New version program summaryProgram Title: Hyper-Fractal Analysis v05CPC Library link to program files:https://doi.org/10.17632/z9knmny56p.1Licensing provisions: GPLv2Programming language: C# 7.3 /.Net Framework 4.7.1Journal reference of previous version: Computer Physics Communications 184, Issue 7 (2013) 1812–1813Does the new version supersede the previous version?: YesNature of problem: Estimating the fractal dimension of images and multi-dimensional objects.Solution method: Optimized version of the fuzzy box-counting algorithm.Reasons for the new version: Migration from VB6 to C# .Net, adding new tools for iso-fractal areas identification.Summary of revisions:•The main goal of this work was creating a highly reusable fractal analysis library. Following this purpose Hyper-Fractal Analysis v04 [1] was migrated from VB6 (an outdated language, with limited OOP capabilities) to C# .Net [2]. This opportunity was also used for refactoring the code to SOLID [3].•The box-counting algorithm [4] implemented in the previous version is limited to 4D objects. In current version the number of dimensions is limited only by the available resources.•The previous fuzzy [5,6] box-counting implementation is accepting only 0-255 shades of gray windows. The new version allows defining parametrizable fuzzy band-pass filters.•An important attention was also paid on developing new tools for identification of iso-fractal areas (the sets of points characterized by the same local fractal dimension [7]). This functionality is expected to be of particular interest in medical image analysis [7-13].•The Windows Forms GUI (Fig. 1) was completely reconsidered on top of the Model View Presenter pattern [14].•For an increased reliability and for minimizing the regression risks during future developments the most relevant functionalities were covered by unit tests [15].

ElectroMicroTransport v2107: Open-source toolbox for paper-based electromigrative separations

electroMicroTransport v2107: Open-source toolbox for paper-based electromigrative separations

Light Intensity Modulated Photo-Electrochemical Methods and Distribution of Relaxation Times As Tools to Investigation Photovoltaic Materials

Electrical Impedance Spectroscopy (EIS, sometimes also Electrochemical Impedance Spectroscopy) found its place in material studies as a routine experiment. EIS proliferation became possible through advances in electronics. Current-voltage characteristics were standard features of potentiostats. Impedance techniques require additional components embedded into a potentiostat, namely signal generator and signal analyzer. Relatively simple devices are available, although dedicated instrumentation with little noise, broad frequency control and perhaps sample temperature control or automated sample handling are particularly desired and command premium in cost.EIS as an electrical technique uses as the perturbing signal alternating electrical signal, either potential of current. In general, impedance studies are intended to study a transfer function and the input and output do not have to be only voltage-current or current-voltage as it would be in EIS. In mechanical engineering are for example useful studies focused on mechanical perturbation. Electrical response by its nature makes EIS amenable to study electrochemical systems.Along with the impedance vs. frequency data collection comes also the need of meaningful data interpretation. This is typically done through data fitting to components of an equivalent electrical circuit. The acquired data -- real and imaginary impedance components gathered over a vast range of frequencies -- are evaluated by fitting to some equivalent circuit using nonlinear regression.Oftentimes both finding the equivalent circuit and the physical component assignment are tedious, if not impossible. Recently a new approach to data reduction has been a treatment via distribution of relaxation times (DRT) [1, 2]. In principle, this approach has been outlined already in 1936 [3]. Its present use was made possible by available computation power of current PCs [4, 5]. While this is not a substitute to equivalent circuits and subsequent physical data interpretation, it is possibly a method, when used appropriately and in automated way, that can provide more information that simple comparison of the Bode or Nyquist graphs. This note is presenting some examples of this approach.In the field of photovoltaic studies of light-gathering materials experimentation on systems responding to illumination makes use of impedance studies in which the perturbing signal is light of certain wavelength (or white light) with variable intensity. Here we will look at the link between EIS and electrical response spectroscopy induced by fluctuating light. This is sometimes for short dubbed as optical spectroscopy although it does make a use of the instruments for optical spectroscopy (such as UV-VIS, IR etc.) as they are known from analytical chemistry literature. Rather, an electric signal in a photoactive material is generated by intensity modulated light at different frequencies (similar range to electrical impedance spectroscopy) and a response (the magnitude and the phase shift) is analyzed. In this respect the magnitude and phase shift correspond to the magnitude and a vector of resistance on the polar plot of electrical impedance measurement. However, the data obtained serve a different purpose. For example, in a photovoltaic shorted system, the time constant associated with the response corresponds to mobility of the photogenerated charge in the material. For an open-circuit potential of such a photovoltaic system the measured time constant response corresponds to the rate of recombination of the generated hole-electron species.Such materials of general structure (ABX3) could be replacements for silicon-based photovoltaic cells. These materials, due to the possible variability of both cations and anions in their composition, offer almost endless possibilities for synthetic preparation. For the active layer of a solar cell is often used the CH3NH3PbI3-2Cl2 perovskite structure, i.e., a structure containing an organic component. We have built and then tested such photovoltaic cells by static and dynamic methods using a unique photoelectric research instrument using the method CIMPS (Controlled Intensity Modulated Photocurrent Spectroscopy).AcknowledgementThis work was supported through the support of The Czech Science Foundation (GACR), contract GA19-23718S.References BioLogic, Distribution of Relaxation Times (DRT): an introduction, BioLogic, Editor. 2017.Weese, J., A reliable and fast method for the solution of Fredhol integral equations of the first kind based on Tikhonov regularization. Computer Physics Communications, 1992. 69(1): p. 99-111.Yager, W.A., The Distribution of Relaxation Times in Typical Dielectrics. Physics, 1936. 7(12): p. 434-450.Shoar Abouzari, M.R., et al., On the physical interpretation of constant phase elements. Solid State Ionics, 2009. 180(14): p. 922-927.Wan, T.H., et al., Influence of the Discretization Methods on the Distribution of Relaxation Times Deconvolution: Implementing Radial Basis Functions with DRTtools. Electrochimica Acta, 2015. 184: p. 483-499.

JaSTA-3: Light scattering simulations for heterogeneous aggregate

JaSTA-3: Light scattering simulations for heterogeneous aggregate

HTR-1.2 solver: Hypersonic Task-based Research solver version 1.2

We present an updated version of the open-source Hypersonics Task-based Research (HTR) solver for hypersonic aerothermodynamics. The solver, whose first version was presented in Di Renzo et al. (2020), is designed for direct numerical simulation (DNS) of canonical hypersonic flows at high Reynolds numbers in which thermo-chemical effects induced by high temperatures are relevant. The solver relies on high-order spatial discretization on structured meshes and efficient time integrators for stiff systems within the Regent/Legion software stack, which makes the code highly portable and scalable in CPU and GPU-based supercomputers. The new version herein presented includes several optimizations and new tools for data analysis, along with novel user option for hybrid skew-symmetric/targeted essentially non-oscillatory numerics, to offer higher computational efficiency and lower numerical dissipation at moderate Mach numbers, inclusion of a new combustion mechanism for methane and oxygen, and new recycling–rescaling inlet boundary conditions targeted to the simulation of fully developed turbulent boundary layers. New version program summaryProgram Title: Hypersonics Task-based Research solverCPC Library link to program files:https://doi.org/10.17632/9zsxjtzfr7.2Developer’s repository link:https://github.com/stanfordhpccenter/HTR-solver.gitLicensing provisions: BSD 2-clauseProgramming language: Regent, C++, and CUDAJournal Reference of previous version: Di Renzo, M., Fu, L., & Urzay, J. (2020). HTR solver: An open-source exascale-oriented task-based multi-GPU high-order code for hypersonic aerothermodynamics. Computer Physics Communications, 107262.Does the new version supersede the previous version?: YesReasons for the new version: Release of new featuresSummary of revisions:•New optional sixth-order hybrid scheme has been implemented (activated by the flag “hybridScheme” in the input file). The scheme combines the energy-preserving properties of a sixth-order skew-symmetric central difference scheme [1] in smooth flow regions with the shock-capturing properties of a sixth-order targeted essentially non-oscillatory (TENO) scheme at points where shocks are involved. The switch between the two schemes is controlled by a TENO sensor whose cutoff value is adapted based on the maximum value of a modified Ducros sensor [2] across the reconstruction stencil. If the flag “hybridScheme” is set to false, the numerical scheme will revert to the TENO6-A scheme released in the previous version of the solver [3];•New recycling–rescaling inflow boundary conditions [4,5] for the simulation of turbulent compressible boundary layers are now available to the user;•Support for Legion tracing, which significantly improves the strong scalability of the solver, has been implemented;•A diagnostic tool to monitor the time evolution of the flow variables in a subvolume of the computational domain is now available;•A single-step chemistry mechanism for methane/oxygen combustion has been added to the mixtures handled by the HTR solver;•Sample scripts for strong scaling have been added to the “testcases” directory;•Unit test and regression test suites have been added to the repository;•The input file scheme has been modified in order to reduce verbosity and increase flexibility in specifying the boundary conditions and type of gas mixture;•Hyperbolic sine stretching functions have been made available to users during the grid generation process;•Computationally intensive tasks have been ported to C++ and CUDA in order to achieve higher efficiency on all hardware;•Several optimizations of the tasks body and mapper have been implemented in order to increase the computational efficiency and reduce the memory footprint.Nature of problem: This code solves the Navier–Stokes equations at hypersonic Mach numbers including finite-rate chemistry for dissociating air and multicomponent transport. The solver is designed for direct numerical simulations (DNS) of transitional and turbulent hypersonic turbulent flows under high-enthalpy conditions, and it accounts for thermochemical effects (vibrational excitation and chemical dissociation).Solution method: This code uses a low-dissipation sixth-order schemes for the spatial discretization of the conservation equations on Cartesian stretched meshes. Time advancement is carried out by either an explicit method if chemistry is slow, hence not introducing additional stiffness, or by an operator-splitting algorithm whereby chemical production rates are handled implicitly.Additional comments including restrictions and unusual features: The HTR solver builds on the runtime Legion [6,7] and is written in the programming language Regent [8,9] developed at Stanford University. Instructions for installation of the components are provided in the README file enclosed with the HTR solver and in the Legion repository [6].

Momentum-space-resolved measurements using oblique electron cyclotron emission for the validation of the quasi-linear theory of electron cyclotron current drive

Electron cyclotron resonance heating (ECRH) can drive large current densities through electron cyclotron current drive (ECCD). ECCD is expected to be crucial for high-performance plasmas in future fusion reactors like ITER and DEMO, making the current drive efficiency of ECCD a critical design parameter for future reactors. In present-day devices, good agreement between measured and predicted current drive efficiency has been found. However, to ensure the reliability in future machines, a direct validation of the electron momentum distribution function is needed.As a first step towards this goal, we present in this paper oblique electron cyclotron emission (ECE) measurements of a low-density plasma in the ASDEX Upgrade tokamak. Two oblique ECE diagnostics are used to allow the simultaneous measurements of electrons streaming co- and counter-directionally with the plasma current. Predictions for the distribution function are computed with the bounce-averaged Fokker-Planck code RELAX. 9 9E Westerhof et al, Rijnhuizen report,1992. To allow direct comparison with the measurements, synthetic radiation temperatures are computed with the code ECRad. 10 10S Denk et al, Computer Physics Communications, p. 107175, 2020. Good agreement is found if radial transport occurring predominantly at low electron energies is included.We demonstrate that oblique ECE diagnostics measure the electron distribution function directly at the ECRH deposition site in phase space. Furthermore, they are sensitive to the abundance of pitch-angle scattered electrons that reduce the ECCD efficiency. Limitations and uncertainties of the measurements and the modeling are discussed.

An effective method to calculate RHEED rocking curves from nanoheteroepitaxial systems

An effective method to calculate RHEED rocking curves from nanoheteroepitaxial systems

ADG: Automated generation and evaluation of many-body diagrams II. Particle-number projected Bogoliubov many-body perturbation theory

We describe the second version (v2.0.0) of the code ADG that automatically (1) generates all valid off-diagonal Bogoliubov many-body perturbation theory diagrams at play in particle-number projected Bogoliubov many-body perturbation theory (PNP-BMBPT) and (2) evaluates their algebraic expression to be implemented for numerical applications. This is achieved at any perturbative order p for a Hamiltonian containing both two-body (four-legs) and three-body (six-legs) interactions (vertices). All valid off-diagonal BMBPT diagrams of order p are systematically generated from the set of diagonal, i.e.,unprojected, BMBPT diagrams. The production of the latter were described at length in Arthuis et al. (2019) dealing with the first version of ADG. The automated evaluation of off-diagonal BMBPT diagrams relies both on the application of algebraic Feynman’s rules and on the identification of a powerful diagrammatic rule providing the result of the remaining p-tuple time integral. The new diagrammatic rule generalizes the one already identified in Arthuis et al. (2019) to evaluate diagonal BMBPT diagrams independently of their perturbative order and topology. The code ADG is written in Python3 and uses the graph manipulation package NetworkX. The code is kept flexible enough to be further expanded throughout the years to tackle the diagrammatics at play in various many-body formalisms that already exist or are yet to be formulated. Program summaryProgram Title:ADGCPC Library link to program files:https://doi.org/10.17632/6h4xrydwfb.2Licensing provisions: GPLv3Programming language: Python3Journal reference of previous version: P. Arthuis, T. Duguet, A. Tichai, R.-D. Lasseri and J.-P. Ebran, ”ADG: Automated generation and evaluation of many-body diagrams I. Bogoliubov many-body perturbation theory”, Computer Physics Communications 240 (2019), pp. 202-227.Does the new version supersede the previous version?: Yes.Reasons for the new version: Incorporation of a new formalism into the program.Summary of revisions: Addition of off-diagonal BMBPT to the formalisms for which diagrams can be generated, fix of a wrong symmetry factor, move of the codebase from Python2 to Python3 while maintaining support for Python2, various optimizations to reduce the time and memory necessary to the program.Nature of problem: As formal and numerical developments in many-body-perturbation-theory-based ab initio methods make higher orders reachable, manually producing and evaluating all the diagrams becomes rapidly untractable as both their number and complexity grow quickly, making it prone to mistakes and oversights.Solution method: Diagonal BMBPT diagrams are encoded as square matrices known as oriented adjacency matrices in graph theory, and then turned into graph objects using the NetworkX package. Off-diagonal BMBPT diagrams can then be generated from those graphs. Checks on the diagrams and evaluation of their time-integrated expression are eventually done on a purely diagrammatic basis. HF-MBPT diagrams are produced and evaluated as well using the same principle.

Chaos Many-Body Engine module for estimating pentaquark production in proton–proton collisions at CBM energies

In Grossu et al., (2019) we proposed a Chaos Many-Body Engine (CMBE) quark toy-model for the Compressed Baryonic Matter (CBM) energies. We started from the following assumptions: (1) the system can be decomposed into a set of two or three-body quark “elementary systems”, i.e. “white” color charged, mesonic or baryonic systems; (2) the bi-particle force is limited to the domain of each elementary system; (3) the physical solution conforms to the minimum potential energy requirement. In the present work we used graph theory for identifying those sets (clusters) of elementary systems close enough to form a bound system (through the exchange of same color charged quarks). In this approach, the cluster production could be understood as an effect of the chaotic behavior of the system. As a direct application, we estimated the pentaquark production probability obtained in p + p collisions, at a center-of-mass energy between 10 and 100 GeV New version program summaryProgram Title: Chaos Many-Body Engine (CMBE)CPC Library link to program files:http://dx.doi.org/10.17632/rh5txj3n4g.2Licensing provisions: GPLv2Programming language: C# 7.3.External routines: BigRational structure (Microsoft).Journal reference of previous version: Computer Physics Communications 239C (2019) 149-152Does the new version supersede the previous version?: Yes.Nature of problem: Estimate the Pentaquark production in nuclear relativistic collisions.Solution method: Clustering algorithm for identifying all five-body quark white systems.Reasons for the new version: Added the Pentaquark identification new feature.Summary of revisions:•Migration from .Net Framework 4.0 to .Net Framework 4.7.1•In [1] we implemented a quark confinement algorithm for decomposing the system into a set of two or three-body quark “elementary systems”, i.e. “white” color charged, mesonic or respectively, baryonic systems, in agreement with the minimum potential energy requirement. In this work we added a new O(n3) algorithm (QcdQuarkBagPentaQuarkAlgorithm class), developed in agreement with the SOLID principles, for the identification of those sets (clusters) of elementary systems close enough to form a bound system (through the exchange of same color charged quarks). Taking this into account the system was associated with an undirected graph G [2,3], whose nodes are the elementary systems. Two nodes are connected if the distance between their geometrical centers is lower than the sum of their radii. Thus, each cluster could be associated to a maximal connected subgraph of G.•Unit tests (QcdQuarkBagPentaQuarkAlgorithmTests class) for checking the new algorithm.•In [1] we proposed a quark toy-model for proton–proton collisions at CBM energies [4]. The model was extended for estimating the pentaquark [5] production probability, as seen in Fig. 1. In this approach, the pentaquark production could be understood as a direct effect of the chaotic behavior of the system [6,7]. Thus, for each center-of-mass energy (s∈10,100GeV) we simulated 2,000 events by choosing Simulations\\Quark Collision from the menu and storing the pentaquark multiplicity at t=300 Fm/c in the quark.pentaquark.log.csv file, generated into the simulation output folder. The collision parameter was given random values in the [0.2, 1.1] Fm range.ReferencesI.V. Grossu, C. Besliu, Al. Jipa, D. Felea, E. Stan, T. Esanu, Implementation of quark confinement and retarded interactions algorithms for Chaos Many-Body Engine, Computer Physics Communications 239C (2019) pp. 149-152, DOI: https://doi.org/10.1016/j.cpc.2019.01.023Alan Gibbons, Algorithmic Graph Theory, Cambridge University Press, 1985.I.V. Grossu, C. Besliu, Al. Jipa, C.C. Bordeianu, D. Felea, E. Stan, T. Esanu, Code C# for chaos analysis of relativistic many-body systems, Computer Physics Communications 181 (2010) 1464–1470, https://doi.org/10.1016/j.cpc.2010.04.015A. Abuhoza et al, The CBM Collaboration, Nuclear Physics A, Volumes 904--905, 2 May 2013, Pages 1059c-1062cA. Abdivaliev, C. Besliu et al., Yad.Fiz.29, v.6, 1979.St. Grosu, Quelques effets des fluctuations de la barrière de potentiel a la surface des conducteurs, Studii si cercetari de fizica, 1, XI, 1960D. Felea, C.C. Bordeianu, I.V. Grossu, C. Besliu, Al. Jipa, A.-A. Radu and E. Stan, Intermittency route to chaos for the nuclear Billiard, EPL, 93 (2011) 42001, DOI: 10.1209∕0295−5075∕93∕42001

VarRhoTurbVOF 2: Modified OpenFOAM volume of fluid solvers with advanced turbulence modeling capability

varRhoTurbVOF contains a set of OpenFOAM volume of fluid (VOF) solvers for turbulent isothermal multiphase flows, which are variable-density incompressible. Unlike their official counterparts, where Favre-averaged and Reynolds-averaged velocities coexist in different equations, new solvers use Favre-averaged velocities consistently in all equations. This major update introduces three main improvements to the previous version of varRhoTurbVOF. First, the implementation is extended to VOF solvers for isothermal and non-isothermal phase change two-phase flows, where the flow is no longer incompressible. Second, in order to introduce backward compatibility and to avoid code duplication, the turbulence model construction procedure is redesigned such that solvers can determine whether the variable-density effect is considered or not in the turbulence modeling part based on the input file at run time. Third, the Egorov turbulence damping model for ω-based turbulence models is implemented with its most recent developments. Plus, an extension to ϵ-based turbulence models is developed and implemented. Program summaryProgram Title: varRhoTurbVOFCPC Library link to program files:http://dx.doi.org/10.17632/4t8z8vzyvs.2Developer’s repository link:https://github.com/wenyuan-fan/varRhoTurbVOF_2Licensing provisions: GPLv3Programming language: C++Journal reference of previous version: Wenyuan Fan, Henryk Anglart, 2020. varRhoTurbVOF: A new set of volume of fluid solvers for turbulent isothermal multiphase flows in OpenFOAM, Computer Physics Communications, 247, 106876.Does the new version supersede the previous version?: YesReasons for the new version: Implementation of new solvers and models; code redesign to maintain backward compatibility.Summary of revisions: The updated version mainly has three new features in comparison with the previous version: •The issue addressed in [1] exists in all OpenFOAM VOF solvers where the variable-density effect is neglected in the turbulence modeling. In order to solve this issue for phase change VOF solvers, two additional solvers, namely, varRhoInterPhaseChangeFoam and varRhoInterCondensatingEvaporatingFoam, are implemented based on their official counterparts.•The procedure of constructing the turbulence model inside individual solvers is redesigned. Using the new design, solvers will construct turbulence models with or without the variable-density effect being considered based on the user input. When the latter mode is activated, solvers will behave just like their official counterparts. In such a way, the backward compatibility of the code is achieved without introducing code duplications.•The Egorov turbulence damping model [2], which is a popular approach to phenomenologically modify turbulence behaviors near a two-phase interface, is implemented. In addition to the original model, recent modifications, e.g., a new length scale calculation and asymmetric damping treatments [3], are implemented as well. The Egorov model only works with ω-based turbulence models. In the new release, the model is extended to ϵ-based turbulence models following the idea proposed in [4] such that the turbulence damping model is applicable to commonly used turbulence models.Nature of problem: Within the VOF framework, the mixture has a density changing with phases. However, in quite a few OpenFOAM VOF solvers, the variable-density effect is considered in the momentum equation but neglected in the turbulence modeling. As a result, Favre-averaged velocities and Reynolds-averaged velocities are used in the momentum equation and turbulence quantity equations, respectively and simultaneously, introducing a severe self-inconsistency. Solutions to this problem have been proposed in [1] with no backward compatibility support. Also, no phase change solvers are supported in the previous version. Another issue with turbulent VOF simulations is that, even with the correct implementation, certain modifications are still needed to correct the turbulence behavior around the interface.Solution method: The extension to phase change VOF solvers is straightforward according to [1]. Regarding the construction of turbulence models, two additional fields, namely, rhoTurb and rhoPhiTurb, are created and used in turbulence models. If the user chooses to use the variable-density turbulence models, real two-phase density rho and mass flux rhoPhi will be assigned to rhoTurb and rhoPhiTurb, respectively. Otherwise, rhoTurb and rhoPhiTurb will use unity and volume flux phi, respectively. This run time selection design enables solvers to operate in two modes, which guarantees the backward compatibility and avoids code duplications. The implementation of the turbulence damping model also takes the backward compatibility into account. At run time, the code will automatically determine whether the mixture density should be included in the turbulence damping model by checking the dimension of the corresponding equation.Additional comments: A manual is provided to introduce how the turbulence damping modeled is formulated and how to use the model. It should be emphasized that the Egorov damping model is a phenomenological model. Model parameters should be carefully selected for any given mesh and flow condition.