Radiation merge3 fixes

.gitignore

@@ -8,3 +8,4 @@ spack-build*
 .clangd/*
 .cache
 compile_commands.json
+build_octotiger.sh

CMakeLists.txt

@@ -259,6 +259,7 @@ set(header_files
     octotiger/grid_scf.hpp
     octotiger/interaction_types.hpp
     octotiger/lane_emden.hpp
+    octotiger/matrix.hpp
     octotiger/node_client.hpp
     octotiger/node_location.hpp
     octotiger/node_registry.hpp

frontend/init_methods.cpp

@@ -358,6 +358,23 @@ void init_problem(void) {
         //		opts().gravity = false;
         set_problem(radiation_test_problem);
         set_refine_test(radiation_test_refine);
+    } else if (opts().problem == RADIATION_TEST) {
+		assert(opts().radiation);
+//		opts().gravity = false;
+		set_problem(radiation_test_problem);
+		set_refine_test(radiation_test_refine);
+	} else if (opts().problem == RADIATION_DIFFUSION) {
+		assert(opts().radiation);
+			//		opts().gravity = false;
+		set_problem (radiation_diffusion_test_problem);
+		set_refine_test(radiation_test_refine);
+		set_analytic(radiation_diffusion_analytic);
+	} else if (opts().problem == RADIATION_COUPLING) {
+		assert(opts().radiation);
+			//		opts().gravity = false;
+		set_problem (radiation_coupling_test_problem);
+		set_refine_test(radiation_test_refine);
+//				set_analytic(radiation_diffusion_analytic);
     } else if (opts().problem == DWD) {
         opts().n_species = 5;
         set_problem(scf_binary);

frontend/main.cpp

@@ -281,3 +281,4 @@ int main(int argc, char* argv[]) {
 #endif
 }
 #endif
+

octotiger/const.hpp

@@ -13,7 +13,7 @@ const double pi = 3.1415926535897932384e0, eulercon = 0.577215664901532861e0, a2
 
 // physical constants
 const double g = 6.6742867e-8, h = 6.6260689633e-27, hbar = 0.5e0 * h / pi, qe = 4.8032042712e-10,
-		avo = 6.0221417930e23, clight = 2.99792458e10, kerg = 1.380650424e-16, ev2erg = 1.60217648740e-12, kev = kerg
+		avo = 6.0221417930e23,  = 2.99792458e10, kerg = 1.380650424e-16, ev2erg = 1.60217648740e-12, kev = kerg
 				/ ev2erg, amu = 1.66053878283e-24, mn = 1.67492721184e-24, mp = 1.67262163783e-24,
 		me = 9.1093821545e-28, rbohr = hbar * hbar / (me * qe * qe), fine = qe * qe / (hbar * clight), hion =
 				13.605698140e0;

octotiger/defs.hpp

@@ -37,7 +37,7 @@ constexpr int REFINE_BW = 2;
 
 
 //#define EXPERIMENT
-constexpr int NRF = 7;
+constexpr int NRF = 4;
 
 #define abort_error() printf( "Error in %s on line %i\n", __FILE__, __LINE__); abort()
 
@@ -92,7 +92,7 @@ constexpr real DEFAULT_OMEGA = 0.0;
 
 #include "hydro_defs.hpp"
 
-constexpr integer R_BW = 2;
+constexpr integer R_BW = 3;
 
 constexpr integer HS_NX = (2 * H_BW + INX/2);
 constexpr integer G_NX = INX;
@@ -117,8 +117,12 @@ constexpr integer G_DNZ = 1;
 constexpr integer G_DN[NDIM] = { G_NX * G_NX, G_NX, 1 };
 
 // Radiation {{{
+constexpr integer RAD_NF = 4;
 constexpr integer RAD_BW = 3;
 constexpr integer RAD_NX = INX + 2 * RAD_BW;
+constexpr integer RAD_INTERIOR_BEGIN = RAD_BW;
+constexpr integer RAD_INTERIOR_END = (RAD_NX - RAD_BW);
+constexpr integer RAD_INTERIOR = (RAD_INTERIOR_END - RAD_INTERIOR_BEGIN);
 constexpr integer RAD_N3 = RAD_NX * RAD_NX * RAD_NX;
 
 constexpr inline integer rindex(integer x, integer y, integer z)

octotiger/diagnostics.hpp

@@ -59,13 +59,17 @@ struct diagnostics_t {
 	std::array<safe_real, NDIM> lsum;
 	safe_real nonvacj;
 	safe_real nonvacjlz;
+	std::vector<std::pair<real,std::vector<real>>> xline;
+	safe_real Trad0;
+	safe_real Tgas0;
 	diagnostics_t() {
 		failed = false;
 		stage = 1;
 		omega = -1.0;
 		grid_com = 0.0;
 		munbound1 = 0.0;
 		munbound2 = 0.0;
+		Trad0 = Tgas0 =0.0;
 		for (integer f = 0; f != opts().n_fields; ++f) {
 			grid_sum[f] = 0.0;
 			grid_out[f] = 0.0;
@@ -153,6 +157,13 @@ struct diagnostics_t {
 				}
 			}
 		}
+		xline.insert(xline.end(),other.xline.begin(),other.xline.end());
+		if( other.Tgas0 ) {
+			Tgas0 = other.Tgas0;
+		}
+		if( other.Trad0 ) {
+			Trad0 = other.Trad0;
+		}
 		munbound1 += other.munbound1;
 		munbound2 += other.munbound2;
 		lsum[0] += other.lsum[0];