Merge pull request #948 from slacgismo/develop

Release beauharnois-25

.gitignore

@@ -83,6 +83,7 @@ Makefile
 config.log
 config.status
 configure
+configure~
 gldcore/gridlabd
 gridlabd.spec
 libtool
@@ -173,3 +174,4 @@ gldcore/converters/autotest/table2glm_input_noclass.glm
 gldcore/converters/autotest/table2glm_input_noname.glm
 gldcore/converters/autotest/table2glm_input.glm
 gldcore/geodata/docs/geodata.conf
+gldcore/gridlabd.in

docs/GLM/General/Inherit.md

@@ -0,0 +1,38 @@
+[[/GLM/General/Inherit.md]] -- Inherit property values
+
+# Synopsis
+
+GLM:
+
+~~~
+object <class>
+{
+    <property-name> inherit;
+    <property_name> @<object-name>;
+}
+~~~
+
+# Description
+
+The `inherit` property value causes the value to be copied from the parent object. This requires that the parent already be defined.  The alternative syntax `@<name>` inherits the property value from the named object.
+
+# Example
+
+This example copies the value of `x` from object `test1` into object `test2`:
+
+~~~
+class test
+{
+    double x;
+}
+object test
+{
+    name test1;
+    x 1.23;
+}
+object test
+{
+    name test2;
+    x @test1;
+}
+~~~

docs/Subcommand/Json-get.md

@@ -5,13 +5,15 @@
 GLM:
 
 ~~~
-#json-get [keys ...] [-k|--keys] [-j|--json|-r|--raw|-r|--csv] <INPUT >OUTPUT 
+#json-get [keys ...] [-k|--keys] [-j|--json|-r|--raw|-c|--csv] [-i|--input INPUT]
+ [-o|--output OUTPUT]
 ~~~
 
 Shell:
 
 ~~~
-$ gridlabd json-get [keys ...] [-k|--keys] [-j|--json|-r|--raw|-r|--csv] [<INPUT] [>OUTPUT] 
+$ gridlabd json-get [keys ...] [-k|--keys] [-j|--json|-r|--raw|-r|--csv] [-i|--input INPUT]
+ [-o|--output OUTPUT] 
 ~~~
 
 # Description

gldcore/autotest/test_inherit.glm

@@ -0,0 +1,21 @@
+module assert;
+
+class test
+ {
+     double x;
+ }
+ object test
+ {
+     name test1;
+     x 8.88;
+ }
+ object test
+ {
+     name test2;
+     x @test1;
+     object assert{
+        target x;
+        relation "==";
+        value 8.88;
+     };
+ }

gldcore/cmdarg.cpp

@@ -2249,7 +2249,7 @@ DEPRECATED static CMDARG main_commands[] = {
 	{"pidfile",		NULL,	pidfile,		"[=<filename>]", "Set the process ID file (default is gridlabd.pid)" },
 	{"threadcount", "T",	threadcount,	"<n>", "Set the maximum number of threads allowed" },
 	{"job",			NULL,	job,			"...", "Start a job"},
-    {"nprocs",      NULL,   nprocs,         "Display the number of processors available to run jobs"},
+    {"nprocs",      NULL,   nprocs,         NULL, "Display the number of processors available to run jobs"},
 
 	{NULL,NULL,NULL,NULL, "System options"},
 	{"avlbalance",	NULL,	avlbalance,		NULL, "Toggles automatic balancing of object index" },

gldcore/converters/json2html.py

@@ -146,7 +146,7 @@ def help():
             if not icon:
                 warning(f"object '{name}' has no known icon (class '{oclass})'")
                 icon = folium.Icon(color=color)
-            elif pos:
+            if pos:
                 obj = folium.Marker(pos,icon=icon,popup=popup,name=name)
         obj.add_to(cluster)
     if mouseposition:

gldcore/geodata/geodata_powerline.py

@@ -163,21 +163,12 @@ def get_air_properties(temp):
     """
     # calculate the constant parameters
     # reference: http://home.eng.iastate.edu/~jdm/wind/TransmissionLineLoadingDesignCriteriaAndHTS.pdf
-    temp_x = np.arange(0.0, 110.0, 10) # unit: DegC
-    specific_mass = [1.29, 1.25, 1.20, 1.17, 1.13, 1.09, 1.06, 1.03, 1.0, 0.97, 0.95] # unit: kg/m3
-    thermal_conductivity = [0.0243, 0.025, 0.0257, 0.0265, 0.0272, 0.028, 0.0287,0.0294, 0.0301, 0.0309, 0.0316] # unit: W/(m*DegC)
-    dynamic_viscosity = [0.175e-4, 0.18e-4, 0.184e-4, 0.189e-4, 0.194e-4, 0.199e-4, 0.203e-4, 0.208e-4, 0.213e-4, 0.217e-4, 0.222e-4] # unit: N*s/m2
 
-    f_specific_mass = interp1d(temp_x, specific_mass, fill_value = "extrapolate")
-    specific_mass_temp = f_specific_mass(temp)
+    specific_mass = -0.00342727*temp + 1.275 
+    thermal_conductivity = 7.32727273e-05*temp + 0.02428182
+    dynamic_viscosity = 4.7e-08*temp + 1.75045455e-05
 
-    f_thermal_conductivity = interp1d(temp_x, thermal_conductivity, fill_value = "extrapolate")
-    thermal_conductivity_temp = f_thermal_conductivity(temp)
-
-    f_dynamic_viscosity = interp1d(temp_x, dynamic_viscosity, fill_value = "extrapolate")
-    dynamic_viscosity_temp = f_dynamic_viscosity(temp)
-
-    return specific_mass_temp, thermal_conductivity_temp, dynamic_viscosity_temp
+    return specific_mass, thermal_conductivity, dynamic_viscosity
 
 hold_values = {}
 def hold0(name,value=None):
@@ -221,8 +212,8 @@ def get_line_tension_coefficient(d_hori):
     """
     # reference: IEC 60826:2017
     # fit of data is used
-    # 0.2106718346 + 0.0003126614987*d_hori
-    return 0.50723325 - 0.29656325/(0.05*d_hori+1)
+    # 0.50723325 - 0.29656325/(0.05*d_hori+1)
+    return 0.2106718346 + 0.0003126614987*d_hori
 
 def get_distance(pos1, pos2):
     """Compute haversine distance between two locations
@@ -366,28 +357,29 @@ def linesag(data):
                 p0 = [line['latitude'],line['longitude']]
                 d0 = 0.0
                 z0 = elevation + line['pole_height']
-                sag = z0 - elevation
+                sag = line['pole_height']
                 result[id] = sag
 
             else: # end of a line segment
 
                 p1 = [line['latitude'],line['longitude']]
                 d1 = get_distance(p0,p1)
                 z1 = elevation + line['pole_height']
-                sag = z1 - elevation
+                sag = line['pole_height']
                 result[id] = sag
 
                 # compute linesag at waypoints
                 for n,l in ld.items():
                     p = [l['latitude'],l['longitude']]
-                    d_hori = get_distance(p0,p)
+                    d_hori = get_distance(p0,p1)
                     try:
                         elevation = l['elevation']
                     except:
                         pass
                     sag = get_sag_value(d_hori,line,cable,p0,p,z0,z1,
                         power_flow,global_horizontal_irradiance,ground_reflectance,
                         ice_thickness,wind_direction,air_temperature,wind_speed,ice_density)
+                    # the line sag is defined as the distance between the interested point on the line and the ground
                     sag = sag - elevation
                     result[n] = round(sag,OPTIONS["precision"]["linesag"])
 
@@ -423,52 +415,54 @@ def get_sag_value(d_hori,line,cable,p0,p1,z0,z1,
     span = sqrt(d_hori*d_hori + d_vert*d_vert)
     rts = cable['rated_tensile_strength']
     unit_weight = cable['unit_weight']
-    # if d_vert > 30:
-    #     k_init = min(k_init, 3.0*unit_weight*span*span/(2*d_vert*rts))
+    diameter = cable['diameter']
+    air_mass, k_f, air_viscosity = get_air_properties(air_temperature)
+    # calculate the line heading and the angle between line heading and wind direction
+    try:
+        line_angle = line['heading']
+    except:
+        line_angle = 180*atan2(p1[0]-p0[0],p1[1]-p0[1])/np.pi
+    phi = (wind_direction - line_angle)*np.pi/180
+    # calculate the new line sag at loaded condition
+    ice_unit_weight = ice_density*np.pi*ice_thickness*(diameter+ice_thickness)*g
+    wind_unit_weight = 0.5*air_mass*(wind_speed*sin(phi))**2 *(diameter+2*ice_thickness)
+    total_unit_weight = sqrt(wind_unit_weight**2+(unit_weight+ice_unit_weight)**2)
+    if d_vert/d_hori > 0.1:
+        k_init = min(k_init, total_unit_weight*d_hori*d_hori/(2*d_vert*rts))
     H_init = rts*k_init
     sag_init = unit_weight*span*span/(8*H_init)
+    # for Q_I
     P_rated = power_flow
     Vll_rated = cable['voltage_rating']
     Irms =  P_rated/(sqrt(3)*Vll_rated)
     R_20C = cable['nominal_resistance']
     coeff_Al = cable['resistivity']
     Q_I_coeff_first = Irms*Irms*R_20C*coeff_Al
     Q_I_coeff_constant = Irms*Irms*R_20C*(1-coeff_Al*(20.0+273.0))
+    # for Q_S
     k_a = 1.0 - cable['reflectivity'] # solar radiation absorption coefficient
     GHI = global_horizontal_irradiance # unit: W/m2
     k_g = ground_reflectance # ground reflect
-    diameter = cable['diameter']
     Q_S_constant = k_a*(diameter+2*ice_thickness)*(1+k_g)*GHI
-
-    try:
-        line_angle = line['heading']
-    except:
-        line_angle = 180*atan2(p1[0]-p0[0],p1[1]-p0[1])/np.pi
-
-    phi = (wind_direction - line_angle)*np.pi/180
+    # for Q_C
     k_angle = 1.194 - cos(phi) + 0.194*cos(2*phi) + 0.368*sin(2*phi)
-    air_mass, k_f, air_viscosity = get_air_properties(air_temperature)
     Nre = wind_speed * (diameter + 2*ice_thickness) * air_mass / air_viscosity
-    if wind_speed > 1.1:
+    if wind_speed < 0.82:
         Q_C_coeff_first = k_angle*(1.01+1.35*Nre**0.52)*k_f
         Q_C_constant = -k_angle*(1.01+1.35*Nre**0.52)*k_f*(air_temperature+273.0)
     else:
         Q_C_coeff_first = k_angle*0.754*(Nre**0.6)*k_f
         Q_C_constant = -k_angle*0.754*(Nre**0.6)*k_f*(air_temperature+273.0)
     # for Q_R
     k_e = cable['emissivity']
-    Q_R_constant = -5.6704e-8*k_e*(diameter+2.0*ice_thickness)*(air_temperature+273.0)**4
-    Q_R_coeff_fourth = 5.6704e-8*k_e*(diameter+2.0*ice_thickness)
+    Q_R_constant = -5.6704e-8*k_e*(diameter+2.0*ice_thickness)*np.pi*(air_temperature+273.0)**4
+    Q_R_coeff_fourth = 5.6704e-8*k_e*(diameter+2.0*ice_thickness)*np.pi
     # for new conductor temp under loading
     coef_sag = [-Q_R_coeff_fourth,0.0,0.0,Q_I_coeff_first-Q_C_coeff_first,Q_I_coeff_constant+Q_S_constant-Q_C_constant-Q_R_constant]
     r = np.roots(coef_sag)
     r = r[~np.iscomplex(r)]
     temp_load = np.absolute(r[r > 0.0]) # unit: K
     temp_load = temp_load - 273.0 # unit: DegC
-    # calculate the new line sag at loaded condition
-    ice_unit_weight = ice_density*np.pi*ice_thickness*(diameter+ice_thickness)*g
-    wind_unit_weight = 0.5*air_mass*(wind_speed*sin(phi))**2 *(diameter+2*ice_thickness)
-    total_unit_weight = sqrt(wind_unit_weight**2+(unit_weight+ice_unit_weight)**2)
 
     area = cable['conductor_crosssection_area']
     elasticity = cable['elasticity']
@@ -480,18 +474,20 @@ def get_sag_value(d_hori,line,cable,p0,p1,z0,z1,
     r = np.roots(coef_H)
     r = r[~np.iscomplex(r)]
     H_load = np.absolute(r[r > 0.0])
-    #
+
     sag_load = total_unit_weight*span*span/(8*H_load)
     sag_angle = atan(wind_unit_weight/(ice_unit_weight+unit_weight))
     C_catenary = H_load / (ice_unit_weight+unit_weight)
     if z0 > z1:
-        d0_hori = d_hori*(1+d_vert/(4*sag_load))/2
+        d0_hori = d_hori/2 + H_load*d_vert/(total_unit_weight*d_hori)
+        d1_hori = d_hori - d0_hori
         sag0 = total_unit_weight*d0_hori**2 /(2*H_load)
         dt = get_distance(p0,p1)
         sag0_cosh = sag0 - C_catenary*(np.cosh((dt-d0_hori)/C_catenary)-1)
         sag_elevation = z0 - sag0_cosh*cos(sag_angle)
     else:
-        d0_hori = d_hori*(1-d_vert/(4*sag_load))/2
+        d0_hori = d_hori/2 - H_load*d_vert/(total_unit_weight*d_hori)
+        d1_hori = d_hori - d0_hori
         sag0 = total_unit_weight*d0_hori**2 /(2*H_load)
         dt = get_distance(p0,p1)
         sag0_cosh = sag0 - C_catenary*(np.cosh((dt-d0_hori)/C_catenary)-1)