Some fixes

mobipy/__init__.py

@@ -70,12 +70,25 @@ class Mobius_Entity_Metadata(ctypes.Structure) :
 	]
 
 
-# NOTE: Must match Reg_Type: We should find a way to auto-generate this from reg_types.incl instead!
-MODULE_TYPE = 1
-COMPONENT_TYPE = 2
-PARAMETER_TYPE = 3
-FLUX_TYPE = 4
-PAR_GROUP_TYPE = 6
+# The below recreates the Reg_Type enum using global variables, equivalent to
+# It is auto-generated so that it is guaranteed to match the C++ code.
+# MODULE_TYPE = 1
+# COMPONENT_TYPE = 2
+# PARAMETER_TYPE = 3
+# FLUX_TYPE = 4
+# PAR_GROUP_TYPE = 6
+# ...
+
+with open('../src/reg_types.incl', 'r') as f :
+	t = f.read()
+	idx = 1
+	for line in t.split('\n') :
+		if line.startswith('//') : continue
+		glob_name = '%s_TYPE' % line[line.find('(')+1:line.rfind(')')].upper()
+		globals()[glob_name] = idx
+		idx += 1
+
+
 
 def load_dll() :
 	shared_ext = 'dll'
@@ -431,10 +444,10 @@ def name(self) :
 
 	def min(self) :
 		if self.entity_id.reg_type != PARAMETER_TYPE :
-			raise ValueError("This entity doesn't have a min value.")
+			raise ValueError("This entity type doesn't have a min value.")
 		type = dll.mobius_get_value_type(self.data_ptr, entity_id)
 		if type > 1 :
-			raise ValueError("This parameter does not have a min value.")
+			raise ValueError("This parameter type does not have a min value.")
 		data = dll.mobius_get_entity_metadata(self.data_ptr, self.entity_id)
 		_check_for_errors()
 		if type == 0 :
@@ -444,10 +457,10 @@ def min(self) :
 
 	def max(self) :
 		if self.entity_id.reg_type != PARAMETER_TYPE :
-			raise ValueError("This entity doesn't have a max value.")
+			raise ValueError("This entity type doesn't have a max value.")
 		type = dll.mobius_get_value_type(self.data_ptr, self.entity_id)
 		if type > 1 :
-			raise ValueError("This parameter does not have a max value.")
+			raise ValueError("This parameter type does not have a max value.")
 		data = dll.mobius_get_entity_metadata(self.data_ptr, self.entity_id)
 		_check_for_errors()
 		if type == 0 :
@@ -459,14 +472,14 @@ def max(self) :
 
 	def description(self) :
 		if self.entity_id.reg_type != PARAMETER_TYPE :
-			raise ValueError("This entity doesn't have a description.")
+			raise ValueError("This entity type doesn't have a description.")
 		data = dll.mobius_get_entity_metadata(self.data_ptr, self.entity_id)
 		_check_for_errors()
 		return data.description.decode('utf-8')
 
 	def unit(self) :
 		if self.entity_id.reg_type != PARAMETER_TYPE :
-			raise ValueError("This entity doesn't have a unit.")
+			raise ValueError("This entity type doesn't have a unit.")
 		data = dll.mobius_get_entity_metadata(self.data_ptr, self.entity_id)
 		_check_for_errors()
 		return data.unit.decode('utf-8')
@@ -552,7 +565,7 @@ def __setitem__(self, indexes, values) :
 			raise ValueError("Slices not yet supported for setting input series")
 
 		if not isinstance(values, pd.Series) :
-			raise ValueError("Expected a pandas Series object for the values")
+			raise ValueError("Expected a pandas.Series object for the values")
 
 		time_steps = len(values)
 		dates = (ctypes.c_int64 * time_steps)(*(ts.astype('datetime64[s]').astype(np.int64) for ts in values.index.values))

mobipy/optim.py

@@ -10,7 +10,7 @@
 # Note: we can't use multiprocessing for this, only multithreading, since a Model_Application object that is allocated from C++ on one process
 # can't be accessed from a different python process.
 
-# Oooops, this doesn't work as intended yet.
+# Oooops, this doesn't work as intended. Why??
 #class Thread_Pool :
 
 #    def map(self, fn, args) :

models/nivafjord_simplycnp_isolated_basin_model.txt

@@ -2,16 +2,6 @@
 
 model("NIVAFjord-SimplyCNP-IsolatedBasin") {
 
-	extend("simplycnp_model.txt") @exclude(
-		downstream : connection, 
-		air : compartment, 
-		water : quantity,
-		heat : quantity,
-		temp : property,
-		precip : property,
-		sol : solver
-	)
-
 	basin_idx     : index_set("Basin index")
 	layer_idx     : index_set("Layer index") @sub(basin_idx)
 
@@ -25,6 +15,13 @@ model("NIVAFjord-SimplyCNP-IsolatedBasin") {
 	salt  : quantity("Salt")
 	heat  : quantity("Heat energy")
 
+	sed : quantity("Sediments")
+	oc : quantity("Organic carbon")
+	on : quantity("Organic nitrogen")
+	op : quantity("Organic phosphorous")
+	din : quantity("Nitrate")
+	phos : quantity("Phosphate")
+
 	temp     : property("Temperature")
 	salinity : property("Salinity")
 	precip   : property("Precipitation")
@@ -59,6 +56,8 @@ model("NIVAFjord-SimplyCNP-IsolatedBasin") {
 	load("modules/atmospheric.txt",
 		module("Atmospheric", air, temp, wind, g_rad, pressure, a_hum, rho, lwd, cos_z, a_vap, s_vap))
 
+
+	#TODO: Replace with the one that has FPV covering.
 	load("modules/airsea.txt",
 		module("AirSea", "AirSea Fjord", air, basin, ice, heat, temp, precip, wind, g_rad, pressure, rho, a_hum, lwd, sw, attn, indicator,
 							evap, cos_z, loc(basin.freeze_t), loc(basin.area), loc(layer.water[vert.top]), loc(layer.water.heat[sw_vert.top])))
@@ -67,8 +66,6 @@ model("NIVAFjord-SimplyCNP-IsolatedBasin") {
 		dims : preamble("NIVAFjord dimensions", basin_idx, layer_idx, layer),
 		module("NIVAFjord basin", air, basin, layer, water, salt, heat, temp, salinity, pressure, wind, g_rad, rho, attn, z, dz, h, area, freeze_t, sw, vert, sw_vert, loc(sediment.heat), dims))
 
-	# TODO: Need a module that positions the horizontal flux from the river.
-
 	sol : solver("NIVAFjord solver", inca_dascru, solver_h, solver_hmin)
 	solve(sol, layer.water, basin.ice)
 
@@ -84,29 +81,48 @@ model("NIVAFjord-SimplyCNP-IsolatedBasin") {
 	load("modules/nivafjord/sediment.txt",
 		module("NIVAFjord sediments", layer, sediment, water, o2, sed, oc, din, on, phos, op, heat, area, temp, dz, vert, chem_par, dims))
 
-	load("modules/easychem.txt",
-		module("Simple River O₂", river, water, o2, temp))
-
 	solve(sol, sediment.sed, sediment.heat)
 	solve(sol, air.cos_z)
 
-	wb_index : index_set("Water body") @union(sc, basin_idx)
-
-	downstream : connection("Downstream") @directed_graph {
-		river+
-	} @no_cycles
-
-	horz : connection("Fjord horizontal") @directed_graph { 
-		(river? layer+) | river
-	} @no_cycles
-
-
-	module("River-Basin simple flux", version(0, 0, 0)) {
+	module("Basin inflow", version(0, 0, 0)) {
+
+		var(air.precip, [m m, day-1], "Precipitation")
+		var(air.temp, [deg_c], "Air temperature")
+
+
+		load("stdlib/seawater.txt", library("Sea oxygen"))
+		load("stdlib/physiochemistry.txt", library("Chemistry"), library("Water utils"))
+
+		# Water
+
+		flow : property("Basin inflow") # Input series
+		var(basin.flow, [m 3, s-1])
 
 		# Just directing the catchment runoff to the top layer instead of doing the density check.
-		flux(river.water, layer.water[horz.below, vert.top], [m 3, s-1], "River discharge to fjord") {
-			flow->>
+		flux(out, layer.water[vert.top], [m 3, s-1], "Discharge from land to basin") { basin.flow->> }
+
+		# O2
+
+		par_group("Inflow oxygen") {
+			f_o2sat    : par_real("Inflow O₂ saturation fraction", [], 0.9, 0, 1)
+		}
+
+		flux(out, layer.water.o2[vert.top], [k g, day-1], "Oxygen from land") {
+			inflow_t := air.temp, # Assume inflow temperature = air temperature (see also below)
+			land_conc := f_o2sat*o2_saturation(inflow_t, 0) * o2_mol_mass,
+			basin.flow*land_conc->>
 		}
+
+		# Heat / temperature
+
+		flux(out, layer.water.heat[vert.top], [J, day-1], "Heat from land") {
+			inflow_t := air.temp, # Assume inflow temperature = air temperature (see also above)
+			(water_temp_to_heat(basin.flow => [m 3], inflow_t) => [J, s-1]) ->>
+		}
+
+		# TODO: Add other substances that may be washed in.
+
+		# TODO: Also outflow?
 	}
 }
 

src/support/residual_structure.cpp

@@ -78,12 +78,6 @@ add_random_error(double* series, s64 time_steps, double *err_param, LL_Type ll_t
 				double std_dev = err_param[0]*series[ts];
 				std::normal_distribution<double> distr(series[ts], std_dev);
 				series[ts] = distr(gen);
-				/*
-				std::normal_distribution<double> distr(0.0, std_dev);
-				double draw = distr(gen);
-				warning_print("Par: ", err_param[0], " series: ", series[ts], " stddev: ", std_dev, " draw: ", draw, "\n");
-				series[ts] = series[ts] + draw;
-				*/
 			} break;
 
 			case LL_Type::wls : {