Change the behaviour in the absence of unit system defined in the IC file. We...

Change the behaviour in the absence of unit system defined in the IC file. We do not default to CGS any more, but us the internal units instead.

src/common_io.c

@@ -250,13 +250,18 @@ void io_write_attribute_s(hid_t grp, const char* name, const char* str) {
 
 /**
  * @brief Reads the Unit System from an IC file.
+ *
+ * If the 'Units' group does not exist in the ICs, we will use the internal
+ * system of units.
+ *
  * @param h_file The (opened) HDF5 file from which to read.
- * @param us The unit_system to fill.
+ * @param ic_units The unit_system to fill.
+ * @param internal_units The internal system of units to copy if needed.
  * @param mpi_rank The MPI rank we are on.
- *
- * If the 'Units' group does not exist in the ICs, cgs units will be assumed
  */
-void io_read_unit_system(hid_t h_file, struct unit_system* us, int mpi_rank) {
+void io_read_unit_system(hid_t h_file, struct unit_system* ic_units,
+                         const struct unit_system* internal_units,
+                         int mpi_rank) {
 
   /* First check if it exists as this is *not* required. */
   const htri_t exists = H5Lexists(h_file, "/Units", H5P_DEFAULT);
@@ -264,16 +269,12 @@ void io_read_unit_system(hid_t h_file, struct unit_system* us, int mpi_rank) {
   if (exists == 0) {
 
     if (mpi_rank == 0)
-      message("'Units' group not found in ICs. Assuming CGS unit system.");
+      message("'Units' group not found in ICs. Assuming internal unit system.");
 
-    /* Default to CGS */
-    us->UnitMass_in_cgs = 1.;
-    us->UnitLength_in_cgs = 1.;
-    us->UnitTime_in_cgs = 1.;
-    us->UnitCurrent_in_cgs = 1.;
-    us->UnitTemperature_in_cgs = 1.;
+    units_copy(ic_units, internal_units);
 
     return;
+
   } else if (exists < 0) {
     error("Serious problem with 'Units' group in ICs. H5Lexists gives %d",
           exists);
@@ -284,15 +285,15 @@ void io_read_unit_system(hid_t h_file, struct unit_system* us, int mpi_rank) {
 
   /* Ok, Read the damn thing */
   io_read_attribute(h_grp, "Unit length in cgs (U_L)", DOUBLE,
-                    &us->UnitLength_in_cgs);
+                    &ic_units->UnitLength_in_cgs);
   io_read_attribute(h_grp, "Unit mass in cgs (U_M)", DOUBLE,
-                    &us->UnitMass_in_cgs);
+                    &ic_units->UnitMass_in_cgs);
   io_read_attribute(h_grp, "Unit time in cgs (U_t)", DOUBLE,
-                    &us->UnitTime_in_cgs);
+                    &ic_units->UnitTime_in_cgs);
   io_read_attribute(h_grp, "Unit current in cgs (U_I)", DOUBLE,
-                    &us->UnitCurrent_in_cgs);
+                    &ic_units->UnitCurrent_in_cgs);
   io_read_attribute(h_grp, "Unit temperature in cgs (U_T)", DOUBLE,
-                    &us->UnitTemperature_in_cgs);
+                    &ic_units->UnitTemperature_in_cgs);
 
   /* Clean up */
   H5Gclose(h_grp);

src/common_io.h

@@ -75,7 +75,9 @@ void io_write_attribute_s(hid_t grp, const char* name, const char* str);
 void io_write_code_description(hid_t h_file);
 void io_write_engine_policy(hid_t h_file, const struct engine* e);
 
-void io_read_unit_system(hid_t h_file, struct unit_system* us, int mpi_rank);
+void io_read_unit_system(hid_t h_file, struct unit_system* ic_units,
+                         const struct unit_system* internal_units,
+                         int mpi_rank);
 void io_write_unit_system(hid_t h_grp, const struct unit_system* us,
                           const char* groupName);
 

src/parallel_io.c

@@ -709,7 +709,7 @@ void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
   struct unit_system* ic_units =
       (struct unit_system*)malloc(sizeof(struct unit_system));
   if (ic_units == NULL) error("Unable to allocate memory for IC unit system");
-  io_read_unit_system(h_file, ic_units, mpi_rank);
+  io_read_unit_system(h_file, ic_units, internal_units, mpi_rank);
 
   /* Tell the user if a conversion will be needed */
   if (mpi_rank == 0) {

src/serial_io.c

@@ -503,7 +503,7 @@ void read_ic_serial(char* fileName, const struct unit_system* internal_units,
 
     /* Read the unit system used in the ICs */
     if (ic_units == NULL) error("Unable to allocate memory for IC unit system");
-    io_read_unit_system(h_file, ic_units, mpi_rank);
+    io_read_unit_system(h_file, ic_units, internal_units, mpi_rank);
 
     if (units_are_equal(ic_units, internal_units)) {
 

src/single_io.c

@@ -410,7 +410,7 @@ void read_ic_single(const char* fileName,
   struct unit_system* ic_units =
       (struct unit_system*)malloc(sizeof(struct unit_system));
   if (ic_units == NULL) error("Unable to allocate memory for IC unit system");
-  io_read_unit_system(h_file, ic_units, 0);
+  io_read_unit_system(h_file, ic_units, internal_units, 0);
 
   /* Tell the user if a conversion will be needed */
   if (units_are_equal(ic_units, internal_units)) {

src/units.c

@@ -128,6 +128,21 @@ void units_init_default(struct unit_system* us, struct swift_params* params,
       parser_get_opt_param_double(params, buffer, def->UnitTemperature_in_cgs);
 }
 
+/**
+ * @brief Copy the content of a #unit_system to another one.
+ *
+ * @param dest The destination of the copy.
+ * @param src The source of the copy.
+ */
+void units_copy(struct unit_system* dest, const struct unit_system* src) {
+
+  dest->UnitMass_in_cgs = src->UnitMass_in_cgs;
+  dest->UnitLength_in_cgs = src->UnitLength_in_cgs;
+  dest->UnitTime_in_cgs = src->UnitTime_in_cgs;
+  dest->UnitCurrent_in_cgs = src->UnitCurrent_in_cgs;
+  dest->UnitTemperature_in_cgs = src->UnitTemperature_in_cgs;
+}
+
 /**
  * @brief Returns the base unit conversion factor for a given unit system
  * @param us The unit_system used

src/units.h

@@ -103,6 +103,7 @@ void units_init_from_params(struct unit_system*, struct swift_params*,
 void units_init_default(struct unit_system* us, struct swift_params* params,
                         const char* category, const struct unit_system* def);
 
+void units_copy(struct unit_system* dest, const struct unit_system* src);
 int units_are_equal(const struct unit_system* a, const struct unit_system* b);
 
 /* Base units */