From 324f780e450bd99c8f95d7dd2a98a6bb4403ad0b Mon Sep 17 00:00:00 2001 From: Malik Date: Wed, 1 Oct 2025 11:32:28 -0600 Subject: [PATCH 01/17] handle missing d.gas --- bird/postprocess/post_quantities.py | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/bird/postprocess/post_quantities.py b/bird/postprocess/post_quantities.py index 15a91515..45d00d22 100644 --- a/bird/postprocess/post_quantities.py +++ b/bird/postprocess/post_quantities.py @@ -912,6 +912,7 @@ def compute_instantaneous_kla( species_names: list[str], n_cells: int | None = None, volume_time: str | None = None, + constant_bubble_diameter: float | None = None, field_dict: dict | None = None, ) -> tuple[dict, dict, dict]: r""" @@ -979,6 +980,9 @@ def compute_instantaneous_kla( volume_time : str | None Time folder to read to get the cell volumes. If None, finds volume time automatically + constant_bubble_diameter : float | None + Constant bubble diameter value to handle missing d.gas + If None, reads d.gas field_dict : dict Dictionary of fields used to avoid rereading the same fields to calculate different quantities @@ -1048,9 +1052,12 @@ def compute_instantaneous_kla( U_liq, field_dict = read_field( field_name="U.liquid", field_dict=field_dict, **kwargs ) - d_gas, field_dict = read_field( - field_name="d.gas", field_dict=field_dict, **kwargs - ) + try: + d_gas, field_dict = read_field( + field_name="d.gas", field_dict=field_dict, **kwargs + ) + except FileNotFoundError: + d_gas = constant_diameter mu_liq, field_dict = read_field( field_name="thermo:mu.liquid", field_dict=field_dict, **kwargs ) From 28cb572485580733d280a8e99b16a8d8d92349a1 Mon Sep 17 00:00:00 2001 From: Malik Date: Wed, 1 Oct 2025 11:34:35 -0600 Subject: [PATCH 02/17] specify units --- bird/postprocess/post_quantities.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/bird/postprocess/post_quantities.py b/bird/postprocess/post_quantities.py index 45d00d22..65d5876a 100644 --- a/bird/postprocess/post_quantities.py +++ b/bird/postprocess/post_quantities.py @@ -981,7 +981,7 @@ def compute_instantaneous_kla( Time folder to read to get the cell volumes. If None, finds volume time automatically constant_bubble_diameter : float | None - Constant bubble diameter value to handle missing d.gas + Constant bubble diameter value (in m) to handle missing d.gas If None, reads d.gas field_dict : dict Dictionary of fields used to avoid rereading the same fields to calculate different quantities From bba709368d06f70f049615c84d7dc82316ffff28 Mon Sep 17 00:00:00 2001 From: Malik Date: Wed, 1 Oct 2025 11:35:33 -0600 Subject: [PATCH 03/17] fix argument passed --- bird/postprocess/post_quantities.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/bird/postprocess/post_quantities.py b/bird/postprocess/post_quantities.py index 65d5876a..d9a2243b 100644 --- a/bird/postprocess/post_quantities.py +++ b/bird/postprocess/post_quantities.py @@ -1057,7 +1057,7 @@ def compute_instantaneous_kla( field_name="d.gas", field_dict=field_dict, **kwargs ) except FileNotFoundError: - d_gas = constant_diameter + d_gas = constant_bubble_diameter mu_liq, field_dict = read_field( field_name="thermo:mu.liquid", field_dict=field_dict, **kwargs ) From 6b0b1288c672173c1fd43f6b9f32b8ecf59de534 Mon Sep 17 00:00:00 2001 From: Malik Date: Wed, 1 Oct 2025 16:53:44 -0600 Subject: [PATCH 04/17] force outputting of mu.liquid --- .../bubble_column_pbe_20L/system/controlDict | 6 ++++++ .../FlatPanel_250L_ASU/system/controlDict | 7 +++++++ tutorial_cases/airlift_40m/system/controlDict | 16 ++++++++++++++++ .../bubble_column_20L/system/controlDict | 2 ++ .../loop_reactor_mixing/system/controlDict | 7 +++---- .../loop_reactor_reacting/system/controlDict | 12 ++++++------ tutorial_cases/stirred_tank/system/controlDict | 11 +++++++++++ 7 files changed, 51 insertions(+), 10 deletions(-) diff --git a/experimental_cases/disengagement/bubble_column_pbe_20L/system/controlDict b/experimental_cases/disengagement/bubble_column_pbe_20L/system/controlDict index 0f3246f2..64dfc486 100644 --- a/experimental_cases/disengagement/bubble_column_pbe_20L/system/controlDict +++ b/experimental_cases/disengagement/bubble_column_pbe_20L/system/controlDict @@ -53,6 +53,12 @@ maxDeltaT 0.0005; functions { + + #includeFunc writeObjects(d.gas) + #includeFunc writeObjects(thermo:rho.gas) + #includeFunc writeObjects(thermo:rho.liquid) + #includeFunc writeObjects(thermo:mu.gas) + #includeFunc writeObjects(thermo:mu.liquid) disengagement { type disengagement; diff --git a/tutorial_cases/FlatPanel_250L_ASU/system/controlDict b/tutorial_cases/FlatPanel_250L_ASU/system/controlDict index baa693c3..40167ce4 100644 --- a/tutorial_cases/FlatPanel_250L_ASU/system/controlDict +++ b/tutorial_cases/FlatPanel_250L_ASU/system/controlDict @@ -50,8 +50,15 @@ maxCo 0.3; maxDeltaT 1; + functions { + + #includeFunc writeObjects(d.gas) + #includeFunc writeObjects(thermo:rho.gas) + #includeFunc writeObjects(thermo:rho.liquid) + #includeFunc writeObjects(thermo:mu.liquid) + #includeFunc writeObjects(thermo:mu.gas) #includeFunc fieldAverage(U.gas, U.liquid, alpha.gas, p) } diff --git a/tutorial_cases/airlift_40m/system/controlDict b/tutorial_cases/airlift_40m/system/controlDict index b9bba7f4..8748e2a7 100644 --- a/tutorial_cases/airlift_40m/system/controlDict +++ b/tutorial_cases/airlift_40m/system/controlDict @@ -51,6 +51,22 @@ maxCo 0.3; maxDeltaT 1; +functions +{ + + #includeFunc writeObjects(d.gas) + #includeFunc writeObjects(thermo:rho.gas) + #includeFunc writeObjects(thermo:rho.liquid) + #includeFunc writeObjects(thermo:mu.liquid) + #includeFunc writeObjects(thermo:mu.gas) + +} +//functions +//{ +// #includeFunc fieldAverage(U.air, U.water, alpha.air, p) +//} + + /*functions { writeFields diff --git a/tutorial_cases/bubble_column_20L/system/controlDict b/tutorial_cases/bubble_column_20L/system/controlDict index ef313841..55951723 100644 --- a/tutorial_cases/bubble_column_20L/system/controlDict +++ b/tutorial_cases/bubble_column_20L/system/controlDict @@ -57,6 +57,8 @@ functions #includeFunc writeObjects(d.gas) #includeFunc writeObjects(thermo:rho.gas) #includeFunc writeObjects(thermo:rho.liquid) + #includeFunc writeObjects(thermo:mu.liquid) + #includeFunc writeObjects(thermo:mu.gas) fieldAverage { diff --git a/tutorial_cases/loop_reactor_mixing/system/controlDict b/tutorial_cases/loop_reactor_mixing/system/controlDict index 204a8cbe..f4665ed8 100644 --- a/tutorial_cases/loop_reactor_mixing/system/controlDict +++ b/tutorial_cases/loop_reactor_mixing/system/controlDict @@ -57,11 +57,10 @@ functions #includeFunc writeObjects(d.gas) #includeFunc writeObjects(thermo:rho.gas) #includeFunc writeObjects(thermo:rho.liquid) + #includeFunc writeObjects(thermo:mu.liquid) + #includeFunc writeObjects(thermo:mu.gas) + #includeFunc fieldAverage(U.air, U.water, alpha.air, p) } -//functions -//{ -// #includeFunc fieldAverage(U.air, U.water, alpha.air, p) -//} // ************************************************************************* // diff --git a/tutorial_cases/loop_reactor_reacting/system/controlDict b/tutorial_cases/loop_reactor_reacting/system/controlDict index 99054d1c..5bf125c6 100644 --- a/tutorial_cases/loop_reactor_reacting/system/controlDict +++ b/tutorial_cases/loop_reactor_reacting/system/controlDict @@ -50,17 +50,17 @@ maxCo 0.5; maxDeltaT 0.01; + functions { - - //#includeFunc writeObjects(d.gas) + + #includeFunc writeObjects(d.gas) #includeFunc writeObjects(thermo:rho.gas) #includeFunc writeObjects(thermo:rho.liquid) + #includeFunc writeObjects(thermo:mu.liquid) + #includeFunc writeObjects(thermo:mu.gas) + #includeFunc fieldAverage(U.air, U.water, alpha.air, p) } -//functions -//{ -// #includeFunc fieldAverage(U.air, U.water, alpha.air, p) -//} // ************************************************************************* // diff --git a/tutorial_cases/stirred_tank/system/controlDict b/tutorial_cases/stirred_tank/system/controlDict index 0f981f30..b207a82e 100644 --- a/tutorial_cases/stirred_tank/system/controlDict +++ b/tutorial_cases/stirred_tank/system/controlDict @@ -51,6 +51,17 @@ maxCo 0.5; maxDeltaT 1; + +functions +{ + + #includeFunc writeObjects(d.gas) + #includeFunc writeObjects(thermo:rho.gas) + #includeFunc writeObjects(thermo:rho.liquid) + #includeFunc writeObjects(thermo:mu.liquid) + #includeFunc writeObjects(thermo:mu.gas) +} + /*functions { writeFields From cf61a60ccf79df3170d37c4f4b9166cd25b66b6f Mon Sep 17 00:00:00 2001 From: Malik Date: Wed, 1 Oct 2025 16:56:29 -0600 Subject: [PATCH 05/17] read bubble diameter from phaseProperties is needed --- bird/postprocess/post_quantities.py | 22 ++++---- bird/utilities/ofio.py | 78 +++++++++++++++++++++++++++++ tests/io/test_read_foam_fields.py | 56 ++++++++++++++++++++- 3 files changed, 142 insertions(+), 14 deletions(-) diff --git a/bird/postprocess/post_quantities.py b/bird/postprocess/post_quantities.py index d9a2243b..515907e5 100644 --- a/bird/postprocess/post_quantities.py +++ b/bird/postprocess/post_quantities.py @@ -884,7 +884,7 @@ def compute_ave_bubble_diam( alpha_liq, field_dict = read_field( field_name="alpha.liquid", field_dict=field_dict, **kwargs ) - d_gas, field_dict = read_field( + d_gas, field_dict = read_bubble_diameter( field_name="d.gas", field_dict=field_dict, **kwargs ) ind_liq, field_dict = _get_ind_liq(field_dict=field_dict, **kwargs) @@ -909,10 +909,9 @@ def compute_ave_bubble_diam( def compute_instantaneous_kla( case_folder: str, time_folder: str, - species_names: list[str], + species_names: str | list[str], n_cells: int | None = None, volume_time: str | None = None, - constant_bubble_diameter: float | None = None, field_dict: dict | None = None, ) -> tuple[dict, dict, dict]: r""" @@ -972,7 +971,7 @@ def compute_instantaneous_kla( Path to case folder time_folder: str Name of time folder to analyze - species_names: list[str] + species_names: str | list[str] List of species name for which to compute kla n_cells : int | None Number of cells in the domain. @@ -980,9 +979,6 @@ def compute_instantaneous_kla( volume_time : str | None Time folder to read to get the cell volumes. If None, finds volume time automatically - constant_bubble_diameter : float | None - Constant bubble diameter value (in m) to handle missing d.gas - If None, reads d.gas field_dict : dict Dictionary of fields used to avoid rereading the same fields to calculate different quantities @@ -1002,6 +998,9 @@ def compute_instantaneous_kla( if field_dict is None: field_dict = {} + if isinstance(species_names, "str"): + species_names = [species_names] + # Read relevant fields kwargs = { "case_folder": case_folder, @@ -1052,12 +1051,9 @@ def compute_instantaneous_kla( U_liq, field_dict = read_field( field_name="U.liquid", field_dict=field_dict, **kwargs ) - try: - d_gas, field_dict = read_field( - field_name="d.gas", field_dict=field_dict, **kwargs - ) - except FileNotFoundError: - d_gas = constant_bubble_diameter + d_gas, field_dict = read_bubble_diameter( + field_name="d.gas", field_dict=field_dict, **kwargs + ) mu_liq, field_dict = read_field( field_name="thermo:mu.liquid", field_dict=field_dict, **kwargs ) diff --git a/bird/utilities/ofio.py b/bird/utilities/ofio.py index 5a72dbc2..fbbcb0f7 100644 --- a/bird/utilities/ofio.py +++ b/bird/utilities/ofio.py @@ -967,6 +967,84 @@ def read_cell_volumes( return cell_volumes, field_dict +def read_bubble_diameter( + case_folder: str, + time_folder: str | None = None, + n_cells: int | None = None, + field_dict: dict | None = None, +) -> tuple[np.ndarray | float, dict]: + """ + Read bubble diameter at a given time and store it in dictionary for later reuse. + A specific function is constructed so that if d.gas is not available, the bubble + diameter is read from phaseProperties. + + Parameters + ---------- + case_folder: str + Path to case folder + time_folder: str | None + Name of time folder to analyze. + If None, it will be found automatically + n_cells : int | None + Number of cells in the domain. + If None, it will deduced from the field reading + field_dict : dict | None + Dictionary of fields used to avoid rereading the same fields to calculate different quantities + + Returns + ---------- + cell_volumes : np.ndarray | float + Field of cell volumes + field_dict : dict + Dictionary of fields read + """ + + if field_dict is None: + field_dict = {} + + kwargs = { + "case_folder": case_folder, + "time_folder": time_folder, + "n_cells": n_cells, + } + + if not ("d.gas" in field_dict) or field_dict["d.gas"] is None: + try: + d_gas, field_dict = read_field( + field_name="d.gas", field_dict=field_dict, **kwargs + ) + + except FileNotFoundError as err: + logger.debug( + "Could not find d.gas, checking if bubble size model is constant" + ) + # d.gas does not exist, it might be because a constant bubble diameter model is used + phaseProperties_dict = read_openfoam_dict( + os.path.join(case_folder, "constant", "phaseProperties") + ) + # If the bubble size model is not constant, raise original error + if ( + not phaseProperties_dict["gas"]["diameterModel"].lower() + == "constant" + ): + logger.error( + f"Bubble size model is not constant ({phaseProperties_dict['gas']['diameterModel']}), yet could not find d.gas" + ) + raise FileNotFoundError(err) + else: + # Get the constant bubble diameter + logger.debug("Reading bubble diameter from phaseProperties") + d_gas = float( + phaseProperties_dict["gas"]["constantCoeffs"]["d"] + ) + + else: + # Get field from dict if it has been read before + d_gas = field_dict["d.gas"] + + return d_gas, field_dict + + def _cross_reference_global_vars( globalVars_dict: dict, ) -> dict: diff --git a/tests/io/test_read_foam_fields.py b/tests/io/test_read_foam_fields.py index a991b8d2..c0beed96 100644 --- a/tests/io/test_read_foam_fields.py +++ b/tests/io/test_read_foam_fields.py @@ -1,9 +1,15 @@ import os +import shutil from pathlib import Path import numpy as np -from bird.utilities.ofio import _readOF, _readOFScal, _readOFVec +from bird.utilities.ofio import ( + _readOF, + _readOFScal, + _readOFVec, + read_bubble_diameter, +) def test_read_nonunif_scal(): @@ -132,3 +138,51 @@ def test_read_unif_vec(): ) assert np.linalg.norm(data_dict["field"] - [0.0, 0.1, 0.0]) < 1e-6 assert data_dict["name"] == "U_unif_dummy" + + +def test_read_bubble_diameter(): + """ + Test for reading bubble diameter + """ + case_folder = os.path.join( + Path(__file__).parent, + "..", + "..", + "bird", + "postprocess", + "data_conditional_mean", + ) + # Read bubble diameter + d_gas, _ = read_bubble_diameter(case_folder=case_folder, time_folder="80") + assert abs(d_gas[2] - 0.00966694) < 1e-8 + # Read bubble diameter if no d.gas is available + # Manufacture a case without d.gas + shutil.move( + os.path.join(case_folder, "80", "d.gas"), + os.path.join(case_folder, "d.gas_tmp"), + ) + shutil.move( + os.path.join(case_folder, "constant", "phaseProperties"), + os.path.join(case_folder, "phaseProperties_tmp"), + ) + shutil.copy( + os.path.join(case_folder, "constant", "phaseProperties_constantD"), + os.path.join(case_folder, "constant", "phaseProperties"), + ) + d_gas, _ = read_bubble_diameter(case_folder=case_folder, time_folder="80") + assert abs(d_gas - 0.003) < 1e-8 + shutil.move( + os.path.join(case_folder, "d.gas_tmp"), + os.path.join(case_folder, "80", "d.gas"), + ) + shutil.move( + os.path.join(case_folder, "phaseProperties_tmp"), + os.path.join(case_folder, "constant", "phaseProperties"), + ) + # Read bubble diameter + d_gas, _ = read_bubble_diameter(case_folder=case_folder, time_folder="80") + assert abs(d_gas[2] - 0.00966694) < 1e-8 + + +if __name__ == "__main__": + test_read_bubble_diameter() From db47a171f1426cc9383a41ddc4d5ef97b43cb16b Mon Sep 17 00:00:00 2001 From: Malik Date: Wed, 1 Oct 2025 17:17:33 -0600 Subject: [PATCH 06/17] test passing species as string --- tests/postprocess/test_post_quantities.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/tests/postprocess/test_post_quantities.py b/tests/postprocess/test_post_quantities.py index 904a5bc3..812abab1 100644 --- a/tests/postprocess/test_post_quantities.py +++ b/tests/postprocess/test_post_quantities.py @@ -307,6 +307,17 @@ def test_instantaneous_kla(): abs(cstar_spec4["CO2"] - cstar_spec1["CO2"]) / cstar_spec1["CO2"] > 1e-6 ) + kla_spec5, cstar_spec5, _ = compute_instantaneous_kla( + species_names="CO2", + case_folder=case_folder, + time_folder="80", + ) + # Make sure values change over time + assert abs(kla_spec5["CO2"] - kla_spec1["CO2"]) / kla_spec1["CO2"] > 1e-6 + assert ( + abs(cstar_spec5["CO2"] - cstar_spec1["CO2"]) / cstar_spec1["CO2"] + > 1e-6 + ) if __name__ == "__main__": From 368a9c95d3aa3bd1b3a4999be2bf19935f28d784 Mon Sep 17 00:00:00 2001 From: Malik Date: Wed, 1 Oct 2025 17:33:15 -0600 Subject: [PATCH 07/17] add phase properties for test --- .../constant/phaseProperties | 303 ++++++++++++++++++ .../constant/phaseProperties_constantD | 237 ++++++++++++++ bird/postprocess/post_quantities.py | 10 +- bird/version.py | 2 +- 4 files changed, 544 insertions(+), 8 deletions(-) create mode 100644 bird/postprocess/data_conditional_mean/constant/phaseProperties create mode 100644 bird/postprocess/data_conditional_mean/constant/phaseProperties_constantD diff --git a/bird/postprocess/data_conditional_mean/constant/phaseProperties b/bird/postprocess/data_conditional_mean/constant/phaseProperties new file mode 100644 index 00000000..586de4fd --- /dev/null +++ b/bird/postprocess/data_conditional_mean/constant/phaseProperties @@ -0,0 +1,303 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 9 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + format ascii; + class dictionary; + object phaseProperties; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#include "$FOAM_CASE/constant/globalVars" + +//type populationBalanceMultiphaseSystem; + +type interfaceCompositionPhaseChangePopulationBalanceMultiphaseSystem; + +phases (gas liquid); + +populationBalances (bubbles); + +gas +{ + type multiComponentPhaseModel;//pureIsothermalPhaseModel; + + diameterModel velocityGroup; + + velocityGroupCoeffs + { + populationBalance bubbles; + + shapeModel spherical; + + sizeGroups + ( + f1 {dSph 1e-3; value 0.0;} + f2 {dSph 1.5e-3; value 0.0;} + f3 {dSph 2e-3; value 1.0;} + f4 {dSph 2.5e-3; value 0.0;} + f5 {dSph 3e-3; value 0.0;} + f6 {dSph 3.5e-3; value 0.0;} + f7 {dSph 4e-3; value 0.0;} + f8 {dSph 4.5e-3; value 0.0;} + f9 {dSph 5e-3; value 0.0;} + f10 {dSph 5.5e-3; value 0.0;} + f11 {dSph 6e-3; value 0.0;} + f12 {dSph 6.5e-3; value 0.0;} + f13 {dSph 7e-3; value 0.0;} + f14 {dSph 7.5e-3; value 0.0;} + f15 {dSph 8e-3; value 0.0;} + f16 {dSph 8.5e-3; value 0.0;} + f17 {dSph 9e-3; value 0.0;} + f18 {dSph 9.5e-3; value 0.0;} + f19 {dSph 10e-3; value 0.0;} + f20 {dSph 10.5e-3; value 0.0;} + f21 {dSph 11e-3; value 0.0;} + ); + } + + residualAlpha 1e-6; +} + +liquid +{ + type multiComponentPhaseModel;//pureIsothermalPhaseModel; + + diameterModel constant; + + constantCoeffs + { + d 1e-4; + } + Sc #codeStream + { + code + #{ + os << ($LeLiqMix * $CpMixLiq * $muMixLiq / $kThermLiq); + #}; + }; + + residualAlpha 1e-6; +} + +populationBalanceCoeffs +{ + bubbles + { + continuousPhase liquid; + + coalescenceModels + ( + LehrMilliesMewes{ + uCrit 0.08; + alphaMax 0.6; + } + ); + + binaryBreakupModels + ( + // LehrMilliesMewes{} + ); + + breakupModels + ( + Laakkonen{ + daughterSizeDistributionModel Laakkonen; + } + ); + + driftModels + ( + densityChange{} + ); + + nucleationModels + (); + } +} + +blending +{ + default + { + type linear; + minFullyContinuousAlpha.gas 0.7; + minPartlyContinuousAlpha.gas 0.3; + minFullyContinuousAlpha.liquid 0.7; + minPartlyContinuousAlpha.liquid 0.3; + } + heatTransfer + { + type linear; + minFullyContinuousAlpha.gas 1; + minPartlyContinuousAlpha.gas 0; + minFullyContinuousAlpha.liquid 1; + minPartlyContinuousAlpha.liquid 0; + } + massTransfer + { + $heatTransfer; + } +} + +surfaceTension +( + (gas and liquid) + { + type constant; + sigma $sigmaLiq; + } +); + +interfaceCompression +(); + +aspectRatio +( + (gas in liquid) + { + type Wellek; + } +); + +drag +( + (gas in liquid) + { + type IshiiZuber; + + residualRe 1e-3; + swarmCorrection + { + type none; + } + } +); + +virtualMass +( + (gas in liquid) + { + type constantCoefficient; + Cvm 0.5; + } +); + +// heatTransfer +// (); + +heatTransfer.gas +( + (gas in liquid) + { + type spherical; + residualAlpha 1e-4; + } + + (liquid in gas) + { + type RanzMarshall; + residualAlpha 1e-4; + } +); + +heatTransfer.liquid +( + (gas in liquid) + { + type RanzMarshall; + residualAlpha 1e-4; + } + + (liquid in gas) + { + type spherical; + residualAlpha 1e-4; + } +); + +interfaceComposition.gas +(); + +interfaceComposition.liquid +( + (liquid and gas) + { + type Henry; + species ( CO2 CO H2); + k ( 6.4e-1 2.0e-2 1.8e-2); + Le $LeLiqMix; + } +); + +diffusiveMassTransfer.gas +(); + +diffusiveMassTransfer.liquid +( + (gas in liquid) + { + type Higbie; + Le $LeLiqMix; + } + + (liquid in gas) + { + type spherical; + Le 1.0; //not used for spherical + } +); + +phaseTransfer +(); + +lift +( + (gas in liquid) + { + type wallDamped; + + wallDamping + { + type cosine; + Cd 3.0; + } + + lift + { + type Tomiyama; + + swarmCorrection + { + type none; + } + } + } + +); + +wallLubrication +( + (gas in liquid) + { + type Antal; + Cw1 -0.01; + Cw2 0.05; + } +); + +turbulentDispersion +( + (gas in liquid) + { + type Burns; + sigma 0.9; + } +); + +// ************************************************************************* // diff --git a/bird/postprocess/data_conditional_mean/constant/phaseProperties_constantD b/bird/postprocess/data_conditional_mean/constant/phaseProperties_constantD new file mode 100644 index 00000000..78650285 --- /dev/null +++ b/bird/postprocess/data_conditional_mean/constant/phaseProperties_constantD @@ -0,0 +1,237 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 9 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + format ascii; + class dictionary; + object phaseProperties; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#include "$FOAM_CASE/constant/globalVars" + +type interfaceCompositionPhaseChangeMultiphaseSystem; + +phases (gas liquid); + +gas +{ + type multiComponentPhaseModel;//pureIsothermalPhaseModel; + + diameterModel constant; + + constantCoeffs + { + d 3e-3; + } + residualAlpha 1e-6; + Sc 0.7; +} + + +liquid +{ + type multiComponentPhaseModel;//pureIsothermalPhaseModel; + + diameterModel constant; + + constantCoeffs + { + d 1e-4; + } + Sc #codeStream + { + code + #{ + os << ($LeLiqMix * $CpMixLiq * $muMixLiq / $kThermLiq); + #}; + }; + + residualAlpha 1e-6; +} + +blending +{ + default + { + type linear; + minFullyContinuousAlpha.gas 0.7; + minPartlyContinuousAlpha.gas 0.3; + minFullyContinuousAlpha.liquid 0.7; + minPartlyContinuousAlpha.liquid 0.3; + } + heatTransfer + { + type linear; + minFullyContinuousAlpha.gas 1; + minPartlyContinuousAlpha.gas 0; + minFullyContinuousAlpha.liquid 1; + minPartlyContinuousAlpha.liquid 0; + } + massTransfer + { + $heatTransfer; + } +} + +surfaceTension +( + (gas and liquid) + { + type constant; + sigma $sigmaLiq; + } +); + +interfaceCompression +(); + +aspectRatio +( + (gas in liquid) + { + type Wellek; + } +); + +drag +( + (gas in liquid) + { + type IshiiZuber; + + residualRe 1e-3; + swarmCorrection + { + type none; + } + } +); + +virtualMass +( + (gas in liquid) + { + type constantCoefficient; + Cvm 0.5; + } +); + +// heatTransfer +// (); + +heatTransfer.gas +( + (gas in liquid) + { + type spherical; + residualAlpha 1e-4; + } + + (liquid in gas) + { + type RanzMarshall; + residualAlpha 1e-4; + } +); + +heatTransfer.liquid +( + (gas in liquid) + { + type RanzMarshall; + residualAlpha 1e-4; + } + + (liquid in gas) + { + type spherical; + residualAlpha 1e-4; + } +); + +interfaceComposition.gas +(); + +interfaceComposition.liquid +( + (liquid and gas) + { + type Henry; + species ( CO2 CO H2); + k ( 6.4e-1 2.0e-2 1.8e-2); + Le $LeLiqMix; + } +); + +diffusiveMassTransfer.gas +(); + +diffusiveMassTransfer.liquid +( + (gas in liquid) + { + type Higbie; + Le $LeLiqMix; + } + + (liquid in gas) + { + type spherical; + Le 1.0; //not used for spherical + } +); + +phaseTransfer +(); + +lift +( + (gas in liquid) + { + type wallDamped; + + wallDamping + { + type cosine; + Cd 3.0; + } + + lift + { + type Tomiyama; + + swarmCorrection + { + type none; + } + } + } + +); + +wallLubrication +( + (gas in liquid) + { + type Antal; + Cw1 -0.01; + Cw2 0.05; + } +); + +turbulentDispersion +( + (gas in liquid) + { + type Burns; + sigma 0.9; + } +); + +// ************************************************************************* // diff --git a/bird/postprocess/post_quantities.py b/bird/postprocess/post_quantities.py index 515907e5..279403a9 100644 --- a/bird/postprocess/post_quantities.py +++ b/bird/postprocess/post_quantities.py @@ -884,9 +884,7 @@ def compute_ave_bubble_diam( alpha_liq, field_dict = read_field( field_name="alpha.liquid", field_dict=field_dict, **kwargs ) - d_gas, field_dict = read_bubble_diameter( - field_name="d.gas", field_dict=field_dict, **kwargs - ) + d_gas, field_dict = read_bubble_diameter(field_dict=field_dict, **kwargs) ind_liq, field_dict = _get_ind_liq(field_dict=field_dict, **kwargs) cell_volume, field_dict = read_cell_volumes( @@ -998,7 +996,7 @@ def compute_instantaneous_kla( if field_dict is None: field_dict = {} - if isinstance(species_names, "str"): + if isinstance(species_names, str): species_names = [species_names] # Read relevant fields @@ -1051,9 +1049,7 @@ def compute_instantaneous_kla( U_liq, field_dict = read_field( field_name="U.liquid", field_dict=field_dict, **kwargs ) - d_gas, field_dict = read_bubble_diameter( - field_name="d.gas", field_dict=field_dict, **kwargs - ) + d_gas, field_dict = read_bubble_diameter(field_dict=field_dict, **kwargs) mu_liq, field_dict = read_field( field_name="thermo:mu.liquid", field_dict=field_dict, **kwargs ) diff --git a/bird/version.py b/bird/version.py index 7b1b5e3f..0461d1e7 100644 --- a/bird/version.py +++ b/bird/version.py @@ -1,3 +1,3 @@ """Bio reactor design version""" -__version__ = "0.0.40" +__version__ = "0.0.42" From dc19482ca2506a0146d52d966bd364fb24d12a05 Mon Sep 17 00:00:00 2001 From: Malik Date: Thu, 2 Oct 2025 06:16:00 -0600 Subject: [PATCH 08/17] add kla fitting function --- applications/compute_kla_uq.py | 11 +- bird/postprocess/kla_utils.py | 56 +++--- bird/postprocess/post_quantities.py | 162 ++++++++++++++++-- docs/source/python_interface_tut.rst | 1 - .../postprocess/python_interface_tut.py | 1 - 5 files changed, 176 insertions(+), 55 deletions(-) diff --git a/applications/compute_kla_uq.py b/applications/compute_kla_uq.py index 6afec4f7..3eb46291 100644 --- a/applications/compute_kla_uq.py +++ b/applications/compute_kla_uq.py @@ -1,6 +1,8 @@ import argparse import os +import numpy as np + from bird import BIRD_KLA_DATA_DIR from bird.postprocess.kla_utils import compute_kla, print_res_dict @@ -55,10 +57,13 @@ else: bootstrap = True +data = np.loadtxt(args.data_file) +data_t = data[:, args.time_index] +data_c = data[:, args.conc_index] + res_dict = compute_kla( - filename=args.data_file, - time_ind=args.time_index, - conc_ind=args.conc_index, + data_t, + data_c, bootstrap=bootstrap, max_chop=args.max_chop, ) diff --git a/bird/postprocess/kla_utils.py b/bird/postprocess/kla_utils.py index 28095d0d..31dcece6 100644 --- a/bird/postprocess/kla_utils.py +++ b/bird/postprocess/kla_utils.py @@ -1,7 +1,9 @@ -import logging - import corner import jax +from jax._src import config + +config.update("jax_platforms", "cpu") + import jax.numpy as jnp import jax.random as random import numpy as np @@ -19,24 +21,14 @@ from numpyro.infer import MCMC, NUTS from prettyPlot.plotting import * -logger = logging.getLogger(__name__) - - -def read_data( - filename: str, time_ind: int, conc_ind: int, ind_start=0 -) -> tuple[np.ndarray, np.ndarray]: - """ - Read observation data - """ - tmp = np.loadtxt(filename) - return tmp[ind_start:, time_ind], tmp[ind_start:, conc_ind] +from bird import logger def c_model( theta: list, t: np.ndarray | jnp.ndarray, t0: float, c0: float ) -> np.ndarray | jnp.ndarray: """ - Concentration vs time function involving kLa + Model of concentration vs time function involving kLa """ cstar, kla = theta return (cstar - c0) * (1.0 - jnp.exp(-kla * (t - t0))) + c0 @@ -47,7 +39,7 @@ def bayes_step( concentration_obs: np.ndarray | jnp.ndarray, ): """ - Rule for MCMC sampling + Prior and likelihood definition for MCMC sampling """ # define parameters (incl. prior ranges) cstar = numpyro.sample("cstar", dist.Uniform(1e-8, 1000)) @@ -61,7 +53,7 @@ def bayes_step( # concentration predictions c_pred = c_model([cstar, kla], t=time_obs, t0=t0, c0=c0) - # likelihood uncertainty + # likelihood uncertainty (assumed homoskedastic) std_c = jnp.ones(c_pred.shape[0]) * sigma # MCMC sampling with multivariate normal @@ -70,7 +62,7 @@ def bayes_step( def hmc_samples_2_np_samples(hmc_samples: dict) -> np.ndarray: """ - go from hmc samples generated by numpyro to numpy + Convert hmc samples generated by numpyro to numpy """ labels = list(hmc_samples.keys()) nsamples = len(hmc_samples[labels[0]]) @@ -96,7 +88,8 @@ def np_samples_2_pred( concentration_obs: np.ndarray, ) -> dict: """ - go from parameter samples to predictions + Go from parameter samples to predictions. + Useful for accuracy calculation and plotting. """ # Uncertainty propagation nsamples = np_hmc_samples.shape[0] @@ -130,7 +123,7 @@ def compute_accuracy( hmc_samples: dict, time_obs: np.ndarray, concentration_obs: np.ndarray ) -> float: """ - Compute how well the sampled parameter explain the observed data + Compute how well the sampled parameter explain the observed concentration data """ np_hmc_samples = hmc_samples_2_np_samples(hmc_samples) pred_dict = np_samples_2_pred(np_hmc_samples, time_obs, concentration_obs) @@ -152,11 +145,10 @@ def check_data_size(time_obs: np.ndarray, concentration_obs: np.ndarray): def compute_kla( - filename: str, - time_ind: int, - conc_ind: int, + data_t: np.ndarray, + data_c: np.ndarray, max_chop: int | None = None, - bootstrap: bool = False, + bootstrap: bool = True, ) -> dict: """ Main loop to compute kla @@ -172,15 +164,18 @@ def compute_kla( kla_err = [] cstar_err = [] ind = [] - data_t_tmp, data_c_tmp = read_data(filename, time_ind, conc_ind, 0) - check_data_size(data_t_tmp, data_c_tmp) + check_data_size(data_t, data_c) + + data_t_tmp = data_t.copy() + data_c_tmp = data_c.copy() # Find where to start in the timeseries for ind_start in range(len(data_t_tmp) - 5): if max_chop is not None: max_chop -= 1 logger.debug(f"Chopping index = {ind_start}") - data_t, data_c = read_data(filename, time_ind, conc_ind, ind_start) + data_t = data_t_tmp[ind_start:] + data_c = data_c_tmp[ind_start:] # Hamiltonian Monte Carlo (HMC) with no u turn sampling (NUTS) kernel = NUTS(bayes_step, target_accept_prob=0.9) @@ -287,9 +282,6 @@ def compute_kla( cstar_err_boot = [cstar_err[-1]] return { - "filename": filename, - "time_ind": time_ind, - "conc_ind": conc_ind, "kla": np.mean(np.array(kla_boot)), "kla_err": np.mean(np.array(kla_err_boot)), "cstar": np.mean(np.array(cstar_boot)), @@ -360,9 +352,6 @@ def print_res_dict(res_dict: dict) -> None: """ Log kla and cstar to screen """ - file = res_dict["filename"] - t_ind = res_dict["time_ind"] - c_ind = res_dict["conc_ind"] kla = res_dict["kla"] kla_err = res_dict["kla_err"] cs = res_dict["cstar"] @@ -375,9 +364,6 @@ def print_res_dict(res_dict: dict) -> None: bs = res_dict["bootstrapped"] - logger.info( - f"For {file} with time index: {t_ind}, concentration index: {c_ind}" - ) if bs: logger.info(f"\tkla = {kla:.4g} +/- {kla_err:.4g}") logger.info(f"\tcstar = {cs:.4g} +/- {cs_err:.4g}") diff --git a/bird/postprocess/post_quantities.py b/bird/postprocess/post_quantities.py index 279403a9..f89d33d9 100644 --- a/bird/postprocess/post_quantities.py +++ b/bird/postprocess/post_quantities.py @@ -5,6 +5,8 @@ from bird import logger from bird.utilities.ofio import * +from .kla_utils import compute_kla + def _field_filter( field: float | np.ndarray, ind: np.ndarray, field_type: str @@ -721,7 +723,6 @@ def compute_ave_conc_liq( volume_time: str | None = None, spec_name: str = "CO2", mol_weight: float = 0.04401, - rho_val: float | None = 1000, field_dict: dict | None = None, ) -> tuple[float, dict]: r""" @@ -753,8 +754,6 @@ def compute_ave_conc_liq( Name of the species mol_weight : float Molecular weight of species (kg/mol) - rho_val : float | None - Constant density not available from time folder (kg/m3) field_dict : dict Dictionary of fields used to avoid rereading the same fields to calculate different quantities @@ -797,18 +796,17 @@ def compute_ave_conc_liq( cell_volume, field_dict = read_cell_volumes( field_dict=field_dict, **kwargs_vol ) - - # Density of liquid is not always printed (special case) - if not ("rho_liq" in field_dict) or field_dict["rho_liq"] is None: - if rho_val is not None: - rho_liq = rho_val - field_dict["rho_liq"] = rho_val - else: - rho_liq_file = os.path.join(case_folder, time_folder, "rhom") - rho_liq = _readOFScal(rho_liq_file, n_cells)["field"] - field_dict["rho_liq"] = rho_liq - else: - rho_liq = field_dict["rho_liq"] + try: + rho_liq, field_dict = read_field( + field_name="thermo:rho.liquid", field_dict=field_dict, **kwargs + ) + except FileNotFoundError: + abs_time_path = os.path.join(case_folder, time_folder) + logger.warning( + f"thermo:rho.liquid not found in {abs_time_path}, assuming it is 1000kg/m3" + ) + rho_liq = 1000 + field_dict["rho_liq"] = rho_liq # Only compute over the liquid alpha_liq = _field_filter(alpha_liq, ind=ind_liq, field_type="scalar") @@ -1117,3 +1115,137 @@ def compute_instantaneous_kla( ) / np.sum(cell_volume * alpha_liq) return kla_spec, cstar_spec, field_dict + + +def compute_fitted_kla( + case_folder: str, + species_names: str | list[str], + n_cells: int | None = None, + volume_time: str | None = None, + field_dict: dict | None = None, +) -> tuple[dict, dict, dict]: + r""" + Calculate :math:`kLa_{\rm spec}` and saturation concentration (:math:`C^*_{\rm spec}`) for a list of species from time series data (rather than instantaneously). + + Given a time series of concentration of species, the following expression is fitted + + .. math:: + [spec](t) = [spec]^* (1 - \operatorname{exp}(-{kLa}_{\rm spec} t)). + + where + + - :math:`kLa_{\rm spec}` is the mass transfer rate of species :math:`\rm spec` in :math:`h^{-1}` + - :math:`t` is the time in :math:`s` + - :math:`[spec]^*` is the estimated saturation concentration of species :math:`\rm spec` in :math:`mol/m^3` + - :math:`[spec](t)` is the instantaneous liquid volume averaged concentration of species :math:`\rm spec` in :math:`mol/m^3` + + Both :math:`[spec]^*` and :math:`kLa_{\rm spec}` are fitted. + The fit is done with Markov Chain Monte Carlo which outputs samples of the posterior PDF of :math:`[spec]^*` and :math:`kLa_{\rm spec}`. + + Parameters + ---------- + case_folder: str + Path to case folder + species_names: str | list[str] + List of species name for which to compute kla + n_cells : int | None + Number of cells in the domain. + If None, it will deduced from the field reading + volume_time : str | None + Time folder to read to get the cell volumes. + If None, finds volume time automatically + field_dict : dict + Dictionary of fields used to avoid rereading the same fields to calculate different quantities + + Returns + ---------- + kla_spec: dict + Instantaneous volume averaged kLa for each species + Keys are species names + Values are the kLa values + cstar_spec: dict + Instantaneous volume averaged cstar for each species + Keys are species names + Values are the cstar values + field_dict : dict + Dictionary of fields read + """ + if field_dict is None: + field_dict = {} + + if isinstance(species_names, "str"): + species_names = [species_names] + + # Read relevant fields + kwargs = { + "case_folder": case_folder, + "n_cells": n_cells, + "volume_time": volume_time, + } + + # Get all the time folders + time_float_sorted, time_str_sorted = get_case_times(case_folder) + + # Read globarVars into a python dict + # Replace all the #calc entries with their numeral values + globalVars = read_global_vars(case_folder=case_folder, cross_ref=True) + + # Check that global vars has the values we want and provide a useful error message otherwise + for species_name in species_names: + if not f"Mw_{species_name}" in globalVars: + err_msg = f"Mw_{species_name} was not found in globalVars." + err_msg += "\nIf you add it, it should be [kg/mol]" + raise KeyError(err_msg) + + # Initialize the mesh fields + mesh_field_dict = {} + if "cell_centers" in field_dict: + mesh_field_dict["cell_centers"] = field_dict["cell_centers"] + else: + mesh_field_dict["cell_centers"], _ = read_cell_centers(case_folder) + if "V" in field_dict: + mesh_field_dict["V"] = field_dict["V"] + else: + mesh_field_dict["V"], _ = read_cell_volumes(case_folder) + + # Initialize the data structure for concentration + c_history = {} + for species_name in species_names: + c_history[species_name] = np.zeros(len(time_str_sorted)) + + # Read concentration + for itime, time_folder in enumerate(time_str_sorted): + + # Reinitialize kla field dict + kla_field_dict = {} + for key in mesh_field_dict: + kla_field_dict[key] = mesh_field_dict[key] + + # Compute reactor averaged liquid concentration for all the species + for species_name in species_names: + c_liq, kla_field_dict = compute_ave_conc_liq( + time_folder=time_folder, + spec_name=species_name, + mol_weight=globalVars[f"Mw_{species_name}"], + field_dict=kla_field_dict, + **kwargs, + ) + c_history[species_name][itime] = c_liq + + # Compute kla + kla_spec = {} + cstar_spec = {} + for species_name in species_names: + kla_res = compute_kla( + np.array(time_float_sorted), c_history[species_name] + ) + kla_spec[species_name] = { + "mean": kla_res["kla"], + "std": kla_res["kla_err"], + } + cstar_spec[species_name] = { + "mean": kla_res["cstar"], + "std": kla_res["cstar_err"], + } + + return kla_spec, cstar_spec, field_dict diff --git a/docs/source/python_interface_tut.rst b/docs/source/python_interface_tut.rst index 69e05dbd..5cf6d15b 100644 --- a/docs/source/python_interface_tut.rst +++ b/docs/source/python_interface_tut.rst @@ -138,7 +138,6 @@ Several of these quantities, will require reading and processing the same fields c_ave_co2, field_dict = compute_ave_conc_liq( spec_name="CO2", mol_weight=0.04401, - rho_val=1000, field_dict=field_dict, **kwargs, ) diff --git a/tutorial_cases/postprocess/python_interface_tut.py b/tutorial_cases/postprocess/python_interface_tut.py index 085713db..98bea46a 100644 --- a/tutorial_cases/postprocess/python_interface_tut.py +++ b/tutorial_cases/postprocess/python_interface_tut.py @@ -69,7 +69,6 @@ c_ave_co2, field_dict = compute_ave_conc_liq( spec_name="CO2", mol_weight=0.04401, - rho_val=1000, field_dict=field_dict, **kwargs, ) From b7ce945ba676433a76f52005a34d866f6bf97ddd Mon Sep 17 00:00:00 2001 From: Malik Date: Thu, 2 Oct 2025 06:25:44 -0600 Subject: [PATCH 09/17] fix the kla test --- tests/postprocess/test_kla.py | 23 +++++++++++++++++------ 1 file changed, 17 insertions(+), 6 deletions(-) diff --git a/tests/postprocess/test_kla.py b/tests/postprocess/test_kla.py index 66ee5e80..2cb192d9 100644 --- a/tests/postprocess/test_kla.py +++ b/tests/postprocess/test_kla.py @@ -1,6 +1,7 @@ import os from pathlib import Path +import numpy as np import pytest from bird.postprocess.kla_utils import compute_kla, print_res_dict @@ -13,15 +14,25 @@ (False, 2), ], ) -def test_kla(bootstrap, max_chop): - BIRD_KLA_DATA_DIR = os.path.join( - Path(__file__).parent, "..", "..", "bird", "postprocess", "data_kla" +def test_fitted_kla(bootstrap, max_chop): + BIRD_KLA_DATA_FILE = os.path.join( + Path(__file__).parent, + "..", + "..", + "bird", + "postprocess", + "data_kla", + "volume_avg.dat", ) + data = np.loadtxt(BIRD_KLA_DATA_FILE) + data_t = data[:, 0] + data_c = data[:, 1] + res_dict = compute_kla( - os.path.join(BIRD_KLA_DATA_DIR, "volume_avg.dat"), - time_ind=0, - conc_ind=1, + data_t, + data_c, bootstrap=bootstrap, max_chop=max_chop, ) + print_res_dict(res_dict) From 2fb5fc69dedadebb79261978fc58b1ddc87e4bca Mon Sep 17 00:00:00 2001 From: Malik Date: Thu, 2 Oct 2025 08:38:33 -0600 Subject: [PATCH 10/17] read molecular weight from case if possible --- bird/postprocess/computeQoI/compute_QOI.py | 15 ++-- bird/postprocess/post_quantities.py | 46 ++++++------ .../read_history.py | 10 +-- bird/utilities/ofio.py | 75 +++++++++++++++++++ docs/source/python_interface_tut.rst | 5 +- tests/postprocess/test_post_quantities.py | 25 +++---- .../loop_reactor_mixing/read_history.py | 10 +-- .../postprocess/python_interface_tut.py | 5 +- 8 files changed, 124 insertions(+), 67 deletions(-) diff --git a/bird/postprocess/computeQoI/compute_QOI.py b/bird/postprocess/computeQoI/compute_QOI.py index 85012408..c8502041 100644 --- a/bird/postprocess/computeQoI/compute_QOI.py +++ b/bird/postprocess/computeQoI/compute_QOI.py @@ -115,7 +115,7 @@ def get_var(case_folder, time_folder, cell_centers, nCells, val_dict, name): case_folder, time_folder, nCells, - spec_name="CO2", + species_name="CO2", field_dict=val_dict, ) elif name.lower() == "co_liq": @@ -123,7 +123,7 @@ def get_var(case_folder, time_folder, cell_centers, nCells, val_dict, name): case_folder, time_folder, nCells, - spec_name="CO", + species_name="CO", field_dict=val_dict, ) elif name.lower() == "h2_liq": @@ -131,7 +131,7 @@ def get_var(case_folder, time_folder, cell_centers, nCells, val_dict, name): case_folder, time_folder, nCells, - spec_name="H2", + species_name="H2", field_dict=val_dict, ) elif name.lower() == "c_co2_liq": @@ -139,8 +139,7 @@ def get_var(case_folder, time_folder, cell_centers, nCells, val_dict, name): case_folder, time_folder, nCells, - spec_name="CO2", - mol_weight=44.00995 * 1e-3, + species_name="CO2", field_dict=val_dict, ) elif name.lower() == "c_co_liq": @@ -148,8 +147,7 @@ def get_var(case_folder, time_folder, cell_centers, nCells, val_dict, name): case_folder, time_folder, nCells, - spec_name="CO", - mol_weight=28.01055 * 1e-3, + species_name="CO", field_dict=val_dict, ) elif name.lower() == "c_h2_liq": @@ -157,8 +155,7 @@ def get_var(case_folder, time_folder, cell_centers, nCells, val_dict, name): case_folder, time_folder, nCells, - spec_name="H2", - mol_weight=2.01594 * 1e-3, + species_name="H2", field_dict=val_dict, ) else: diff --git a/bird/postprocess/post_quantities.py b/bird/postprocess/post_quantities.py index f89d33d9..f7897526 100644 --- a/bird/postprocess/post_quantities.py +++ b/bird/postprocess/post_quantities.py @@ -636,9 +636,9 @@ def compute_superficial_gas_velocity( def compute_ave_y_liq( case_folder: str, time_folder: str, + species_name: str = "CO2", n_cells: int | None = None, volume_time: str | None = None, - spec_name: str = "CO2", field_dict: dict | None = None, ) -> tuple[float, dict]: r""" @@ -665,7 +665,7 @@ def compute_ave_y_liq( volume_time : str | None Time folder to read to get the cell volumes. If None, finds volume time automatically - spec_name : str + species_name : str Name of the species field_dict : dict | None Dictionary of fields used to avoid rereading the same fields to calculate different quantities @@ -719,10 +719,9 @@ def compute_ave_y_liq( def compute_ave_conc_liq( case_folder: str, time_folder: str, + species_name: str = "CO2", n_cells: int | None = None, volume_time: str | None = None, - spec_name: str = "CO2", - mol_weight: float = 0.04401, field_dict: dict | None = None, ) -> tuple[float, dict]: r""" @@ -744,16 +743,14 @@ def compute_ave_conc_liq( Path to case folder time_folder: str Name of time folder to analyze + species_name : str + Name of the species n_cells : int | None Number of cells in the domain. If None, it will deduced from the field reading volume_time : str | None Time folder to read to get the cell volumes. If None, finds volume time automatically - spec_name : str - Name of the species - mol_weight : float - Molecular weight of species (kg/mol) field_dict : dict Dictionary of fields used to avoid rereading the same fields to calculate different quantities @@ -767,12 +764,12 @@ def compute_ave_conc_liq( if field_dict is None: field_dict = {} + mol_weight = _species_name_to_mw( + case_folder=case_folder, species_name=species_name + ) logger.debug( - f"Computing concentration for {spec_name} with molecular weight {mol_weight:.4g} kg/mol" + f"Computing concentration for {species_name} with molecular weight {mol_weight:.4g} kg/mol" ) - if rho_val is not None: - rho_val = float(rho_val) - logger.debug(f"Assuming liquid density {rho_val} kg/m3") # Read relevant fields kwargs = { @@ -789,7 +786,7 @@ def compute_ave_conc_liq( field_name="alpha.liquid", field_dict=field_dict, **kwargs ) y_liq, field_dict = read_field( - field_name=f"{spec_name}.liquid", field_dict=field_dict, **kwargs + field_name=f"{species_name}.liquid", field_dict=field_dict, **kwargs ) ind_liq, field_dict = _get_ind_liq(field_dict=field_dict, **kwargs) @@ -1014,19 +1011,19 @@ def compute_instantaneous_kla( globalVars = read_global_vars(case_folder=case_folder, cross_ref=True) # Check that global vars has the values we want and provide a useful error message otherwise + mw_species = {} for species_name in species_names: if not f"He_{species_name}" in globalVars: err_msg = f"He_{species_name} was not found in globalVars." err_msg += f'\nIf you add it, it should be looking like #calc "$H_{species_name}_298 * exp($DH_{species_name} *(1. / $T0 - 1./298.15))";' raise KeyError(err_msg) - if not f"Mw_{species_name}" in globalVars: - err_msg = f"Mw_{species_name} was not found in globalVars." - err_msg += "\nIf you add it, it should be [kg/mol]" - raise KeyError(err_msg) if not f"D_{species_name}" in globalVars: err_msg = f"D_{species_name} was not found in globalVars." err_msg += f'\nIf you add it, it should be looking like #calc "1.173e-16 * pow($WC_psi * $WC_M,0.5) * $T0 / $muMixLiq / pow($WC_V_{species_name},0.6)";' raise KeyError(err_msg) + mw_species[species_name] = _species_name_to_mw( + case_folder=case_folder, species_name=species_name + ) # Get liquid domain ind_liq, field_dict = _get_ind_liq(field_dict=field_dict, **kwargs) @@ -1096,7 +1093,7 @@ def compute_instantaneous_kla( rho_gas * species_gas[species_name] * globalVars[f"He_{species_name}"] - ) / globalVars[f"Mw_{species_name}"] + ) / mw_species[species_name] # Do volume average cell_volume, field_dict = read_cell_volumes( @@ -1190,12 +1187,12 @@ def compute_fitted_kla( # Replace all the #calc entries with their numeral values globalVars = read_global_vars(case_folder=case_folder, cross_ref=True) - # Check that global vars has the values we want and provide a useful error message otherwise + # Get Mw of the species + mw_species = {} for species_name in species_names: - if not f"Mw_{species_name}" in globalVars: - err_msg = f"Mw_{species_name} was not found in globalVars." - err_msg += "\nIf you add it, it should be [kg/mol]" - raise KeyError(err_msg) + mw_species[species_name] = _species_name_to_mw( + case_folder=case_folder, species_name=species_name + ) # Initialize the mesh fields mesh_field_dict = {} @@ -1225,8 +1222,7 @@ def compute_fitted_kla( for species_name in species_names: c_liq, kla_field_dict = compute_ave_conc_liq( time_folder=time_folder, - spec_name=species_name, - mol_weight=globalVars[f"Mw_{species_name}"], + species_name=species_name, field_dict=kla_field_dict, **kwargs, ) diff --git a/bird/preprocess/data_case_gen/loop_reactor_pbe_dynmix_nonstat_headbranch_scaleup/read_history.py b/bird/preprocess/data_case_gen/loop_reactor_pbe_dynmix_nonstat_headbranch_scaleup/read_history.py index 8d1bec1e..0ca46649 100644 --- a/bird/preprocess/data_case_gen/loop_reactor_pbe_dynmix_nonstat_headbranch_scaleup/read_history.py +++ b/bird/preprocess/data_case_gen/loop_reactor_pbe_dynmix_nonstat_headbranch_scaleup/read_history.py @@ -65,27 +65,25 @@ co2_history[itime], field_dict = compute_ave_y_liq( case_folder, time_str_sorted[itime], - spec_name="CO2", + species_name="CO2", field_dict=field_dict, ) h2_history[itime], field_dict = compute_ave_y_liq( case_folder, time_str_sorted[itime], - spec_name="H2", + species_name="H2", field_dict=field_dict, ) c_co2_history[itime], field_dict = compute_ave_conc_liq( case_folder, time_str_sorted[itime], - spec_name="CO2", - mol_weight=0.04401, + species_name="CO2", field_dict=field_dict, ) c_h2_history[itime], field_dict = compute_ave_conc_liq( case_folder, time_str_sorted[itime], - spec_name="H2", - mol_weight=0.002016, + species_name="H2", field_dict=field_dict, ) diff --git a/bird/utilities/ofio.py b/bird/utilities/ofio.py index fbbcb0f7..78b4c9e6 100644 --- a/bird/utilities/ofio.py +++ b/bird/utilities/ofio.py @@ -1191,3 +1191,78 @@ def read_global_vars( globalVars_dict = _cross_reference_global_vars(globalVars_dict) return globalVars_dict + + +def _species_name_to_mw(case_folder: str, species_name: str) -> float: + """ + Get species molecular weight from thermophysical properties or globalVars + Return species molecular weight in kg/mol + """ + thermo_gas_filename = os.path.join( + case_folder, "constant", "thermophysicalProperties.gas" + ) + thermo_liq_filename = os.path.join( + case_folder, "constant", "thermophysicalProperties.liquid" + ) + globalVars_filename = os.path.join(case_folder, "constant", "globalVars") + + assert isinstance(species_name, str) + + # Special case for water + if species_name.lower() in ["water", "h2o"]: + species_names = ["water", "h2o", "H2O", "WATER"] + else: + species_names = [ + species_name, + species_name.upper(), + species_name.lower(), + ] + + mw = None + # Try finding the molecular weight from thermo liquid + if os.path.isfile(thermo_liq_filename): + thermo = read_openfoam_dict(thermo_liq_filename) + for name in species_names: + if name in thermo: + mw = float(thermo[name]["specie"]["molWeight"]) * 1e-3 + break + if mw is not None: + logger.debug( + f"Read the {species_name} molecular weight ({mw}) from thermophysicalProperties.liquid" + ) + else: + logger.debug( + f"Could not read the {species_name} molecular weight from thermophysicalProperties.liquid" + ) + + # Try finding the molecular weight from thermo gas + if mw is None and os.path.isfile(thermo_gas_filename): + thermo = read_openfoam_dict(thermo_gas_filename) + for name in species_names: + if name in thermo: + mw = float(thermo[name]["specie"]["molWeight"]) * 1e-3 + break + if mw is not None: + logger.debug( + f"Read the {species_name} molecular weight ({mw}) from thermophysicalProperties.gas" + ) + else: + logger.debug( + f"Could not read the {species_name} molecular weight from thermophysicalProperties.gas" + ) + + # Last resort: try finding the molecular weight from thermo gas + if mw is None and os.path.isfile(globalVars_filename): + globalVars = read_global_vars(case_folder=case_folder, cross_ref=True) + for name in species_names: + if f"Mw_{name}" in globalVars: + mw = float(globalVars[f"Mw_{name}"]) + break + if mw is not None: + mw = globalVars[f"Mw_{name}"] + else: + err_msg = f"Mw_{species_name} was not found in globalVars." + err_msg += "\nIf you add it, it should be [kg/mol]" + raise KeyError(err_msg) + + return mw diff --git a/docs/source/python_interface_tut.rst b/docs/source/python_interface_tut.rst index 5cf6d15b..a1547d41 100644 --- a/docs/source/python_interface_tut.rst +++ b/docs/source/python_interface_tut.rst @@ -129,15 +129,14 @@ Several of these quantities, will require reading and processing the same fields # Compute reactor-averaged CO2 mass fraction y_ave_co2, field_dict = compute_ave_y_liq( - spec_name="CO2", field_dict=field_dict, **kwargs + species_name="CO2", field_dict=field_dict, **kwargs ) print("fields stored = ", list(field_dict.keys())) print(f"Reactor averaged YCO2 = {y_ave_co2:.4g}") # Compute reactor-averaged CO2 concentration c_ave_co2, field_dict = compute_ave_conc_liq( - spec_name="CO2", - mol_weight=0.04401, + species_name="CO2", field_dict=field_dict, **kwargs, ) diff --git a/tests/postprocess/test_post_quantities.py b/tests/postprocess/test_post_quantities.py index 812abab1..05d0f6c1 100644 --- a/tests/postprocess/test_post_quantities.py +++ b/tests/postprocess/test_post_quantities.py @@ -166,26 +166,26 @@ def test_ave_y_liq(): } field_dict = {} ave_y_co2, field_dict = compute_ave_y_liq( - spec_name="CO2", field_dict=field_dict, **kwargs + species_name="CO2", field_dict=field_dict, **kwargs ) ave_y_co, field_dict = compute_ave_y_liq( - spec_name="CO", field_dict=field_dict, **kwargs + species_name="CO", field_dict=field_dict, **kwargs ) ave_y_h2, field_dict = compute_ave_y_liq( - spec_name="H2", field_dict=field_dict, **kwargs + species_name="H2", field_dict=field_dict, **kwargs ) # Make sure None arguments are correctly handled n_cells = len(field_dict["H2.liquid"]) time_folder = kwargs["time_folder"] ave_y_h21, _ = compute_ave_y_liq( - case_folder=case_folder, time_folder=time_folder, spec_name="H2" + case_folder=case_folder, time_folder=time_folder, species_name="H2" ) ave_y_h22, _ = compute_ave_y_liq( case_folder=case_folder, n_cells=n_cells, time_folder=time_folder, - spec_name="H2", + species_name="H2", ) # Results need to be exactly the same @@ -213,20 +213,17 @@ def test_ave_conc_liq(): } field_dict = {} ave_conc_co2, field_dict = compute_ave_conc_liq( - spec_name="CO2", - mol_weight=44.00995 * 1e-3, + species_name="CO2", field_dict=field_dict, **kwargs, ) ave_conc_co, field_dict = compute_ave_conc_liq( - spec_name="CO", - mol_weight=28.01055 * 1e-3, + species_name="CO", field_dict=field_dict, **kwargs, ) ave_conc_h2, field_dict = compute_ave_conc_liq( - spec_name="H2", - mol_weight=2.01594 * 1e-3, + species_name="H2", field_dict=field_dict, **kwargs, ) @@ -236,14 +233,12 @@ def test_ave_conc_liq(): ave_conc_h21, _ = compute_ave_conc_liq( case_folder=case_folder, time_folder=time_folder, - spec_name="H2", - mol_weight=2.01594 * 1e-3, + species_name="H2", ) ave_conc_h22, _ = compute_ave_conc_liq( case_folder=case_folder, time_folder=time_folder, - spec_name="H2", - mol_weight=2.01594 * 1e-3, + species_name="H2", n_cells=n_cells, ) diff --git a/tutorial_cases/loop_reactor_mixing/read_history.py b/tutorial_cases/loop_reactor_mixing/read_history.py index f6211019..c27eae94 100644 --- a/tutorial_cases/loop_reactor_mixing/read_history.py +++ b/tutorial_cases/loop_reactor_mixing/read_history.py @@ -69,27 +69,25 @@ co2_history[itime], field_dict = compute_ave_y_liq( case_path, time_str_sorted[itime], - spec_name="CO2", + species_name="CO2", field_dict=field_dict, ) h2_history[itime], field_dict = compute_ave_y_liq( case_path, time_str_sorted[itime], - spec_name="H2", + species_name="H2", field_dict=field_dict, ) c_co2_history[itime], field_dict = compute_ave_conc_liq( case_path, time_str_sorted[itime], - spec_name="CO2", - mol_weight=0.04401, + species_name="CO2", field_dict=field_dict, ) c_h2_history[itime], field_dict = compute_ave_conc_liq( case_path, time_str_sorted[itime], - spec_name="H2", - mol_weight=0.002016, + species_name="H2", field_dict=field_dict, ) diff --git a/tutorial_cases/postprocess/python_interface_tut.py b/tutorial_cases/postprocess/python_interface_tut.py index 98bea46a..d7ab6451 100644 --- a/tutorial_cases/postprocess/python_interface_tut.py +++ b/tutorial_cases/postprocess/python_interface_tut.py @@ -62,13 +62,12 @@ print("fields stored = ", list(field_dict.keys())) print(f"Superficial velocity = {sup_vel:.4g} m/s") y_ave_co2, field_dict = compute_ave_y_liq( - spec_name="CO2", field_dict=field_dict, **kwargs + species_name="CO2", field_dict=field_dict, **kwargs ) print("fields stored = ", list(field_dict.keys())) print(f"Reactor averaged YCO2 = {y_ave_co2:.4g}") c_ave_co2, field_dict = compute_ave_conc_liq( - spec_name="CO2", - mol_weight=0.04401, + species_name="CO2", field_dict=field_dict, **kwargs, ) From 38868eb5dd685f586f16e409f15b5640dc3d5e03 Mon Sep 17 00:00:00 2001 From: Malik Date: Thu, 2 Oct 2025 08:38:55 -0600 Subject: [PATCH 11/17] provide thermo file in test case --- .../constant/thermophysicalProperties.gas | 165 ++++++++++++++++++ .../constant/thermophysicalProperties.liquid | 149 ++++++++++++++++ 2 files changed, 314 insertions(+) create mode 100644 bird/postprocess/data_conditional_mean/constant/thermophysicalProperties.gas create mode 100644 bird/postprocess/data_conditional_mean/constant/thermophysicalProperties.liquid diff --git a/bird/postprocess/data_conditional_mean/constant/thermophysicalProperties.gas b/bird/postprocess/data_conditional_mean/constant/thermophysicalProperties.gas new file mode 100644 index 00000000..66b5c9b5 --- /dev/null +++ b/bird/postprocess/data_conditional_mean/constant/thermophysicalProperties.gas @@ -0,0 +1,165 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 9 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + format ascii; + class dictionary; + object thermophysicalProperties.gas; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +thermoType +{ + type heRhoThermo; + mixture multiComponentMixture; + transport sutherland;//const;// + thermo janaf; //eConst;// + equationOfState perfectGas; + specie specie; + energy sensibleInternalEnergy; +} + +species +( + CO2 + CO + H2 + water + N2 +); + +defaultSpecie N2; + +H2 +{ + specie + { + molWeight 2.01594; + } + thermodynamics + { + Tlow 200; + Thigh 3500; + Tcommon 1000; + highCpCoeffs ( 3.3372792 -4.94024731e-05 4.99456778e-07 -1.79566394e-10 2.00255376e-14 -950.158922 -3.20502331 ); + lowCpCoeffs ( 2.34433112 0.00798052075 -1.9478151e-05 2.01572094e-08 -7.37611761e-12 -917.935173 0.683010238 ); + } + transport + { + As 6.362e-07; + Ts 72; + } + elements + { + H 2; + } +} + +CO2 +{ + specie + { + molWeight 44.00995; + } + thermodynamics + { + Tlow 200; + Thigh 3500; + Tcommon 1000; + highCpCoeffs ( 3.85746029 0.00441437026 -2.21481404e-06 5.23490188e-10 -4.72084164e-14 -48759.166 2.27163806 ); + lowCpCoeffs ( 2.35677352 0.00898459677 -7.12356269e-06 2.45919022e-09 -1.43699548e-13 -48371.9697 9.90105222 ); + } + transport + { + As 1.572e-06; + Ts 240; + } + elements + { + C 1; + O 2; + } +} + +CO +{ + specie + { + molWeight 28.01055; + } + thermodynamics + { + Tlow 200; + Thigh 3500; + Tcommon 1000; + highCpCoeffs ( 2.71518561 0.00206252743 -9.98825771e-07 2.30053008e-10 -2.03647716e-14 -14151.8724 7.81868772 ); + lowCpCoeffs ( 3.57953347 -0.00061035368 1.01681433e-06 9.07005884e-10 -9.04424499e-13 -14344.086 3.50840928 ); + } + transport + { + As 1.512e-06; + Ts 120; + } + elements + { + C 1; + O 1; + } +} + +water +{ + specie + { + molWeight 18.01534; + } + thermodynamics + { + Tlow 200; + Thigh 3500; + Tcommon 1000; + highCpCoeffs ( 3.03399249 0.00217691804 -1.64072518e-07 -9.7041987e-11 1.68200992e-14 -30004.2971 4.9667701 ); + lowCpCoeffs ( 4.19864056 -0.0020364341 6.52040211e-06 -5.48797062e-09 1.77197817e-12 -30293.7267 -0.849032208 ); + } + transport + { + As 1.512e-06; + Ts 120; + } + elements + { + H 2; + O 1; + } +} + +N2 +{ + specie + { + molWeight 28.0134; + } + thermodynamics + { + Tlow 250; + Thigh 5000; + Tcommon 1000; + highCpCoeffs ( 2.92664 0.0014879768 -5.68476e-07 1.0097038e-10 -6.753351e-15 -922.7977 5.980528 ); + lowCpCoeffs ( 3.298677 0.0014082404 -3.963222e-06 5.641515e-09 -2.444854e-12 -1020.8999 3.950372 ); + } + transport + { + As 1.512e-06; + Ts 120; + } + elements + { + N 2; + } +} +// ************************************************************************* // diff --git a/bird/postprocess/data_conditional_mean/constant/thermophysicalProperties.liquid b/bird/postprocess/data_conditional_mean/constant/thermophysicalProperties.liquid new file mode 100644 index 00000000..b13f203c --- /dev/null +++ b/bird/postprocess/data_conditional_mean/constant/thermophysicalProperties.liquid @@ -0,0 +1,149 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 9 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + format ascii; + class dictionary; + object thermophysicalProperties.liquid; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#include "$FOAM_CASE/constant/globalVars" + +thermoType +{ + type heRhoThermo; + mixture multiComponentMixture; + transport const; + thermo eConst;//hConst; + equationOfState rhoConst; + specie specie; + energy sensibleInternalEnergy; +} + +species +( + CO2 + CO + H2 + water +); + +inertSpecie water; + +"(mixture|water)" +{ + specie + { + molWeight 18.0153; + } + equationOfState + { + rho 997; + } + thermodynamics + { + Cv $CpMixLiq; + Hf -1.5879e+07; + } + transport + { + mu $muMixLiq; + Pr 2.289; + } +} + +CO2 +{ + specie + { + molWeight 44.00995; + } + equationOfState + { + rho 997; + } + thermodynamics + { + Cv $CpMixLiq; + Hf -1.5879e+07; + } + transport + { + mu $muMixLiq; + Pr $PrCO2; + } +} + +CO +{ + specie + { + molWeight 28.0106; + } + equationOfState + { + rho 997; + } + thermodynamics + { + Cv $CpMixLiq; + Hf -1.5879e+07;//-9402451; + } + transport + { + mu $muMixLiq; + Pr $PrCO; + } +} + +H2 +{ + specie + { + molWeight 2.01594; + } + equationOfState + { + rho 997; + } + thermodynamics + { + Cv $CpMixLiq; + Hf -1.5879e+07;//-9402451; + } + transport + { + mu $muMixLiq; + Pr $PrH2; + } +} + +// mixture +// { +// specie +// { +// molWeight 18; +// } +// equationOfState +// { +// rho 997; +// } +// thermodynamics +// { +// Cp 4195; +// Hf 0; +// } +// transport +// { +// mu 3.645e-4; +// Pr 2.289; +// } +// } + +// ************************************************************************* // From 3165e0b6e5032a50c627d997846f02ccae781e0e Mon Sep 17 00:00:00 2001 From: Malik Date: Thu, 2 Oct 2025 08:51:05 -0600 Subject: [PATCH 12/17] automatic reading of mol weight --- bird/postprocess/post_quantities.py | 8 ++++---- bird/utilities/ofio.py | 30 +++++++++++++++++++++++++---- 2 files changed, 30 insertions(+), 8 deletions(-) diff --git a/bird/postprocess/post_quantities.py b/bird/postprocess/post_quantities.py index f7897526..868dc83d 100644 --- a/bird/postprocess/post_quantities.py +++ b/bird/postprocess/post_quantities.py @@ -695,7 +695,7 @@ def compute_ave_y_liq( field_name="alpha.liquid", field_dict=field_dict, **kwargs ) y_liq, field_dict = read_field( - field_name=f"{spec_name}.liquid", field_dict=field_dict, **kwargs + field_name=f"{species_name}.liquid", field_dict=field_dict, **kwargs ) ind_liq, field_dict = _get_ind_liq(field_dict=field_dict, **kwargs) @@ -764,7 +764,7 @@ def compute_ave_conc_liq( if field_dict is None: field_dict = {} - mol_weight = _species_name_to_mw( + mol_weight = species_name_to_mw( case_folder=case_folder, species_name=species_name ) logger.debug( @@ -1021,7 +1021,7 @@ def compute_instantaneous_kla( err_msg = f"D_{species_name} was not found in globalVars." err_msg += f'\nIf you add it, it should be looking like #calc "1.173e-16 * pow($WC_psi * $WC_M,0.5) * $T0 / $muMixLiq / pow($WC_V_{species_name},0.6)";' raise KeyError(err_msg) - mw_species[species_name] = _species_name_to_mw( + mw_species[species_name] = species_name_to_mw( case_folder=case_folder, species_name=species_name ) @@ -1190,7 +1190,7 @@ def compute_fitted_kla( # Get Mw of the species mw_species = {} for species_name in species_names: - mw_species[species_name] = _species_name_to_mw( + mw_species[species_name] = species_name_to_mw( case_folder=case_folder, species_name=species_name ) diff --git a/bird/utilities/ofio.py b/bird/utilities/ofio.py index 78b4c9e6..e3551040 100644 --- a/bird/utilities/ofio.py +++ b/bird/utilities/ofio.py @@ -1193,10 +1193,25 @@ def read_global_vars( return globalVars_dict -def _species_name_to_mw(case_folder: str, species_name: str) -> float: - """ - Get species molecular weight from thermophysical properties or globalVars - Return species molecular weight in kg/mol +def species_name_to_mw(case_folder: str, species_name: str) -> float: + r""" + Get molecular weight in :math:`kg/mol` from the species name. + In order of availability, the molecular weight is read from the thermophysicalProperties.XXX, + and globalVars. + If nothing is found in globalVars (last resort) an error is raised. + This function is primarily useful to compute concentrations + + Parameters + ---------- + case_folder: str + Path to case folder + species_name: str + The name of the species for which molecular weight is desired + + Returns + ---------- + mw : float + The species molecular weight """ thermo_gas_filename = os.path.join( case_folder, "constant", "thermophysicalProperties.gas" @@ -1206,6 +1221,13 @@ def _species_name_to_mw(case_folder: str, species_name: str) -> float: ) globalVars_filename = os.path.join(case_folder, "constant", "globalVars") + # Handle corner case of len 1 list of strings, and interpret it as string + if ( + isinstance(species_name, list) + and len(species_name) == 1 + and isinstance(species_name[0], str) + ): + species_name = species_name[0] assert isinstance(species_name, str) # Special case for water From 2c5f6b6984ee2f146c6228ad2e0648c650adf890 Mon Sep 17 00:00:00 2001 From: Malik Date: Thu, 2 Oct 2025 09:16:16 -0600 Subject: [PATCH 13/17] adding test for molecular weight reading --- .../data_conditional_mean/constant/globalVars | 1 + .../species_gen/setup_thermo_prop.py | 28 ++++---- bird/utilities/ofio.py | 26 ++++++++ tests/io/test_chem.py | 65 +++++++++++++++++++ tests/preprocess/test_spec_thermo_gen.py | 16 ++--- 5 files changed, 116 insertions(+), 20 deletions(-) create mode 100644 tests/io/test_chem.py diff --git a/bird/postprocess/data_conditional_mean/constant/globalVars b/bird/postprocess/data_conditional_mean/constant/globalVars index 2d944d37..c1b3557d 100644 --- a/bird/postprocess/data_conditional_mean/constant/globalVars +++ b/bird/postprocess/data_conditional_mean/constant/globalVars @@ -10,6 +10,7 @@ He_CO2 0.64; Mw_CO2 0.044;//kg/mol Mw_CO 0.028; Mw_H2 0.002; +Mw_H2O 0.018; //****Liquid properties************** CpMixLiq 4181; muMixLiq 0.001; // 1 cP diff --git a/bird/preprocess/species_gen/setup_thermo_prop.py b/bird/preprocess/species_gen/setup_thermo_prop.py index 0cb0d981..5e5485e3 100644 --- a/bird/preprocess/species_gen/setup_thermo_prop.py +++ b/bird/preprocess/species_gen/setup_thermo_prop.py @@ -27,13 +27,13 @@ def check_phase_name(phase: str): raise NotImplementedError(error_msg) -def get_species_name(case_dir: str, phase: str = "gas") -> list[str]: +def get_species_name(case_folder: str, phase: str = "gas") -> list[str]: """ Get list of species name in a phase Parameters ---------- - case_dir: str + case_folder: str Path to OpenFOAM case phase: str Name of phase where to find the species @@ -47,7 +47,9 @@ def get_species_name(case_dir: str, phase: str = "gas") -> list[str]: logger.debug(f"Finding species in phase '{phase}'") thermo_properties = read_openfoam_dict( - os.path.join(case_dir, "constant", f"thermophysicalProperties.{phase}") + os.path.join( + case_folder, "constant", f"thermophysicalProperties.{phase}" + ) ) try: @@ -447,22 +449,24 @@ def update_liq_thermo_prop( return thermo_properties -def write_species_properties(case_dir: str, phase: str = "gas") -> None: +def write_species_properties(case_folder: str, phase: str = "gas") -> None: """ Write thermo properties open foam dict Parameters ---------- - case_dir: str + case_folder: str Path to OpenFOAM case phase: str Name of phase where to find the species """ - logger.debug(f"Writing properties for phase '{phase}' in case {case_dir}") - species_name = get_species_name(case_dir=case_dir, phase=phase) + logger.debug( + f"Writing properties for phase '{phase}' in case {case_folder}" + ) + species_name = get_species_name(case_folder=case_folder, phase=phase) species_prop = get_species_properties(species_name) thermo_properties_file = os.path.join( - case_dir, "constant", f"thermophysicalProperties.{phase}" + case_folder, "constant", f"thermophysicalProperties.{phase}" ) thermo_properties = read_openfoam_dict(thermo_properties_file) pair_species_keys = get_species_key_pair( @@ -477,7 +481,7 @@ def write_species_properties(case_dir: str, phase: str = "gas") -> None: thermo_properties, species_prop, pair_species_keys ) filename = os.path.join( - case_dir, "constant", f"thermophysicalProperties.{phase}" + case_folder, "constant", f"thermophysicalProperties.{phase}" ) write_openfoam_dict(thermo_properties_update, filename=filename) @@ -485,8 +489,8 @@ def write_species_properties(case_dir: str, phase: str = "gas") -> None: if __name__ == "__main__": from bird import BIRD_DIR - case_dir = os.path.join(BIRD_DIR, "../experimental_cases/deckwer17") - write_species_properties(case_dir, phase="gas") - write_species_properties(case_dir, phase="liquid") + case_folder = os.path.join(BIRD_DIR, "../experimental_cases/deckwer17") + write_species_properties(case_folder, phase="gas") + write_species_properties(case_folder, phase="liquid") # fill_global_prop(os.path.join(BIRD_DIR,"../experimental_cases_new/disengagement/bubble_column_pbe_20L/")) # fill_global_prop(os.path.join(BIRD_DIR, "../experimental_cases_new/deckwer17")) diff --git a/bird/utilities/ofio.py b/bird/utilities/ofio.py index e3551040..62da5d5a 100644 --- a/bird/utilities/ofio.py +++ b/bird/utilities/ofio.py @@ -1248,6 +1248,19 @@ def species_name_to_mw(case_folder: str, species_name: str) -> float: if name in thermo: mw = float(thermo[name]["specie"]["molWeight"]) * 1e-3 break + # Handle situation of type ("species1|species2") for the key + if mw is None: + for name in species_names: + for key in thermo: + if name in key: + try: + mw = ( + float(thermo[name]["specie"]["molWeight"]) + * 1e-3 + ) + break + except: + pass if mw is not None: logger.debug( f"Read the {species_name} molecular weight ({mw}) from thermophysicalProperties.liquid" @@ -1264,6 +1277,19 @@ def species_name_to_mw(case_folder: str, species_name: str) -> float: if name in thermo: mw = float(thermo[name]["specie"]["molWeight"]) * 1e-3 break + # Handle situation of type ("species1|species2") for the key + if mw is None: + for name in species_names: + for key in thermo: + if name in key: + try: + mw = ( + float(thermo[name]["specie"]["molWeight"]) + * 1e-3 + ) + break + except: + pass if mw is not None: logger.debug( f"Read the {species_name} molecular weight ({mw}) from thermophysicalProperties.gas" diff --git a/tests/io/test_chem.py b/tests/io/test_chem.py new file mode 100644 index 00000000..cda84820 --- /dev/null +++ b/tests/io/test_chem.py @@ -0,0 +1,65 @@ +import os +import shutil +from pathlib import Path + +import numpy as np + +from bird.preprocess.species_gen.setup_thermo_prop import ( + get_species_name, +) +from bird.utilities.ofio import species_name_to_mw + + +def test_species_mw(): + """ + Test for listing all time folders in a case + """ + case_folder = os.path.join( + Path(__file__).parent, + "..", + "..", + "bird", + "postprocess", + "data_conditional_mean", + ) + species_names = get_species_name(case_folder=case_folder, phase="liquid") + mw_species = {} + for species_name in species_names: + mw_species[species_name] = species_name_to_mw( + case_folder=case_folder, species_name=species_name + ) + assert abs(mw_species["CO2"] - 44.00995 * 1e-3) < 1e-6 + assert abs(mw_species["CO"] - 28.01055 * 1e-3) < 1e-6 + assert abs(mw_species["H2"] - 2.01594 * 1e-3) < 1e-6 + assert abs(mw_species["water"] - 18.01534 * 1e-3) < 1e-6 + + shutil.move( + os.path.join( + case_folder, "constant", "thermophysicalProperties.liquid" + ), + os.path.join(case_folder, "thermophysicalProperties.liquid_tmp"), + ) + shutil.move( + os.path.join(case_folder, "constant", "thermophysicalProperties.gas"), + os.path.join(case_folder, "thermophysicalProperties.gas_tmp"), + ) + + for species_name in species_names: + mw_species[species_name] = species_name_to_mw( + case_folder=case_folder, species_name=species_name + ) + assert abs(mw_species["CO2"] - 0.044) < 1e-6 + assert abs(mw_species["CO"] - 0.028) < 1e-6 + assert abs(mw_species["H2"] - 0.002) < 1e-6 + assert abs(mw_species["water"] - 0.018) < 1e-6 + + shutil.move( + os.path.join(case_folder, "thermophysicalProperties.liquid_tmp"), + os.path.join( + case_folder, "constant", "thermophysicalProperties.liquid" + ), + ) + shutil.move( + os.path.join(case_folder, "thermophysicalProperties.gas_tmp"), + os.path.join(case_folder, "constant", "thermophysicalProperties.gas"), + ) diff --git a/tests/preprocess/test_spec_thermo_gen.py b/tests/preprocess/test_spec_thermo_gen.py index 41bf7df7..428e4bad 100644 --- a/tests/preprocess/test_spec_thermo_gen.py +++ b/tests/preprocess/test_spec_thermo_gen.py @@ -17,7 +17,7 @@ def test_species_thermo_write(): Artificially modify the thermo properties and make sure they are written back """ - case_dir = os.path.join( + case_folder = os.path.join( Path(__file__).parent, "..", "..", @@ -28,7 +28,7 @@ def test_species_thermo_write(): # Output to temporary directory and delete when done with tempfile.TemporaryDirectory() as tmpdirname: shutil.copytree( - os.path.join(case_dir, "constant"), + os.path.join(case_folder, "constant"), os.path.join(tmpdirname, "constant"), ) @@ -68,7 +68,7 @@ def test_species_names(): Make sure the species names of all the phases can be identified """ - case_dir = os.path.join( + case_folder = os.path.join( Path(__file__).parent, "..", "..", @@ -76,14 +76,14 @@ def test_species_names(): "bubble_column_20L", ) - gas_spec_names = get_species_name(case_dir, phase="gas") + gas_spec_names = get_species_name(case_folder, phase="gas") assert len(gas_spec_names) == 3 assert "O2" in gas_spec_names assert "N2" in gas_spec_names assert "water" in gas_spec_names - liq_spec_names = get_species_name(case_dir, phase="liquid") + liq_spec_names = get_species_name(case_folder, phase="liquid") assert len(liq_spec_names) == 2 assert "O2" in liq_spec_names @@ -150,7 +150,7 @@ def test_species_key_pair(): Make sure the species names are linked to the correct openfoam dicts """ - case_dir = os.path.join( + case_folder = os.path.join( Path(__file__).parent, "..", "..", @@ -158,9 +158,9 @@ def test_species_key_pair(): "deckwer17", ) - liq_spec_names = get_species_name(case_dir, phase="liquid") + liq_spec_names = get_species_name(case_folder, phase="liquid") thermo_file = os.path.join( - case_dir, "constant", "thermophysicalProperties.liquid" + case_folder, "constant", "thermophysicalProperties.liquid" ) foam_dict = read_openfoam_dict(thermo_file) From 23171315c4c9ff373dd614a41525e2871986cefa Mon Sep 17 00:00:00 2001 From: Malik Date: Thu, 2 Oct 2025 10:01:16 -0600 Subject: [PATCH 14/17] made sure kla runs with will's case --- bird/postprocess/kla_utils.py | 8 ++++---- bird/postprocess/post_quantities.py | 17 +++++++++++------ docs/source/uq.rst | 9 ++++----- 3 files changed, 19 insertions(+), 15 deletions(-) diff --git a/bird/postprocess/kla_utils.py b/bird/postprocess/kla_utils.py index 31dcece6..7720ee27 100644 --- a/bird/postprocess/kla_utils.py +++ b/bird/postprocess/kla_utils.py @@ -365,11 +365,11 @@ def print_res_dict(res_dict: dict) -> None: bs = res_dict["bootstrapped"] if bs: - logger.info(f"\tkla = {kla:.4g} +/- {kla_err:.4g}") - logger.info(f"\tcstar = {cs:.4g} +/- {cs_err:.4g}") + logger.info(f"\tkla = {kla*3600:.4g} +/- {kla_err*3600:.4g} [h-1]") + logger.info(f"\tcstar = {cs:.4g} +/- {cs_err:.4g} [mol/m3]") logger.info(f"Without data bootstrap") - logger.info(f"\tkla = {kla_nb:.4g} +/- {kla_err_nb:.4g}") - logger.info(f"\tcstar = {cs_nb:.4g} +/- {cs_err_nb:.4g}") + logger.info(f"\tkla = {kla_nb*3600:.4g} +/- {kla_err_nb*3600:.4g} [h-1]") + logger.info(f"\tcstar = {cs_nb:.4g} +/- {cs_err_nb:.4g} [mol/m3]") if __name__ == "__main__": diff --git a/bird/postprocess/post_quantities.py b/bird/postprocess/post_quantities.py index 868dc83d..2be1d17d 100644 --- a/bird/postprocess/post_quantities.py +++ b/bird/postprocess/post_quantities.py @@ -40,7 +40,7 @@ def _field_filter( filtered_field = field else: err_msg = f"Got field type {type(field)}." - err_msg += " Expected float of np.ndarray for scalar field" + err_msg += " Expected float or np.ndarray for scalar field" raise TypeError(err_msg) elif field_type.lower() == "vector": @@ -802,7 +802,7 @@ def compute_ave_conc_liq( logger.warning( f"thermo:rho.liquid not found in {abs_time_path}, assuming it is 1000kg/m3" ) - rho_liq = 1000 + rho_liq = 1000.0 field_dict["rho_liq"] = rho_liq # Only compute over the liquid @@ -1170,7 +1170,7 @@ def compute_fitted_kla( if field_dict is None: field_dict = {} - if isinstance(species_names, "str"): + if isinstance(species_names, str): species_names = [species_names] # Read relevant fields @@ -1205,6 +1205,8 @@ def compute_fitted_kla( else: mesh_field_dict["V"], _ = read_cell_volumes(case_folder) + logger.info("Reading the species concentration history") + # Initialize the data structure for concentration c_history = {} for species_name in species_names: @@ -1212,7 +1214,7 @@ def compute_fitted_kla( # Read concentration for itime, time_folder in enumerate(time_str_sorted): - + logger.debug(f"Reading {time_folder}") # Reinitialize kla field dict kla_field_dict = {} for key in mesh_field_dict: @@ -1228,6 +1230,9 @@ def compute_fitted_kla( ) c_history[species_name][itime] = c_liq + breakpoint() + + logger.info("Doing kla fit") # Compute kla kla_spec = {} cstar_spec = {} @@ -1236,8 +1241,8 @@ def compute_fitted_kla( np.array(time_float_sorted), c_history[species_name] ) kla_spec[species_name] = { - "mean": kla_res["kla"], - "std": kla_res["kla_err"], + "mean": kla_res["kla"] * 3600, + "std": kla_res["kla_err"] * 3600, } cstar_spec[species_name] = { "mean": kla_res["cstar"], diff --git a/docs/source/uq.rst b/docs/source/uq.rst index 02c0d2e8..1661bb39 100644 --- a/docs/source/uq.rst +++ b/docs/source/uq.rst @@ -63,12 +63,11 @@ Generates scenario 2 scenario 3 For bird/postprocess/data_kla/volume_avg.dat with time index: 0, concentration index: 1 - kla = 0.09005 +/- 0.0006387 - cstar = 0.3107 +/- 0.0006122 + kla = 324 +/- 2.446 [h-1] + cstar = 0.3107 +/- 0.0006522 [mol/m3] Without data bootstrap - kla = 0.09014 +/- 0.0005957 - cstar = 0.3105 +/- 0.0005472 - + kla = 324.5 +/- 2.075 [h-1] + cstar = 0.3105 +/- 0.0005444 [mol/m3] Compute mean statistics with uncertainty ------------ From d281035ebc7c1aacc886392cf537b76f4ffd4bd5 Mon Sep 17 00:00:00 2001 From: Malik Date: Thu, 2 Oct 2025 10:04:58 -0600 Subject: [PATCH 15/17] more explicit doc string --- bird/postprocess/post_quantities.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/bird/postprocess/post_quantities.py b/bird/postprocess/post_quantities.py index 2be1d17d..943d966d 100644 --- a/bird/postprocess/post_quantities.py +++ b/bird/postprocess/post_quantities.py @@ -1159,11 +1159,11 @@ def compute_fitted_kla( kla_spec: dict Instantaneous volume averaged kLa for each species Keys are species names - Values are the kLa values + Values are dictionaries with key 'mean' (mean kLa value) and 'std' (1 standard deviation for the kLa value) cstar_spec: dict Instantaneous volume averaged cstar for each species Keys are species names - Values are the cstar values + Values are dictionaries with key 'mean' (mean cstar value) and 'std' (1 standard deviation for the cstar value) field_dict : dict Dictionary of fields read """ From f214e3e40bf7c19ccb4155e16cb4fb38acf58ad0 Mon Sep 17 00:00:00 2001 From: Malik Date: Thu, 2 Oct 2025 10:51:56 -0600 Subject: [PATCH 16/17] handle missing mu liq --- bird/postprocess/post_quantities.py | 4 +- .../species_gen/setup_thermo_prop.py | 75 +------- bird/utilities/ofio.py | 171 +++++++++++++++++- bird/version.py | 2 +- tests/io/test_chem.py | 32 +++- tests/io/test_read_foam_fields.py | 32 +++- tests/preprocess/test_spec_thermo_gen.py | 30 +-- 7 files changed, 237 insertions(+), 109 deletions(-) diff --git a/bird/postprocess/post_quantities.py b/bird/postprocess/post_quantities.py index 943d966d..a7fbf0a6 100644 --- a/bird/postprocess/post_quantities.py +++ b/bird/postprocess/post_quantities.py @@ -1045,9 +1045,7 @@ def compute_instantaneous_kla( field_name="U.liquid", field_dict=field_dict, **kwargs ) d_gas, field_dict = read_bubble_diameter(field_dict=field_dict, **kwargs) - mu_liq, field_dict = read_field( - field_name="thermo:mu.liquid", field_dict=field_dict, **kwargs - ) + mu_liq, field_dict = read_mu_liquid(field_dict=field_dict, **kwargs) species_gas = {} for species_name in species_names: species_gas[species_name], field_dict = read_field( diff --git a/bird/preprocess/species_gen/setup_thermo_prop.py b/bird/preprocess/species_gen/setup_thermo_prop.py index 5e5485e3..95ff2138 100644 --- a/bird/preprocess/species_gen/setup_thermo_prop.py +++ b/bird/preprocess/species_gen/setup_thermo_prop.py @@ -6,79 +6,14 @@ from ruamel.yaml import YAML from bird import BIRD_CONST_DIR, logger -from bird.utilities.ofio import read_openfoam_dict, write_openfoam_dict +from bird.utilities.ofio import ( + get_species_name, + read_openfoam_dict, + write_openfoam_dict, +) from bird.utilities.parser import parse_yaml -def check_phase_name(phase: str): - """ - Check that phase name is valid - - Parameters - ---------- - phase: str - Name of phase where to find the species - """ - try: - assert phase in ["gas", "liquid"] - except AssertionError: - error_msg = f"Phase name ('{phase}') is not in ['gas', 'liquid']" - logger.error(error_msg) - raise NotImplementedError(error_msg) - - -def get_species_name(case_folder: str, phase: str = "gas") -> list[str]: - """ - Get list of species name in a phase - - Parameters - ---------- - case_folder: str - Path to OpenFOAM case - phase: str - Name of phase where to find the species - - Returns - ---------- - species_name: list[str] - List of species name in the phase - """ - check_phase_name(phase) - logger.debug(f"Finding species in phase '{phase}'") - - thermo_properties = read_openfoam_dict( - os.path.join( - case_folder, "constant", f"thermophysicalProperties.{phase}" - ) - ) - - try: - species = thermo_properties["species"] - if not isinstance(species, list): - assert isinstance(species, str) - species = [species] - except KeyError: - species = [] - try: - defaultSpecie = thermo_properties["defaultSpecie"] - if not isinstance(defaultSpecie, list): - assert isinstance(defaultSpecie, str) - defaultSpecie = [defaultSpecie] - except KeyError: - defaultSpecie = [] - try: - inertSpecie = thermo_properties["inertSpecie"] - if not isinstance(inertSpecie, list): - assert isinstance(inertSpecie, str) - inertSpecie = [inertSpecie] - except KeyError: - inertSpecie = [] - - species_name = list(set(species + defaultSpecie + inertSpecie)) - logger.debug(f"Species in phase '{phase}' are {species_name}") - return species_name - - def get_species_properties(species_name: list[str]) -> dict: """ Get thermo properties of species diff --git a/bird/utilities/ofio.py b/bird/utilities/ofio.py index 62da5d5a..bb506256 100644 --- a/bird/utilities/ofio.py +++ b/bird/utilities/ofio.py @@ -7,6 +7,23 @@ from bird import logger +def _check_phase_name(phase: str): + """ + Check that phase name is valid + + Parameters + ---------- + phase: str + Name of phase where to find the species + """ + try: + assert phase in ["gas", "liquid"] + except AssertionError: + error_msg = f"Phase name ('{phase}') is not in ['gas', 'liquid']" + logger.error(error_msg) + raise NotImplementedError(error_msg) + + def _read_mesh(filename: str) -> np.ndarray: """ Reads cell center location from meshCellCentres_X.obj @@ -1037,7 +1054,7 @@ def read_bubble_diameter( d_gas = float( phaseProperties_dict["gas"]["constantCoeffs"]["d"] ) - + field_dict["d.gas"] = d_gas else: # Get field from dict if it has been read before d_gas = field_dict["d.gas"] @@ -1045,6 +1062,158 @@ def read_bubble_diameter( return d_gas, field_dict +def get_species_name(case_folder: str, phase: str = "gas") -> list[str]: + """ + Get list of species name in a phase + + Parameters + ---------- + case_folder: str + Path to OpenFOAM case + phase: str + Name of phase where to find the species + + Returns + ---------- + species_name: list[str] + List of species name in the phase + """ + _check_phase_name(phase) + logger.debug(f"Finding species in phase '{phase}'") + + thermo_properties = read_openfoam_dict( + os.path.join( + case_folder, "constant", f"thermophysicalProperties.{phase}" + ) + ) + + try: + species = thermo_properties["species"] + if not isinstance(species, list): + assert isinstance(species, str) + species = [species] + except KeyError: + species = [] + try: + defaultSpecie = thermo_properties["defaultSpecie"] + if not isinstance(defaultSpecie, list): + assert isinstance(defaultSpecie, str) + defaultSpecie = [defaultSpecie] + except KeyError: + defaultSpecie = [] + try: + inertSpecie = thermo_properties["inertSpecie"] + if not isinstance(inertSpecie, list): + assert isinstance(inertSpecie, str) + inertSpecie = [inertSpecie] + except KeyError: + inertSpecie = [] + + species_name = list(set(species + defaultSpecie + inertSpecie)) + logger.debug(f"Species in phase '{phase}' are {species_name}") + return species_name + + +def read_mu_liquid( + case_folder: str, + time_folder: str | None = None, + n_cells: int | None = None, + field_dict: dict | None = None, +) -> tuple[np.ndarray | float, dict]: + """ + Read mu_liquid at a given time and store it in dictionary for later reuse. + A specific function is constructed so that if thermo:mu.liquid is not available, the liquid viscosity is read from globalVars + + Parameters + ---------- + case_folder: str + Path to case folder + time_folder: str | None + Name of time folder to analyze. + If None, it will be found automatically + n_cells : int | None + Number of cells in the domain. + If None, it will deduced from the field reading + field_dict : dict | None + Dictionary of fields used to avoid rereading the same fields to calculate different quantities + + Returns + ---------- + cell_volumes : np.ndarray | float + Field of cell volumes + field_dict : dict + Dictionary of fields read + """ + + if field_dict is None: + field_dict = {} + + kwargs = { + "case_folder": case_folder, + "time_folder": time_folder, + "n_cells": n_cells, + } + + if ( + not ("thermo:mu.liquid" in field_dict) + or field_dict["thermo:mu.liquid"] is None + ): + try: + mu_liq, field_dict = read_field( + field_name="thermo:mu.liquid", field_dict=field_dict, **kwargs + ) + + except FileNotFoundError as err: + logger.debug( + "Could not find thermo:mu.liquid, checking if it can be read from globalVars" + ) + # thermo:mu.liquid does not exist + mu_liq = None + globalVars = read_global_vars(case_folder, cross_ref=True) + thermo = read_openfoam_dict( + os.path.join( + case_folder, "constant", "thermophysicalProperties.liquid" + ) + ) + liquid_species = get_species_name(case_folder, phase="liquid") + main_liq_species = None + for species_name in liquid_species: + if "water" in liquid_species: + main_liq_species = "water" + break + if "WATER" in liquid_species: + main_liq_species = "WATER" + break + if "h2o" in liquid_species: + main_liq_species = "h2o" + break + if "H2O" in liquid_species: + main_liq_species = "H2O" + break + for key in thermo: + if main_liq_species in key: + liq_species_prop = thermo[key] + break + if liq_species_prop["transport"]["mu"] == "$muMixLiq": + # You are using a constant mu + mu_liq = globalVars["muMixLiq"] + if mu_liq is None: + logger.error( + f"Liquid viscosity is not constant, yet could not find thermo:mu.liquid" + ) + raise FileNotFoundError(err) + else: + # Get the constant bubble diameter + logger.debug("Reading liquid viscosity from globalVars") + field_dict["thermo:mu.liquid"] = mu_liq + + else: + # Get field from dict if it has been read before + mu_liq = field_dict["thermo:mu.liquid"] + + return mu_liq, field_dict + + def _cross_reference_global_vars( globalVars_dict: dict, ) -> dict: diff --git a/bird/version.py b/bird/version.py index 0461d1e7..9154fcd5 100644 --- a/bird/version.py +++ b/bird/version.py @@ -1,3 +1,3 @@ """Bio reactor design version""" -__version__ = "0.0.42" +__version__ = "0.0.43" diff --git a/tests/io/test_chem.py b/tests/io/test_chem.py index cda84820..d366bfef 100644 --- a/tests/io/test_chem.py +++ b/tests/io/test_chem.py @@ -4,10 +4,7 @@ import numpy as np -from bird.preprocess.species_gen.setup_thermo_prop import ( - get_species_name, -) -from bird.utilities.ofio import species_name_to_mw +from bird.utilities.ofio import get_species_name, species_name_to_mw def test_species_mw(): @@ -63,3 +60,30 @@ def test_species_mw(): os.path.join(case_folder, "thermophysicalProperties.gas_tmp"), os.path.join(case_folder, "constant", "thermophysicalProperties.gas"), ) + + +def test_species_names(): + """ + Make sure the species names of all the phases can be identified + """ + + case_folder = os.path.join( + Path(__file__).parent, + "..", + "..", + "tutorial_cases", + "bubble_column_20L", + ) + + gas_spec_names = get_species_name(case_folder, phase="gas") + + assert len(gas_spec_names) == 3 + assert "O2" in gas_spec_names + assert "N2" in gas_spec_names + assert "water" in gas_spec_names + + liq_spec_names = get_species_name(case_folder, phase="liquid") + + assert len(liq_spec_names) == 2 + assert "O2" in liq_spec_names + assert "water" in liq_spec_names diff --git a/tests/io/test_read_foam_fields.py b/tests/io/test_read_foam_fields.py index c0beed96..44afb59e 100644 --- a/tests/io/test_read_foam_fields.py +++ b/tests/io/test_read_foam_fields.py @@ -9,6 +9,7 @@ _readOFScal, _readOFVec, read_bubble_diameter, + read_mu_liquid, ) @@ -184,5 +185,34 @@ def test_read_bubble_diameter(): assert abs(d_gas[2] - 0.00966694) < 1e-8 +def test_read_mu_liquid(): + """ + Test for reading bubble diameter + """ + case_folder = os.path.join( + Path(__file__).parent, + "..", + "..", + "bird", + "postprocess", + "data_conditional_mean", + ) + # Read bubble diameter + mu_liq, _ = read_mu_liquid(case_folder=case_folder, time_folder="80") + assert abs(mu_liq - 0.001) < 1e-8 + # Read bubble diameter if no d.gas is available + # Manufacture a case without d.gas + shutil.move( + os.path.join(case_folder, "80", "thermo:mu.liquid"), + os.path.join(case_folder, "thermo:mu.liquid_tmp"), + ) + mu_liq, _ = read_mu_liquid(case_folder=case_folder, time_folder="80") + assert abs(mu_liq - 0.001) < 1e-8 + shutil.move( + os.path.join(case_folder, "thermo:mu.liquid_tmp"), + os.path.join(case_folder, "80", "thermo:mu.liquid"), + ) + + if __name__ == "__main__": - test_read_bubble_diameter() + test_read_mu_liquid() diff --git a/tests/preprocess/test_spec_thermo_gen.py b/tests/preprocess/test_spec_thermo_gen.py index 428e4bad..8f9fd932 100644 --- a/tests/preprocess/test_spec_thermo_gen.py +++ b/tests/preprocess/test_spec_thermo_gen.py @@ -5,11 +5,10 @@ from bird.preprocess.species_gen.setup_thermo_prop import ( get_species_key_pair, - get_species_name, get_species_properties, write_species_properties, ) -from bird.utilities.ofio import read_openfoam_dict +from bird.utilities.ofio import get_species_name, read_openfoam_dict def test_species_thermo_write(): @@ -63,33 +62,6 @@ def test_species_thermo_write(): ) -def test_species_names(): - """ - Make sure the species names of all the phases can be identified - """ - - case_folder = os.path.join( - Path(__file__).parent, - "..", - "..", - "tutorial_cases", - "bubble_column_20L", - ) - - gas_spec_names = get_species_name(case_folder, phase="gas") - - assert len(gas_spec_names) == 3 - assert "O2" in gas_spec_names - assert "N2" in gas_spec_names - assert "water" in gas_spec_names - - liq_spec_names = get_species_name(case_folder, phase="liquid") - - assert len(liq_spec_names) == 2 - assert "O2" in liq_spec_names - assert "water" in liq_spec_names - - def test_species_prop(): """ Make sure the species properties are correctly read From fc63a38ffc51958732643fca767f0b24b3c0dd14 Mon Sep 17 00:00:00 2001 From: Malik Date: Thu, 2 Oct 2025 10:54:19 -0600 Subject: [PATCH 17/17] fix docstring --- bird/utilities/ofio.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/bird/utilities/ofio.py b/bird/utilities/ofio.py index bb506256..776c1204 100644 --- a/bird/utilities/ofio.py +++ b/bird/utilities/ofio.py @@ -1121,7 +1121,7 @@ def read_mu_liquid( field_dict: dict | None = None, ) -> tuple[np.ndarray | float, dict]: """ - Read mu_liquid at a given time and store it in dictionary for later reuse. + Read liquid viscosity at a given time and store it in dictionary for later reuse. A specific function is constructed so that if thermo:mu.liquid is not available, the liquid viscosity is read from globalVars Parameters