More work on karstic groundwater, allow logging in mobipy

dev_notes/todo.txt

@@ -2,6 +2,15 @@
 
 	*** High pri ***
 
+		Make functionality to patch in data from a sub-model in a data set.
+			E.g. you have a setup for Simply+NIVAFjord and you want to patch in from a different file that just has the Simply part of the setup.
+			Should not be an "include" on file level, because that causes too much ambiguity about how this should be saved if overwritten.
+			Instead, make it an explicit operation accessible from API (thus MobiView2 and/or mobipy).
+			void patch(Data_Set *target, Data_Set *source); // Or something like that.
+
+			This is useful because often you calibrate the catchment separately (runs much faster alone), then want to combine the new results with the full basin model.
+				Is currently done with text copy pasting, but can be a bit tiresome.
+
 		Store name of expected model file in data sets (?)
 			Then allow loading that file only, and it infers the model it needs.
 				Would probably only work if it is in /Mobius2/models/

mobipy/__init__.py

@@ -131,7 +131,7 @@ def load_dll() :
 	dll.mobius_get_steps.argtypes = [ctypes.c_void_p, ctypes.c_int32]
 	dll.mobius_get_steps.restype = ctypes.c_int64
 
-	dll.mobius_run_model.argtypes = [ctypes.c_void_p, ctypes.c_int64]
+	dll.mobius_run_model.argtypes = [ctypes.c_void_p, ctypes.c_int64, ctypes.CFUNCTYPE(None, ctypes.c_void_p, ctypes.c_double)]
 	dll.mobius_run_model.restype = ctypes.c_bool
 
 	dll.mobius_get_time_step_size.argtypes = [ctypes.c_void_p]
@@ -201,14 +201,14 @@ def load_dll() :
 
 	dll.mobius_index_names.argtypes = [ctypes.c_void_p, Entity_Id, ctypes.POINTER(Mobius_Index_Value)]
 
-	#dll.mobius_allow_logging.argtypes = [ctypes.c_bool]
+	dll.mobius_allow_logging.argtypes = [ctypes.c_bool]
 
 	return dll
 
 dll = load_dll()
 
-#def toggle_logging(allow) :
-#	dll.mobius_allow_logging(allow)
+def allow_logging(allow=True) :
+	dll.mobius_allow_logging(allow)
 
 def _check_for_errors() :
 	buflen = 1024
@@ -390,6 +390,11 @@ def list_all(self, type) :
 		return [(id.decode('utf-8'), n.decode('utf-8')) for id, n in zip(idents, names)]
 
 
+# TODO: Maybe make something nicer here.
+@ctypes.CFUNCTYPE(None, ctypes.c_void_p, ctypes.c_double)
+def _run_logger(_p, percent) :
+	print("Run progress: %g%%"%percent)
+
 class Model_Application(Scope) :
 	def __init__(self, data_ptr, is_main) :
 		super().__init__(data_ptr, invalid_entity_id)
@@ -431,8 +436,12 @@ def copy(self, copy_results = False) :
 		new_ptr = dll.mobius_copy_data(self.data_ptr, copy_results)
 		return Model_Application(new_ptr, False)
 
-	def run(self, ms_timeout=-1) :
-		finished = dll.mobius_run_model(self.data_ptr, ms_timeout)
+	def run(self, ms_timeout=-1, log=False) :
+		if log :
+			finished = dll.mobius_run_model(self.data_ptr, ms_timeout, _run_logger)
+		else :
+			no_cb = ctypes.CFUNCTYPE(None, ctypes.c_void_p, ctypes.c_double)(0)
+			finished = dll.mobius_run_model(self.data_ptr, ms_timeout, no_cb)
 		_check_for_errors()
 		return finished
 

mobius_jl/mobius.jl

@@ -134,9 +134,10 @@ end
 
 finalize!(data::Model_Data) = free(data)
 
+# TODO: Could allow callback here too..
 function run_model(data::Model_Data, ms_timeout::Int64=-1)::Bool
-	result = ccall(run_model_h, Cint, (Ptr{Cvoid}, Clonglong),
-		data.ptr, ms_timeout)
+	result = ccall(run_model_h, Cint, (Ptr{Cvoid}, Clonglong, Ptr{Cvoid}),
+		data.ptr, ms_timeout, C_NULL)
 	check_error()
 	return result
 end

models/modules/simplyq.txt

@@ -105,6 +105,7 @@ module("SimplyQ groundwater karstic", version(0, 0, 0),
 		tc_g    : par_real("Groundwater time constant", [day],          30,  1,  400)
 		gw_ret  : par_real("Groundwater retention volume", [m m],        0,  1, 1e4)
 		gw_stf  : par_real("Groundwater subterranean flow fraction", [], 0, 0, 1)
+		gw_init : par_real("Initial groundwater excess volume", [m m],   0,  0, 1e4)
 	}
 
 	par_group("Groundwater extraction", gw) {
@@ -115,7 +116,7 @@ module("SimplyQ groundwater karstic", version(0, 0, 0),
 
 	flux(gw.water, out, [m m, day-1], "Groundwater extraction rate") { gw_ext*s_response(gw.water, 0, 10[m m], 0, 1) }
 
-	var(gw.water,   [m m], "Groundwater volume") @initial { gw_ret }
+	var(gw.water,   [m m], "Groundwater volume") @initial { gw_ret + gw_init }
 
 	var(gw.water.flow, [m m, day-1], "Groundwater flow") { max(0, water-gw_ret) / tc_g }
 

models/nivafjord_simplycnp_model.txt

@@ -28,7 +28,8 @@ model("NIVAFjord-SimplyCNP") {
 		a_vap : property,
 		s_vap : property,
 		rho : property,
-		sol : solver
+		sol : solver,
+		gw_st : connection,
 	)
 
 	air        : compartment("Atmosphere", wb_index)
@@ -54,5 +55,33 @@ model("NIVAFjord-SimplyCNP") {
 	horz : connection("Fjord horizontal") @directed_graph(e_idx) { 
 		(river? layer+ bnd_layer?) | river
 	} @no_cycles
+
+	option(gw_module.karstic) {
+
+		# Allow connecting subterranean groundwater flow directly to ocean/lake basins.
+		gw_st : connection("Groundwater subterranean") @directed_graph {
+			gw+ layer
+		} @no_cycles
+
+		unit_conversion(gw.water, layer.water) { a_catch->> }
+
+		# TODO: Maybe factor out
+		# Needed since we have direct discharge from gw to lake
+		module("Groundwater temp and O2", version(0,0,0)) {
+
+			load("stdlib/physiochemistry.txt", library("Water utils"))
+
+			var(gw.water.heat, [J, k m-2], "Groundwater heat") @override {
+				t := aggregate(soil.temp),
+				V := 1[k m 2]*gw.water -> [m 3],
+				(water_temp_to_heat(V, t) / 1[k m 2]) ->>
+			}
+
+			var(gw.water.o2, [k g, k m-2], "Groundwater oxygen") @override_conc {
+				conc(river.water.o2) # TODO: do something else here?
+			}
+		}
+
+	}
 }
 

src/c_abi.cpp

@@ -139,10 +139,10 @@ mobius_save_data_set(Model_Data *data, char *data_file) {
 }
 
 DLLEXPORT bool
-mobius_run_model(Model_Data *data, s64 ms_timeout) {
+mobius_run_model(Model_Data *data, s64 ms_timeout, run_callback_type run_callback) {
 
 	try {
-		return run_model(data, ms_timeout);
+		return run_model(data, ms_timeout, false, run_callback);
 	} catch(int) {}
 	return false;
 }

src/c_abi.h

@@ -80,7 +80,7 @@ DLLEXPORT void
 mobius_save_data_set(Model_Data *data, char *data_file);
 
 DLLEXPORT bool
-mobius_run_model(Model_Data *data, s64 ms_timeout);
+mobius_run_model(Model_Data *data, s64 ms_timeout, run_callback_type run_callback);
 
 DLLEXPORT s64
 mobius_get_steps(Model_Data *data, Var_Id::Type type);

src/data_set.cpp

@@ -1297,7 +1297,7 @@ write_component_ref(Data_Set *data_set, Scope_Writer *writer, Compartment_Ref &r
 			auto index_set = data_set->index_sets[index_set_id];
 
 			bool quote;
-			std::string index_name = data_set->index_data.get_index_name(ref.indexes, ref.indexes.indexes[idx], &quote);
+			std::string index_name = data_set->index_data.get_index_name(ref.indexes, ref.indexes.indexes[idx], &quote, true);
 			maybe_quote(index_name, quote);
 
 			writer->write(" %s", index_name.data());

src/index_data.cpp

@@ -596,7 +596,7 @@ Index_Data::get_instance_count(const std::vector<Entity_Id> &index_sets) {
 }
 
 std::string
-Index_Data::get_index_name_base(Index_T index, Index_T index_of_super, bool *is_quotable) {
+Index_Data::get_index_name_base(Index_T index, Index_T index_of_super, bool *is_quotable, bool numeric_plain) {
 	auto &data = index_data[index.index_set.id];
 
 	// TODO: Remove this once we fix MobiView2
@@ -611,8 +611,10 @@ Index_Data::get_index_name_base(Index_T index, Index_T index_of_super, bool *is_
 			if(index.index > 0)
 				from = data.pos_vals[index.index-1];
 			double to = data.pos_vals[index.index];
-			sprintf(buf, "%.15g-%.15g", from, to);
-			//sprintf(buf, ".15g", from);
+			if(numeric_plain)
+				sprintf(buf, "%.15g", from);
+			else
+				sprintf(buf, "%.15g-%.15g", from, to);
 			return buf;
 		} else
 			return std::to_string(index.index);
@@ -630,7 +632,7 @@ Index_Data::get_index_name_base(Index_T index, Index_T index_of_super, bool *is_
 }
 
 std::string
-Index_Data::get_index_name(Indexes &indexes, Index_T index, bool *is_quotable) {
+Index_Data::get_index_name(Indexes &indexes, Index_T index, bool *is_quotable, bool numeric_plain) {
 
 	auto index_set = catalog->index_sets[index.index_set];
 	Index_T index_of_super = Index_T::no_index();
@@ -649,10 +651,10 @@ Index_Data::get_index_name(Indexes &indexes, Index_T index, bool *is_quotable) {
 
 	if(!index_set->union_of.empty() && data.type == Index_Record::Type::named) {
 		Index_T below = lower(index, index_of_super);
-		return get_index_name_base(below, index_of_super, is_quotable);
+		return get_index_name_base(below, index_of_super, is_quotable, numeric_plain);
 	}
 
-	return get_index_name_base(index, index_of_super, is_quotable);
+	return get_index_name_base(index, index_of_super, is_quotable, numeric_plain);
 }
 
 std::string

src/index_data.h

@@ -107,7 +107,7 @@ Index_Data {
 	void check_valid_distribution(std::vector<Entity_Id> &index_sets, Source_Location source_loc);
 	s64 get_instance_count(const std::vector<Entity_Id> &index_sets);
 
-	std::string get_index_name(Indexes &indexes, Index_T index, bool *is_quotable = nullptr);
+	std::string get_index_name(Indexes &indexes, Index_T index, bool *is_quotable = nullptr, bool numeric_plain = false);
 	std::string get_possibly_quoted_index_name(Indexes &indexes, Index_T index, bool quote = true);
 	void get_index_names(Indexes &indexes, std::vector<std::string> &names_out, bool quote = false);
 
@@ -143,7 +143,7 @@ private :
 
 	Index_T find_index_base(Entity_Id index_set, Token *idx_name, Index_T index_of_super = Index_T::no_index());
 	s32   get_count_base(Entity_Id index_set, Index_T index_of_super = Index_T::no_index());
-	std::string get_index_name_base(Index_T index, Index_T index_of_super, bool *is_quotable);
+	std::string get_index_name_base(Index_T index, Index_T index_of_super, bool *is_quotable, bool numeric_plain=false);
 
 	void initialize(Entity_Id index_set_id, Index_T parent_idx, Index_Record::Type type, Source_Location source_loc);
 

src/model_application.cpp

@@ -956,7 +956,8 @@ pre_process_connection_data(Model_Application *app, Data_Set *data_set, Entity_I
 		// NOTE: We still have to allow it to be edge indexed even if max outgoing <= 1 since that is dependent on the data, not on the model.
 	}
 
-	match_regex(app, conn_id, connection_data->source_loc);
+	// TODO: This is broken, needs to be fixed.
+	//match_regex(app, conn_id, connection_data->source_loc);
 }
 
 void

src/model_composition.cpp

@@ -1101,7 +1101,13 @@ get_unit_conversion(Model_Application *app, Var_Location &loc1, Var_Location &lo
 
 	if(!found_conv && need_conv) {
 		// TODO: Better error if it is a connection_aggregate, not a flux.
-		fatal_error(Mobius_Error::model_building, "The units of the source and target of the flux \"", app->vars[flux_id]->name, "\" are not the same, but no unit conversion are provided between them in the model.");
+		begin_error(Mobius_Error::model_building);
+		error_print("A unit conversion is needed between the two following locations since there is a flux between them:\n");
+		error_print_location(model, loc1);
+		error_print("\n");
+		error_print_location(model, loc2);
+		mobius_error_exit();
+		//fatal_error(Mobius_Error::model_building, "The units of the source and target of the flux \"", app->vars[flux_id]->name, "\" are not the same, but no unit conversion are provided between them in the model.");
 	}
 
 	auto ast   = found_conv->code;
@@ -1138,6 +1144,7 @@ add_connection_agg_weights(Model_Application *app, std::vector<Conversion_Data>
 	Conversion_Data data;
 	data.source_id = app->vars.id_of(loc1);
 	data.weight = owns_code(get_aggregation_weight(app, loc1, target_comp, conn_id));
+	// TODO: It passes the target_var_id, but the function expects a flux_id. This only matters for the error message that is printed though.
 	data.unit_conv = owns_code(get_unit_conversion(app, loc1, loc2, target_var_id));
 	if(data.weight || data.unit_conv) conv_data.push_back(std::move(data));
 }