Fix stupid Soret coefficient bug and throw if computation of not-impl…

…emented stuff is attempted.

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<double>& 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);

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)
+    test_soret(models[0], 2, silent=True)