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/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/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/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/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; +// } +// } + +// ************************************************************************* // diff --git a/bird/postprocess/kla_utils.py b/bird/postprocess/kla_utils.py index 28095d0d..7720ee27 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,15 +364,12 @@ 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}") + 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 15a91515..a7fbf0a6 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 @@ -38,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": @@ -634,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""" @@ -663,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 @@ -693,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) @@ -717,11 +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, - rho_val: float | None = 1000, field_dict: dict | None = None, ) -> tuple[float, dict]: r""" @@ -743,18 +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) - 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 @@ -768,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 = { @@ -790,25 +786,24 @@ 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) 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.0 + field_dict["rho_liq"] = rho_liq # Only compute over the liquid alpha_liq = _field_filter(alpha_liq, ind=ind_liq, field_type="scalar") @@ -884,9 +879,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( - 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( @@ -909,7 +902,7 @@ 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, field_dict: dict | None = None, @@ -971,7 +964,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. @@ -998,6 +991,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, @@ -1015,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) @@ -1048,12 +1044,8 @@ 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 - ) - mu_liq, field_dict = read_field( - field_name="thermo:mu.liquid", field_dict=field_dict, **kwargs - ) + d_gas, field_dict = read_bubble_diameter(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( @@ -1099,7 +1091,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( @@ -1118,3 +1110,141 @@ 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 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 dictionaries with key 'mean' (mean cstar value) and 'std' (1 standard deviation for the cstar value) + 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) + + # Get Mw of the species + mw_species = {} + for species_name in species_names: + mw_species[species_name] = species_name_to_mw( + case_folder=case_folder, species_name=species_name + ) + + # 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) + + logger.info("Reading the species concentration history") + + # 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): + logger.debug(f"Reading {time_folder}") + # 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, + species_name=species_name, + field_dict=kla_field_dict, + **kwargs, + ) + c_history[species_name][itime] = c_liq + + breakpoint() + + logger.info("Doing kla fit") + # 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"] * 3600, + "std": kla_res["kla_err"] * 3600, + } + cstar_spec[species_name] = { + "mean": kla_res["cstar"], + "std": kla_res["cstar_err"], + } + + return kla_spec, cstar_spec, field_dict 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/preprocess/species_gen/setup_thermo_prop.py b/bird/preprocess/species_gen/setup_thermo_prop.py index 0cb0d981..95ff2138 100644 --- a/bird/preprocess/species_gen/setup_thermo_prop.py +++ b/bird/preprocess/species_gen/setup_thermo_prop.py @@ -6,77 +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_dir: str, phase: str = "gas") -> list[str]: - """ - Get list of species name in a phase - - Parameters - ---------- - case_dir: 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_dir, "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 @@ -447,22 +384,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 +416,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 +424,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 5a72dbc2..776c1204 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 @@ -967,6 +984,236 @@ 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"] + ) + field_dict["d.gas"] = d_gas + else: + # Get field from dict if it has been read before + d_gas = field_dict["d.gas"] + + 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 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 + ---------- + 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: @@ -1113,3 +1360,126 @@ 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: + 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" + ) + thermo_liq_filename = os.path.join( + case_folder, "constant", "thermophysicalProperties.liquid" + ) + 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 + 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 + # 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" + ) + 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 + # 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" + ) + 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/bird/version.py b/bird/version.py index 7b1b5e3f..9154fcd5 100644 --- a/bird/version.py +++ b/bird/version.py @@ -1,3 +1,3 @@ """Bio reactor design version""" -__version__ = "0.0.40" +__version__ = "0.0.43" diff --git a/docs/source/python_interface_tut.rst b/docs/source/python_interface_tut.rst index 69e05dbd..a1547d41 100644 --- a/docs/source/python_interface_tut.rst +++ b/docs/source/python_interface_tut.rst @@ -129,16 +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, - rho_val=1000, + species_name="CO2", field_dict=field_dict, **kwargs, ) 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 ------------ 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/tests/io/test_chem.py b/tests/io/test_chem.py new file mode 100644 index 00000000..d366bfef --- /dev/null +++ b/tests/io/test_chem.py @@ -0,0 +1,89 @@ +import os +import shutil +from pathlib import Path + +import numpy as np + +from bird.utilities.ofio import get_species_name, 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"), + ) + + +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 a991b8d2..44afb59e 100644 --- a/tests/io/test_read_foam_fields.py +++ b/tests/io/test_read_foam_fields.py @@ -1,9 +1,16 @@ 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, + read_mu_liquid, +) def test_read_nonunif_scal(): @@ -132,3 +139,80 @@ 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 + + +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_mu_liquid() 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) diff --git a/tests/postprocess/test_post_quantities.py b/tests/postprocess/test_post_quantities.py index 904a5bc3..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, ) @@ -307,6 +302,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__": diff --git a/tests/preprocess/test_spec_thermo_gen.py b/tests/preprocess/test_spec_thermo_gen.py index 41bf7df7..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(): @@ -17,7 +16,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 +27,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"), ) @@ -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_dir = os.path.join( - Path(__file__).parent, - "..", - "..", - "tutorial_cases", - "bubble_column_20L", - ) - - gas_spec_names = get_species_name(case_dir, 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") - - 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 @@ -150,7 +122,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 +130,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) 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/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/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/postprocess/python_interface_tut.py b/tutorial_cases/postprocess/python_interface_tut.py index 085713db..d7ab6451 100644 --- a/tutorial_cases/postprocess/python_interface_tut.py +++ b/tutorial_cases/postprocess/python_interface_tut.py @@ -62,14 +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, - rho_val=1000, + species_name="CO2", field_dict=field_dict, **kwargs, ) 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