Better way of avoiding ill-posed problem.

modelica · Jan 21, 2025 · b010701 · b010701
1 parent caf3eb7
commit b010701
Showing 1 changed file with 12 additions and 16 deletions.
diff --git a/Modelica/Fluid/Utilities.mo b/Modelica/Fluid/Utilities.mo
@@ -691,7 +691,7 @@ for a smooth transition from y1 to y2.
           c := 0.1*Delta0;
         end if;
         theta0 := (y0d*mu + y1d*eta)/h0;
-        if abs(theta0 - c) < 1e-6 then
+        if abs(theta0 - c) < 1e-6*abs(c) then
           // Slightly reduce c in order to avoid ill-posed problem
           c := (1 - 1e-6)*theta0;
         end if;
@@ -700,22 +700,18 @@ for a smooth transition from y1 to y2.
         eta_tilde := rho*eta;
         xi1 := x0 + mu_tilde;
         xi2 := x1 - eta_tilde;
-        if  xi1 < x0 or xi2>x1 or abs(xi1-xi2) < 100*Modelica.Constants.eps then
-          // The limits don't make sense, just use linear interpolation
-          y := (y1-y0)*(x-x0)/(x1-x0) + y0;
+
+        a1 := (y0d - c)/max(mu_tilde^2, 100*Modelica.Constants.eps);
+        a2 := (y1d - c)/max(eta_tilde^2, 100*Modelica.Constants.eps);
+        const12 := y0 - a1/3*(x0 - xi1)^3 - c*x0;
+        const3 := y1 - a2/3*(x1 - xi2)^3 - c*x1;
+        // Do actual interpolation
+        if (x < xi1) then
+          y := a1/3*(x - xi1)^3 + c*x + const12;
+        elseif (x < xi2) then
+          y := c*x + const12;
         else
-          a1 := (y0d - c)/max(mu_tilde^2, 100*Modelica.Constants.eps);
-          a2 := (y1d - c)/max(eta_tilde^2, 100*Modelica.Constants.eps);
-          const12 := y0 - a1/3*(x0 - xi1)^3 - c*x0;
-          const3 := y1 - a2/3*(x1 - xi2)^3 - c*x1;
-          // Do actual interpolation
-          if (x < xi1) then
-            y := a1/3*(x - xi1)^3 + c*x + const12;
-          elseif (x < xi2) then
-            y := c*x + const12;
-          else
-            y := a2/3*(x - xi2)^3 + c*x + const3;
-          end if;
+          y := a2/3*(x - xi2)^3 + c*x + const3;
         end if;
       else
         // Cubic S0 is monotonic, use it as is