Skip to content

Commit

Permalink
change saturation convergence criterion:
Browse files Browse the repository at this point in the history 
do not use RSS, instead keep iterating as long as residual OR Newton step are larger than tolerance
  • Loading branch information
@thorade
thorade committed May 3, 2020
1 parent 7753831 commit 89c6b9b
Show file tree
Hide file tree
Showing 3 changed files with 10 additions and 21 deletions.
6 changes: 3 additions & 3 deletions HelmholtzMedia/Examples/ConvergenceTest/setSat.mo
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,8 +12,8 @@ model setSat

Modelica.SIunits.TemperatureDifference T_err_abs;
Modelica.SIunits.PressureDifference p_err_abs;
Modelica.SIunits.Density dl_err_abs;
Modelica.SIunits.Density dv_err_abs;
Modelica.SIunits.Density dl_err_abs(min=-1e3) "Liquid density difference";
Modelica.SIunits.Density dv_err_abs(min=-1e3) "Vapour density difference";

Modelica.Blocks.Sources.Ramp T_ramp(
duration=8,
Expand Down Expand Up @@ -43,7 +43,7 @@ equation
T_err_abs = sat_T.Tsat - sat_p.Tsat;
p_err_abs = sat_T.psat - sat_p.psat;
dl_err_abs = sat_T.liq.d - sat_dl.liq.d;
dv_err_abs = sat_T.vap.d - sat_dl.vap.d;
dv_err_abs = sat_T.vap.d - sat_dv.vap.d;

annotation (experiment(
StopTime=10,
Expand Down
14 changes: 4 additions & 10 deletions HelmholtzMedia/Interfaces/PartialHelmholtzMedium/package.mo
View file
Original file line number Diff line number Diff line change
Expand Up @@ -151,12 +151,11 @@ protected
Real K_vap_delta;

Real RES[2] "residual function vector";
Real RSS "residual sum of squares";
Real Jacobian[2,2] "Jacobian matrix";
Real NS[2] "Newton step vector";

constant Real lambda(min=0.1,max=1) = 1 "convergence speed, default=1";
constant Real tolerance=1e-18 "tolerance for RSS";
constant Real tolerance=1e-7 "tolerance for p and d, needs to be smaller than simulation tolerance";
Integer iter = 0;
constant Integer iter_max = 200;

Expand All @@ -183,9 +182,8 @@ protected
K_vap := delta_vap*fv.rd + fv.r + log(delta_vap);
// residual vector
RES := {(J_vap-J_liq), (K_vap-K_liq)};
RSS := RES*RES/2;

while (RSS>tolerance) and (iter<iter_max) loop
while (iter<iter_max) and (iter<1 or abs(RES[1])>tolerance or abs(RES[2])>tolerance or abs(NS[1])>tolerance or abs(NS[2])>tolerance) loop
iter := iter+1;

// calculate gradients of J and K, set up Jacobian matrix, get Newton step
Expand Down Expand Up @@ -223,7 +221,6 @@ protected
K_vap := delta_vap*fv.rd + fv.r + log(delta_vap);
// residual vector
RES := {(J_vap-J_liq), (K_vap-K_liq)};
RSS := RES*RES/2;

end while;
// Modelica.Utilities.Streams.print("setSat_T total iteration steps " + String(iter), "printlog.txt");
Expand Down Expand Up @@ -297,12 +294,11 @@ protected
EoS.HelmholtzDerivs fv;

Real RES[3] "residual function vector";
Real RSS "residual sum of squares";
Real Jacobian[3,3] "Jacobian matrix";
Real NS[3] "Newton step vector";

constant Real lambda(min=0.1,max=1) = 1 "convergence speed, default=1";
constant Real tolerance=1e-9 "tolerance for RSS";
constant Real tolerance=1e-6 "tolerance for T and d, needs to be smaller than simulation tolerance";
Integer iter = 0;
constant Integer iter_max = 200;

Expand All @@ -321,9 +317,8 @@ protected
fl := EoS.setHelmholtzDerivsSecond(d=sat.liq.d, T=sat.Tsat, phase=1);
fv := EoS.setHelmholtzDerivsSecond(d=sat.vap.d, T=sat.Tsat, phase=1);
RES := {EoS.p(fl)-p, EoS.p(fv)-p, EoS.g(fl)-EoS.g(fv)};
RSS := RES*RES/2;

while (RSS>tolerance) and (iter<iter_max) loop
while (iter<iter_max) and (iter<1 or abs(RES[1])>tolerance or abs(RES[2])>tolerance or abs(NS[1])>tolerance or abs(NS[2])>tolerance or abs(NS[3])>tolerance) loop
iter := iter+1;

// calculate Jacobian matrix and Newton Step vector
Expand All @@ -349,7 +344,6 @@ protected
fl := EoS.setHelmholtzDerivsSecond(d=sat.liq.d, T=sat.Tsat, phase=1);
fv := EoS.setHelmholtzDerivsSecond(d=sat.vap.d, T=sat.Tsat, phase=1);
RES := {EoS.p(fl)-p, EoS.p(fv)-p, EoS.g(fl)-EoS.g(fv)};
RSS := RES*RES/2;
end while;
// if verbose then Modelica.Utilities.Streams.print("setSat_p total iteration steps " + String(iter), "printlog.txt"); end if;
// Modelica.Utilities.Streams.print("setSat_p total iteration steps " + String(iter), "printlog.txt");
Expand Down
11 changes: 3 additions & 8 deletions HelmholtzMedia/Interfaces/PartialHelmholtzMedium/setSat_d.mo
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,12 +20,11 @@ protected
EoS.HelmholtzDerivs fv;

Real RES[2] "residual function vector";
Real RSS "residual sum of squares";
Real Jacobian[2,2] "Jacobian matrix";
Real NS[2] "Newton step vector";

constant Real lambda(min=0.1,max=1) = 1 "convergence speed, default=1";
constant Real tolerance=1e-9 "tolerance for RSS";
constant Real tolerance=1e-6 "tolerance for p and T and d, needs to be smaller than simulation tolerance";
Integer iter = 0;
constant Integer iter_max = 200;

Expand All @@ -43,9 +42,8 @@ algorithm
fl := EoS.setHelmholtzDerivsSecond(d=sat.liq.d, T=sat.Tsat, phase=1);
fv := EoS.setHelmholtzDerivsSecond(d=sat.vap.d, T=sat.Tsat, phase=1);
RES := {EoS.p(fl)-EoS.p(fv), EoS.g(fl)-EoS.g(fv)};
RSS := RES*RES/2;

while (RSS>tolerance) and (iter<iter_max) loop
while (iter<iter_max) and (iter<1 or abs(RES[1])>tolerance or abs(RES[2])>tolerance or abs(NS[1])>tolerance or abs(NS[2])>tolerance) loop
iter := iter+1;

// calculate Jacobian matrix and Newton Step vector
Expand All @@ -67,7 +65,6 @@ algorithm
fl := EoS.setHelmholtzDerivsSecond(d=sat.liq.d, T=sat.Tsat, phase=1);
fv := EoS.setHelmholtzDerivsSecond(d=sat.vap.d, T=sat.Tsat, phase=1);
RES := {EoS.p(fl)-EoS.p(fv), EoS.g(fl)-EoS.g(fv)};
RSS := RES*RES/2;
end while;
sat.liq := setState_dTX(d=sat.liq.d, T=sat.Tsat, phase=1);
sat.vap := setState_dTX(d=sat.vap.d, T=sat.Tsat, phase=1);
Expand All @@ -82,9 +79,8 @@ algorithm
fv := EoS.setHelmholtzDerivsSecond(d=sat.vap.d, T=sat.Tsat, phase=1);
fl := EoS.setHelmholtzDerivsSecond(d=sat.liq.d, T=sat.Tsat, phase=1);
RES := {EoS.p(fv)-EoS.p(fl), EoS.g(fv)-EoS.g(fl)};
RSS := RES*RES/2;

while (RSS>tolerance) and (iter<iter_max) loop
while (iter<iter_max) and (iter<1 or abs(RES[1])>tolerance or abs(RES[2])>tolerance or abs(NS[1])>tolerance or abs(NS[2])>tolerance) loop
iter := iter+1;

// calculate Jacobian matrix and Newton Step vector
Expand All @@ -106,7 +102,6 @@ algorithm
fv := EoS.setHelmholtzDerivsSecond(d=sat.vap.d, T=sat.Tsat, phase=1);
fl := EoS.setHelmholtzDerivsSecond(d=sat.liq.d, T=sat.Tsat, phase=1);
RES := {EoS.p(fv)-EoS.p(fl), EoS.g(fv)-EoS.g(fl)};
RSS := RES*RES/2;
end while;
sat.liq := setState_dTX(d=sat.liq.d, T=sat.Tsat, phase=1);
sat.vap := setState_dTX(d=sat.vap.d, T=sat.Tsat, phase=1);
Expand Down

2 comments on commit 89c6b9b

@thorade
Copy link
Owner Author

@thorade thorade commented on 89c6b9b May 3, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

for #34

@thorade
Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

also see #40

Please sign in to comment.