More improvements

dev_notes/doodle.txt

@@ -1,5 +1,6 @@
 
 
+
 ## How to make index set resizing API (including in mobipy.., where you would most often call it from).
 
 	data = mobipy.Data_Set.from_file('file.dat')

dev_notes/todo.txt

@@ -207,8 +207,7 @@
 				soil.conc_so4 <- something,
 			}
 			This way we could port all of MAGIC to Mobius2 and not call out to a C++ external_computation for it.
-				Although it would be annoying to do the equilibrium constants without array syntax..
-				And we would want multiple return values from functions for this.
+				Although it would be annoying to do the equilibrium constants without array and struct syntax..
 
 		Not possible to scope into a preamble from model scope. Should do this in a good way. Can't always export the preamble to model scope.
 			Maybe we should always scope into preambles using the . syntax? (Also when they are passed to modules?).
@@ -226,9 +225,6 @@
 						@use(a, b2=b, c) ... # (requires new syntax type for the AST, but would also be nice to have for passing load arguments to modules).
 					and also require that when a module takes a preamble as an argument. So that the module can control what the symbols are called in its own scope (and also what symbols it expects to get from what preamble).
 
-		Merge unit_conversion and aggregation_weight to just "weight" (?)
-			Not necessarily a good idea. You may want to look up an aggregate value that is not a flux (and so does not need unit conversion).
-
 		Would be very useful to be able to turn on an option to store out what variables are causing error control to reduce step size.
 		Also, get the average step size through the time step, not just the last one.	
 			Requires a lot of extra infrastructure for this to be possible though.
@@ -237,15 +233,6 @@
 			Turned out to be very difficult due to how things are organized, and caused a totally unrelated series to become NaN.
 				May have messed up indexing of series values?
 
-		Multiple return values from functions (tuples). Needed if we want to implement MAGIC entirely in Mobius2 code for instance.
-			Only annoying thing is that we then have to introduce 'tuple(N)' to the type system and check against it everywhere it is not allowed.
-				A possible solution is to force it to be unpacked immediately at the call location of the function. I.e. if a function returns a tuple, you can only call it using something like
-					a; b := fun(),
-				which de-sugars internally to
-					t := fun()
-					a := unpack(t, 0)
-					b := unpack(t, 1)
-
 		Position_Map
 			There can be problems if the max depth doesn't exactly match a boundary of two widths.
 				See for instance Breiangen in nivafjord_moss
@@ -307,7 +294,7 @@
 
 
 		Data_Set
-			Need a reindex() function on the Data_Set or something like that, which updates all data based on new indexes for the index sets.
+			Need a resize() function on the Data_Set or something like that, which updates all data based on new indexes for the index sets.
 				Must go through everything that is indexed and delete items that are attached to removed indexes.
 					parameters, connections, input series.
 						anything else?
@@ -336,15 +323,14 @@
 			May need a more granular system for specifying what series to store (for non-declared series).
 				Like have a store() declaration in the Data_Set.
 
-			Make @no_store available for @override quantities even on solvers.
-				I.e. for all declared variables that are not is_mass_balance_quantity().
-
 			Determine variables that are going to be run-time constant and just compute them in the initial step, not the rest. (Also goes for stored variables).
-				Not allowed for 'quantity' unless it is @override.
+				Not allowed for ODE quantities.
 					Due to difficulty of determining what it will look like before the code generation step (the way it is set up now).
 				A bit annoying if the @initial code for it is different from the main code..
+					Probably we can't do this in that case.
 
 				Could also be done for variables that are not @no_store in the outset. They could be promoted as no_store.
+					Although that is annoying since then they don't show up in the interface.
 
 
 		Right now you get warning if you load a library that loads another library, but you don't use all the functions of the first library so that the second library is not referenced in active code.

docs/mobius2docs/math_format.md

@@ -212,14 +212,35 @@ These are of the form
 <function-identifier>(<primary-expression>, .., <primary-expression>)
 ```
 
-The function has 0 or more arguments.
+A function takes 0 or more arguments.
 
 The function identifier identifies either a function declaration that is visible in the parent declaration scope, or an intrinsic function.
 
 If it is a declared function, it can have requirements about the units of the arguments, and the result will have the unit of the expression of the body of the function declaration. 
 
 Declared functions are inlined at the site they are evaluated. (This means that you can't have recursive declared functions for now, this may be implemented later).
 
+#### Multiple return values
+
+You can make a function return multiple values by returning a `tuple`. This tuple can then be unpacked into local variables at the call site. For instance, if a function is declared as
+
+```python
+very_fun : function(x) {
+	tuple(x + 1, x + 2)
+}
+```
+
+You can evaluate it in another math block by binding the result to a `;`-separated list of local variable identifiers as follows:
+
+```python
+a;b := very_fun(8),
+# Equivalent to: a := 8+1, b := 8+2,
+```
+
+Any finite number of return values are supported. Note that this is the only intended use of tuples, they are not supported in other contexts at the moment.
+
+#### Intrinsic functions
+
 The following intrinsic functions are visible in every function scope. They are implemented either using [LLVM intrinsics](https://llvm.org/docs/LangRef.html#intrinsic-functions) or [LLVM libc](https://libc.llvm.org/math/index.html).
 
 | Signature  | Description | Units |

models/example_data/NIVAFjord/nivafjord_drammen_different_layer_width.dat

@@ -60,7 +60,7 @@ NOTE: Breiangen is set to be this shallow (12 m) because for boundary basins we
 		[ 2015-01-01 ]
 
 		par_datetime("End date") 
-		[ 2021-12-31 ]
+		[ 2021-11-10 ]
 
 	}
 

models/modules/hbv_snow.txt

@@ -29,8 +29,6 @@ Authors: Magnus D. Norling
 	}
 
 	# Note: We don't give a function body to compute temp and precip, and so if they are not re-declared with a formula in another module, they are assumed to be input series.
-
-	#TODO: The declaration of input series should maybe be in model scope?
 	var(air.temp, [deg_c], "Air temperature")
 	var(air.precip, [m m, day-1])
 
@@ -53,7 +51,6 @@ Authors: Magnus D. Norling
 	flux(snow_box.water, water_target, [m m, day-1], "Melt runoff") {   max(0, water - snow*snow_liq)*1[day-1]  }
 
 	discrete_order {
-		#TODO: we also want to be able to specify special placement of add and sub instructions for these if necessary. But not for this module.
 		p_snow
 		p_rain
 		melt

models/modules/nivafjord/boudary_chem.txt

@@ -32,10 +32,9 @@ module("NIVAFjord boundary chemistry", version(0, 0, 0),
 	dipconc : property("DIP conc")
 	pocconc : property("POC conc")
 
-	# Input series :
+	# In case data was not provided, assume boundary O2 is at saturation
+	# In general, it is better if it is provided as an input series.
 	var(bnd_layer.water.o2conc, [m g, l-1], "Boundary layer O₂ conc") {
-		#In case data was not provided, assume it is at saturation
-		#TODO: This is not a good idea except near surface?
 		o2_saturation_concentration(temp, salinity)*o2_mol_mass->>
 	}
 
@@ -46,49 +45,49 @@ module("NIVAFjord boundary chemistry", version(0, 0, 0),
 	var(bnd_layer.water.pocconc, [m g, l-1], "Boundary layer POC conc")
 
 	var(bnd_layer.water.o2, [k g], [m g, l-1], "Boundary layer O₂")
-		@override_conc { o2conc }
+		@override_conc { o2conc } @no_store
 
 	var(bnd_layer.water.oc, [k g], [m g, l-1], "Boundary layer DOC")
-		@override_conc { docconc }
+		@override_conc { docconc } @no_store
 
 	var(bnd_layer.water.din, [k g], [m g, l-1], "Boundary layer DIN")
-		@override_conc { dinconc }
+		@override_conc { dinconc } @no_store
 
 	var(bnd_layer.water.on, [k g], [m g, l-1], "Boundary layer DON")
-		@override_conc { docconc / cn_molar_to_mass_ratio(bnd_cn)  }
+		@override_conc { docconc / cn_molar_to_mass_ratio(bnd_cn)  } @no_store
 
 	var(bnd_layer.water.phos, [k g], [m g, l-1], "Boundary layer DIP")
-		@override_conc { dipconc }
+		@override_conc { dipconc } @no_store
 
 	var(bnd_layer.water.op, [k g], [m g, l-1], "Boundary layer DOP")
-		@override_conc { docconc / cp_molar_to_mass_ratio(bnd_cp) }
+		@override_conc { docconc / cp_molar_to_mass_ratio(bnd_cp) } @no_store
 
 	var(bnd_layer.water.phyt, [k g], [m g, l-1], "Boundary layer phytoplankton")
-		@override_conc { chl_a / chl_a_f ->> }
+		@override_conc { chl_a / chl_a_f ->> } @no_store
 
 	var(bnd_layer.water.sed, [k g], [m g, l-1], "Boundary layer SS")
-		@override_conc { pocconc / pocfrac }
+		@override_conc { pocconc / pocfrac } @no_store
 
 	var(bnd_layer.water.sed.oc, [k g], "Boundary layer POC")
-		@override_conc { pocfrac } # This is the concentration in water.sed
+		@override_conc { pocfrac } @no_store # This is the concentration in water.sed
 
 	var(bnd_layer.water.sed.on, [k g], "Boundary layer PON")
-		@override_conc { pocfrac / cn_molar_to_mass_ratio(bnd_cn) }
+		@override_conc { pocfrac / cn_molar_to_mass_ratio(bnd_cn) } @no_store
 
 	var(bnd_layer.water.sed.op, [k g], "Boundary layer POP")
-		@override_conc { pocfrac / cp_molar_to_mass_ratio(bnd_cp) }
+		@override_conc { pocfrac / cp_molar_to_mass_ratio(bnd_cp) } @no_store
 
 	var(bnd_layer.water.sed.phos, [k g], "Boudary layer PIP")
-		@override_conc { 0=>> }  #TODO: What should this be?
+		@override_conc { 0=>> } @no_store #TODO: What should this be?
 
 	var(bnd_layer.water.phyt.oc, [k g], "Boundary layer phytoplankton C")
-		@override_conc { 1 }
+		@override_conc { 1 } @no_store
 
 	var(bnd_layer.water.phyt.on, [k g], "Boundary layer phytoplankton N")
-		@override_conc { 1 / cn_molar_to_mass_ratio(phyt_cn) }
+		@override_conc { 1 / cn_molar_to_mass_ratio(phyt_cn) } @no_store
 
 	var(bnd_layer.water.phyt.op, [k g], "Boundary layer phytoplankton P")
-		@override_conc { 1 / cn_molar_to_mass_ratio(phyt_cp) }
+		@override_conc { 1 / cn_molar_to_mass_ratio(phyt_cp) } @no_store
 
 	# TODO: CO2 and CH4
 

models/modules/nivafjord/boundary_basin.txt

@@ -31,8 +31,9 @@ module("NIVAFjord boundary", version(0, 0, 1),
 	} @no_store
 
 	# These are needed for automatic transport with fluxes going from the boundary into a basin.
-	var(bnd_layer.water.heat, [J], "Boundary layer thermal energy")     @override_conc { C_water * (temp -> [K]) * rho_water } #See note in basin.txt about maybe having salinity dependence of C_water
-	var(bnd_layer.water.salt, [k g], [k g, m-3], "Boundary layer salt") @override_conc { 1e-3*salinity*rho_water }
+	# See note in basin.txt about maybe having salinity dependence of C_water
+	var(bnd_layer.water.heat, [J], "Boundary layer thermal energy")     @override_conc { C_water * (temp -> [K]) * rho_water } @no_store
+	var(bnd_layer.water.salt, [k g], [k g, m-3], "Boundary layer salt") @override_conc { 1e-3*salinity*rho_water } @no_store
 
 	sea_level_component : function(amp : [c m], phase : [deg], per : [hr], t : [hr]) {
 		amp * cos( (pi/180[deg])*(t*360[deg]/per - phase) )

models/modules/simply_glacier.txt

@@ -68,8 +68,8 @@ Authors: Leah A. Jackson-Blake, Magnus D. Norling
 		0                      otherwise
 	}
 
-	flux(glacier.water, soil.water, [m m, day-1], "Input to soil water") {
+	flux(glacier.water, soil.water, [m m, day-1], "Glacial input to soil water") {
 		perc_frac*flow     if ice > ice_epsilon,
-		flow                  otherwise
+		flow               otherwise
 	}
 }

src/function_tree.cpp

@@ -786,8 +786,7 @@ resolve_special_directive(Function_Call_AST *ast, Directive directive, Function_
 			fatal_error_trace(scope);
 		}
 		ident->var_id = conc_id;
-		//if(conc_id.type == Var_Id::Type::series)  // Annoying that this is not automatic..
-		//	ident->variable_type = Variable_Type::series;
+
 		result.unit = data->app->vars[conc_id]->unit.standard_form;
 	} else
 		result.unit = std::move(arg_units[var_idx]);
@@ -801,6 +800,8 @@ resolve_tuple(Function_Call_AST *ast, Function_Resolve_Data *data, Function_Scop
 	new_tuple->value_type = Value_Type::tuple;
 	resolve_arguments(new_tuple, ast, data, scope, new_tuple->element_units);
 
+	arguments_must_be_values(new_tuple, scope);
+
 	result.fun = new_tuple;
 }
 

src/model_compilation.cpp

@@ -913,9 +913,10 @@ basic_instruction_solver_configuration(Model_Application *app, std::vector<Model
 				}
 				continue;
 			}
-			if(!app->vars[var_id]->store_series) {
+			auto var = app->vars[var_id];
+			if(!var->store_series && var->is_mass_balance_quantity()) {
 				source_loc.print_error_header(Mobius_Error::model_building);
-				fatal_error("The variable \"", app->vars[var_id]->name, "\" was specified as @no_store, and so can not be put on a solver.");
+				fatal_error("The variable \"", app->vars[var_id]->name, "\" was specified as @no_store, and so can not be put on a solver (unless it is @override).");
 			}
 			instructions[var_id.id].solver = solver_id;
 		}
@@ -1956,8 +1957,17 @@ Model_Application::compile(bool store_code_strings) {
 				else
 					vars.push_back(var);
 			}
-			build_batch_arrays(this, vars,     instructions, batch.arrays,     false);
-			build_batch_arrays(this, vars_ode, instructions, batch.arrays_ode, false);
+
+			if(vars_ode.empty()) {
+				// In some rare cases we could end up with no ODEs in a batch that was assigned a solver. In that case we can demote it to a discrete batch.
+				batch.solver = invalid_entity_id;
+				for(auto instr_id : batch.instrs)
+					instructions[instr_id].solver = invalid_entity_id;
+				build_batch_arrays(this, batch.instrs, instructions, batch.arrays, false);
+			} else {
+				build_batch_arrays(this, vars,     instructions, batch.arrays,     false);
+				build_batch_arrays(this, vars_ode, instructions, batch.arrays_ode, false);
+			}
 		}
 	}
 

src/model_composition.cpp

@@ -701,7 +701,10 @@ create_conc_var(Model_Application *app, Var_Id var_id, Var_Id dissolved_in_id, E
 	auto dissolved_in = app->vars[dissolved_in_id];
 
 	sprintf(varname, "concentration(%s, %s)", var_name.data(), dissolved_in->name.data());
-	Var_Id gen_conc_id = register_state_variable<State_Var::Type::dissolved_conc>(app, invalid_entity_id, false, varname);
+
+	// Default the concentration variable to have the same storage setting as the main variable
+	bool store = app->vars[var_id]->store_series;
+	Var_Id gen_conc_id = register_state_variable<State_Var::Type::dissolved_conc>(app, invalid_entity_id, false, varname, !store);
 	auto conc_var = as<State_Var::Type::dissolved_conc>(app->vars[gen_conc_id]);
 	conc_var->conc_of = var_id;
 	conc_var->conc_in = dissolved_in_id;

test/models/another_simple.txt

@@ -1,27 +1,22 @@
 
 
 model("Another simple") {
-
-	i : index_set("I")
-	j : index_set("J")
-	k : index_set("K") @union(i, j)
 
-	a : compartment("A", k)
-	b : compartment("B", j)
+	a : compartment("A")
+	q : quantity("Q")
 
 	module("So simple", version(2, 0, 0)) {
-
-		par_group("Pars", a) {
-			p : par_real("P", [], 5)
-		}
-
-		t : property("Test 1")
-		var(a.t, []) { p }
 
-		t3 : property("Test 3")
-		var(b.t3, []) { 50 }
+		var(a.q, [], "A Queue!") @override {
+			time.step=>>
+		} @no_store
 
-		t2 : property("Test 2")
-		var(b.t2, []) { a.t }
+		r : property("P")
+		var(a.r, [], "Arrr!") { a.q*2 }
+
 	}
+
+	sol : solver("solver", inca_dascru, [2, hr], 0.01)
+
+	solve(sol, a.q)
 }

test/models/another_simple_data.dat

@@ -4,9 +4,6 @@ data_set {
 The doc string
 """
 
-	i : index_set("I") [ 1 ]
-	j : index_set("J") [ 1 ]
-
 	par_group("System") {
 
 		par_datetime("Start date")