Fix rare error with distribution of quantity index sets.

dev_notes/todo.txt

@@ -2,6 +2,12 @@
 
 	*** High pri ***
 
+		Index set dependency propagation along connections is wrong in some edge cases.
+
+			e.g if it transports a quantity that is indexed over an index set, this index set dependency is not propagated by the connection aggregation.
+
+
+
 		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.
@@ -168,15 +174,6 @@
 
 	*** Intermediate-pri ***
 
-		Got an (internal) error in the following case (nivafjord pipe boundary model)
-			flux(bnd.water, layer.water[vert.specific]) ...
-
-			bnd.water.tox can be distributed over tox_index, but that dependency was ruled out due to lack of differing parametrization
-			layer.water.tox was distributed over tox_index with the actual dependency.
-			Caused problem when the dissolved flux were to be added to the connection aggregate (lack of index set).
-				Probably also a problem for vert.top etc.
-				It is a bit tricky for the framework to predetermine this.. But there is maybe something like this for graph aggregates already that could be applied here too.
-
 		If an input series is provided indexed over "subcatchment", but should be indexed over "water body" and "subcatchment" is a union member of "water body", that should be made to work. (and all similar cases) Just re-map the index.
 
 		Store name of expected model file in data sets (?)

models/demonstration_models/nivafjord_lake_simplycnp_glacier.txt

@@ -27,4 +27,24 @@ model("NIVAFjord-SimplyCNP-Glacier") {
 
 	# Add a module for the lake to discharge out of the system.
 	load("modules/nivafjord/lake_discharge.txt", module("NIVAFjord lake discharge", basin, layer, water, dz, vert, dims))
+
+	ti : index_set("Tracer index")
+	tr : quantity("Tracer", ti)
+
+	module("Extra things", version(0, 0, 0)) {
+
+		par_group("Tracer", tr) {
+			init_tr : par_real("Initial tracer", [m g, l-1], 1)
+		}
+
+		var(soil.water.tr, [k g, k m-2], [m g, l-1], "Soil water tracer") @initial_conc { init_tr }
+
+		var(gw.water.tr, [k g, k m-2], [m g, l-1], "Groundwater tracer")
+
+		var(river.water.tr, [k g], [m g, l-1], "River water tracer")
+
+		var(layer.water.tr, [k g], [m g, l-1], "Layer water tracer")
+
+
+	}
 }

models/example_data/NIVAFjord/glacier_lake_example.dat

@@ -19,6 +19,8 @@ data_set {
 
 	phyt_type : index_set("Phytoplankton type") [ "Phyto" ]
 
+	ti : index_set("Tracer index") [ "Cl" "Ca" ]
+
 	position_map(layer) [!
 		0 : 1
 	]

src/model_compilation.cpp

@@ -18,7 +18,10 @@ std::string
 Model_Instruction::debug_string(Model_Application *app) const {
 	std::stringstream ss;
 
-	if(type == Model_Instruction::Type::compute_state_var)
+	if(type == Model_Instruction::Type::invalid)
+		// NOTE: This sometimes happens when it prints instruction dependencies.
+		ss << "(invalid)";
+	else if(type == Model_Instruction::Type::compute_state_var)
 		ss << "\"" << app->vars[var_id]->name << "\"";
 	else if(type == Model_Instruction::Type::subtract_discrete_flux_from_source)
 		ss << "\"" << app->vars[source_id]->name << "\" -= \"" << app->vars[var_id]->name << "\"";
@@ -336,6 +339,25 @@ insert_var_order_depencencies(Model_Application *app, Model_Instruction *instr,
 	}
 }
 
+Index_Set_Tuple
+get_quantity_index_sets(Model_Application *app, Var_Id var_id) {
+	Index_Set_Tuple result;
+
+	auto var = app->vars[var_id];
+	auto &loc = var->loc1;
+	if(!is_located(loc))
+		fatal_error(Mobius_Error::internal, "Misuse of get_quantity_index_sets");
+	// Iterate from 1 to not include the compartment.
+	for(int idx = 1; idx < loc.n_components; ++idx) {
+		auto quantity = app->model->components[loc.components[idx]];
+		if(quantity->decl_type != Decl_Type::quantity)
+			fatal_error(Mobius_Error::internal, "Misuse of get_quantity_index_sets");
+		for(auto index_set : quantity->index_sets)
+			result.insert(index_set);
+	}
+	return result;
+}
+
 void
 resolve_basic_dependencies(Model_Application *app, std::vector<Model_Instruction> &instructions) {
 
@@ -351,6 +373,13 @@ resolve_basic_dependencies(Model_Application *app, std::vector<Model_Instruction
 			auto var = app->vars[instr.var_id];
 			if(var->specific_target.get())
 				register_dependencies(var->specific_target.get(), &code_depends);
+			if(var->is_mass_balance_quantity()) {
+				// By default always distribute over quantity distributions.
+				// Otherwise there are several failure edge cases where these are not propagated (through connections mostly).
+				// And there is little to gain by optimizing them away
+				auto quantity_index_sets = get_quantity_index_sets(app, instr.var_id);
+				insert_dependecies(app, &instr, quantity_index_sets);
+			}
 		}
 
 		for(auto &dep : code_depends) {
@@ -820,22 +849,21 @@ process_graph_connection_aggregation(Model_Application *app, std::vector<Model_I
 		insert_dependency(app, &instructions[flux_id.id], conn->edge_index_set);
 	}
 
-	if(!agg_var->is_out) { // TODO: should (something like) this also be done for the source aggregate in directed_graph?
-
-		// TODO: Make a better explanation of what is going on in this block and why it is needed (What is the failure case otherwise).
+	if(!agg_var->is_out) {
 
 		// If the target compartment (not just what the connection indexes over) has an index set shared with the source compartment, we must index the target variable over that.
+		// Since the target could get a different value from the connection depending on its own index, we have to force it to be computed per each of these indexes even if it were not to have an index set dependency on this otherwise.
+
 		auto target_comp_id = app->vars[agg_var->agg_for]->loc1.first();
 		auto find_target = app->find_connection_component(agg_var->connection, target_comp_id);
 
 		auto target_comp = model->components[target_comp_id];
-		auto target_index_sets = find_target->index_sets; // vector copy;
+		auto target_index_sets = find_target->index_sets; // intentional copy;
 		for(auto index_set : find_source->index_sets) {
 			if(std::find(target_comp->index_sets.begin(), target_comp->index_sets.end(), index_set) != target_comp->index_sets.end())
 				target_index_sets.push_back(index_set);
 		}
 
-		// Since the target could get a different value from the connection depending on its own index, we have to force it to be computed per each of these indexes even if it were not to have an index set dependency on this otherwise.
 		for(auto index_set : target_index_sets)
 			insert_dependency(app, &instructions[agg_var->agg_for.id], index_set);
 	}

test/models/another_simple.txt

@@ -2,21 +2,34 @@
 
 model("Another simple") {
 
+	i : index_set("Q index")
+
 	a : compartment("A")
-	q : quantity("Q")
+	b : compartment("B")
+	r : quantity("R")
+	q : quantity("Q", i)
+
+	c : connection("Connection") @directed_graph { a b }
 
 	module("So simple", version(2, 0, 0)) {
 
-		var(a.q, [], "A Queue!") @override {
-			time.step=>>
-		} @no_store
+		par_group("Q pars", q) {
+			init_q : par_real("Init Q", [], 1)
+		}
+
+		var(a.r, [k g], "A water") @initial { 10[k g] }
+
+		var(a.r.q, [k g], "A queue") @initial_conc { init_q }
+
+		var(b.r, [k g], "B water")
+
+		var(b.r.q, [k g], "B queue")
 
-		r : property("P")
-		var(a.r, [], "Arrr!") { a.q*2 }
+		flux(a.r, c, [k g, day-1], "A r flux") { a.r*0.01[day-1] }
 
 	}
 
 	sol : solver("solver", inca_dascru, [2, hr], 0.01)
 
-	solve(sol, a.q)
+	solve(sol, a.r, b.r)
 }

test/models/another_simple_data.dat

@@ -4,6 +4,18 @@ data_set {
 The doc string
 """
 
+	i : index_set("Q index") [ "Q1" "Q2" ]
+
+	connection("Connection") {
+		a : compartment("A")
+		b : compartment("B")
+
+		directed_graph [
+			a[] -> b[]
+		]
+	}
+
+
 	par_group("System") {
 
 		par_datetime("Start date")