From 63ceed9817ff36c8e27c44c8e2c1d26957e1c48f Mon Sep 17 00:00:00 2001 From: Vegard Gjeldvik Jervell Date: Tue, 12 Nov 2024 01:54:32 +0100 Subject: [PATCH] Fix stupid Soret coefficient bug and throw if computation of not-implemented stuff is attempted. --- cpp/KineticGas_properties.cpp | 5 ++- tests/test_binary_limit.py | 60 ++++++++++++++++++++++++++++++++++- 2 files changed, 63 insertions(+), 2 deletions(-) diff --git a/cpp/KineticGas_properties.cpp b/cpp/KineticGas_properties.cpp index f0cb1ac..a0779ae 100644 --- a/cpp/KineticGas_properties.cpp +++ b/cpp/KineticGas_properties.cpp @@ -175,13 +175,14 @@ Eigen::MatrixXd KineticGas::interdiffusion(double T, double Vm, const vector1d& if (dependent_idx < 0 && frame_of_reference == FrameOfReference::solvent) dependent_idx = solvent_idx; while (dependent_idx < 0) dependent_idx += Ncomps; if (frame_of_reference == FrameOfReference::zarate_x){ + if (!is_idealgas) throw std::runtime_error("Zarate relation for diffusion only implemented for ideal gas!"); Eigen::MatrixXd D = interdiffusion(T, Vm, x, N, FrameOfReference::CoN, dependent_idx, dependent_idx, true); return D; } else if (frame_of_reference == FrameOfReference::zarate){ Eigen::MatrixXd X = get_zarate_X_matr(x, dependent_idx); Eigen::MatrixXd Dx = interdiffusion(T, Vm, x, N, FrameOfReference::zarate_x, dependent_idx); - return X * Dx * X.inverse(); + return X.inverse() * Dx * X; } else if (frame_of_reference == FrameOfReference::zarate_w){ Eigen::MatrixXd W = get_zarate_W_matr(x, dependent_idx); @@ -257,6 +258,7 @@ Eigen::VectorXd KineticGas::thermal_diffusion_coeff(double T, double Vm, const v DT /= AVOGADRO; if (use_zarate){ + if (!is_idealgas) throw std::runtime_error("Zarate relation for thermal diffusion only implemented for ideal gas!"); Eigen::VectorXd DT_indep(Ncomps - 1); Eigen::MatrixXd D_indep = interdiffusion(T, Vm, x, N, frame_of_reference, dependent_idx, solvent_idx); @@ -319,6 +321,7 @@ Eigen::MatrixXd KineticGas::thermal_diffusion_factor(double T, double Vm, const } Eigen::VectorXd KineticGas::soret_coefficient(double T, double Vm, const std::vector& x, int N, int dependent_idx){ + if (!is_idealgas) throw std::runtime_error("Zarate relation for Soret coefficient only implemented for ideal gas!"); Eigen::MatrixXd D = interdiffusion(T, Vm, x, N, FrameOfReference::zarate, dependent_idx, dependent_idx, false); Eigen::VectorXd DT = thermal_diffusion_coeff(T, Vm, x, N, FrameOfReference::zarate, dependent_idx); return D.partialPivLu().solve(DT); diff --git a/tests/test_binary_limit.py b/tests/test_binary_limit.py index 01fcaec..f446d32 100644 --- a/tests/test_binary_limit.py +++ b/tests/test_binary_limit.py @@ -117,5 +117,63 @@ def test_th_diffusion(model, N, silent=True): assert all((diff / DT_b) < 1e-3) +@pytest.mark.parametrize('model', models) +@pytest.mark.parametrize('N', [2, 3]) +def test_soret(model, N, silent=True): + T = 300 + Vm = 1 / 10 + x_bin = [0.3, 0.7] + bin_12 = model('H2O,O2', is_idealgas=True) + ternary = model('H2O,O2,N2', is_idealgas=True) + + x_bin = [0.6, 0.4] + ST_1b2 = bin_12.soret_coefficient(T, Vm, x_bin, N=N, use_zarate=True)[0] + ST_1t, ST_2t = ternary.soret_coefficient(T, Vm, get_ternary_x(x_bin, 0.5), N=N, use_zarate=True) + diff = abs((ST_1t - ST_2t) - ST_1b2) + x3_lst = np.logspace(-1, -5, 10) + + print(ST_1t - ST_2t, np.diff([ST_1t, ST_2t]), ST_1b2) + + if silent is False: + plt.plot(x3_lst, [- np.diff(ternary.soret_coefficient(T, Vm, get_ternary_x(x_bin, x3), N=N, use_zarate=True)) for x3 in x3_lst]) + plt.plot(x3_lst, ST_1b2 * np.ones_like(x3_lst)) + plt.xscale('log') + plt.show() + + for x3 in x3_lst: + ST_1t, ST_2t = ternary.soret_coefficient(T, Vm, get_ternary_x(x_bin, x3), N=N, use_zarate=True) + new_diff = abs((ST_1t - ST_2t) - ST_1b2) + assert new_diff < diff + diff = new_diff + + bin_23 = model('O2,N2', is_idealgas=True) + def get_ternary_x1(x1): + return [x1, x_bin[0] * (1 - x1), x_bin[1] * (1 - x1)] + + ST_2b3 = bin_23.soret_coefficient(T, Vm, x_bin, N=N, use_zarate=True)[0] + _, ST_2t = ternary.soret_coefficient(T, Vm, get_ternary_x1(0.11), N=N, use_zarate=True) + diff = abs(ST_2t - ST_2b3) + x1_lst = np.logspace(-1, -5, 10) + for x1 in x1_lst: + _, ST_2t = ternary.soret_coefficient(T, Vm, get_ternary_x1(x1), N=N, use_zarate=True) + new_diff = abs(ST_2t - ST_2b3) + assert new_diff < diff + diff = new_diff + + bin_13 = model('H2O,N2', is_idealgas=True) + def get_ternary_x2(x2): + return [x_bin[0] * (1 - x2), x2, x_bin[1] * (1 - x2)] + + ST_1b3 = bin_13.soret_coefficient(T, Vm, x_bin, N=N, use_zarate=True)[0] + ST_1t, _ = ternary.soret_coefficient(T, Vm, get_ternary_x2(0.5), N=N, use_zarate=True) + diff = abs(ST_1t - ST_1b3) + x2_lst = np.logspace(-1, -5, 10) + for x2 in x2_lst: + ST_1t, _ = ternary.soret_coefficient(T, Vm, get_ternary_x2(x2), N=N, use_zarate=True) + new_diff = abs(ST_1t - ST_1b3) + assert new_diff < diff + diff = new_diff + + if __name__ == "__main__": - test_diffusion(models[0], 3) \ No newline at end of file + test_soret(models[0], 2, silent=True) \ No newline at end of file