Merge branch 'main' of https://github.com/Schulik/aiolos into main

Schulik · Oct 15, 2020 · 926b61f · 926b61f
2 parents 88896ad + 0435f43
commit 926b61f
Show file tree

Hide file tree

Showing 22 changed files with 559 additions and 176 deletions.
diff --git a/.gitignore b/.gitignore
@@ -0,0 +1,10 @@
+*.o
+*.dat
+*~
+test_files/*.o
+test_files/*.dat
+test_files/*.png
+test_files/*~
+aiolos
+tests
+.vscode
diff --git a/advection.cpp b/advection.cpp
@@ -1,4 +1,4 @@
-#include "main.h"
+#include "aiolos.h"
 
 void hydro_run::execute() { 
 
@@ -39,13 +39,7 @@ void hydro_run::execute() {
         // Step 1: Boundary values
         //
         //
-
-        if(boundaries_number == 4) {
-            phi[num_cells+1] = get_phi_grav(x_i12[num_cells], planet_mass);
-        } else {
-            phi[num_cells+1] = get_phi_grav(x_i12[num_cells+1], planet_mass);
-        }
-
+
         apply_boundary_left() ;
         apply_boundary_right() ;
 

diff --git a/main.h → aiolos.h b/main.h → aiolos.h
@@ -191,7 +191,7 @@ class hydro_run
     double planet_mass;     //in Earth masses
     double planet_position; //inside the simulation domain
     double rs;
-    double rs_at_moment;
+    double rs_at_moment = 0 ;
     double rs_time;
     int    init_static_atmosphere;
     int    static_atmosphere_tempprofile;
@@ -298,11 +298,6 @@ class hydro_run
     //
     // Boundaries
     //
-
-    void boundaries_const_both(AOS &left_ghost, const AOS &leftval, const AOS &rightval, AOS &right_ghost );
-    void boundaries_open_both(AOS &left_ghost, const AOS &leftval, const AOS &leftval2, const AOS &rightval2, const AOS &rightval, AOS &right_ghost );
-    void boundaries_planet_mdot(AOS &left_ghost, const AOS &leftval, const AOS &rightval, AOS &right_ghost );
-    void boundaries_wall_both(AOS &left_ghost, const AOS &leftval, const AOS &rightval, AOS &right_ghost );
     void add_wave(double time, double tmax, double amplitude);
 
     void apply_boundary_left();
@@ -339,7 +334,7 @@ class hydro_run
 
 public:
 
-    hydro_run(string);
+    hydro_run(string, double debug=0);
     ~hydro_run();
 
     void execute(); //Main loop

diff --git a/helpers.cpp b/helpers.cpp
@@ -1,6 +1,6 @@
 #include <sstream>
 #include <stdexcept>
-#include "main.h"
+#include "aiolos.h"
 
 std::vector<AOS> init_AOS(int num) {   
     return std::vector<AOS>(num);
@@ -240,15 +240,29 @@ void hydro_run::print_AOS_component_tofile(const std::vector<double>& x,
                                           const std::vector<AOS>& data,
                                           const std::vector<AOS>& fluxes,
                                           int timestepnumber) {
-    string filename;
-    stringstream filenamedummy;
-    string truncated_name = stringsplit(simname,".")[0];
-    filenamedummy<<"output_"<<truncated_name<<"_t"<<timestepnumber<<".dat";
-
+
+    // Choose a sensible default output name
+    string filename ;
+    {
+        stringstream filenamedummy;
+        string truncated_name = stringsplit(simname,".")[0];
+        filenamedummy<<"output_"<<truncated_name;
+        filename = filenamedummy.str() ;
+    }
+
+    filename = read_parameter_from_file<string>(simname, "OUTPUT_FILENAME", debug, filename).value ;
+
+    // Add the snap number
+    {
+        stringstream filenamedummy;
+        filenamedummy << filename << "_t" << timestepnumber << ".dat";
+        filename = filenamedummy.str() ;
+    }
+
     if(debug > 1)
-        cout<<"Trying to open file "<<filenamedummy.str()<<endl;
+        cout<<"Trying to open file "<<filename<<endl;
 
-    ofstream outfile(filenamedummy.str(), ios::out);
+    ofstream outfile(filename, ios::out);
 
     if (outfile.is_open())
     {
@@ -262,8 +276,8 @@ void hydro_run::print_AOS_component_tofile(const std::vector<double>& x,
         //Print the domain
         for(int i=1; i<=num_cells; i++) {
 
-            hydrostat = flux[i].u2/dx[i] ; //hydrostat2 + hydrostat3 ; 
-            hydrostat2 = flux[i+1].u2/dx[i];//pressure[i+1] - pressure[i];
+            hydrostat = flux[i-1].u2/dx[i] ; //hydrostat2 + hydrostat3 ; 
+            hydrostat2 = flux[i].u2/dx[i];//pressure[i+1] - pressure[i];
             hydrostat3 = source[i].u2;//0.5 * (data[i].u1 + data[i+1].u1) * (phi[i+1] - phi[i]);
 
 
@@ -276,10 +290,11 @@ void hydro_run::print_AOS_component_tofile(const std::vector<double>& x,
 
         //Print right ghost stuff
         outfile<<x[num_cells+1]<<'\t'<<data[num_cells+1].u1<<'\t'<<data[num_cells+1].u2<<'\t'<<data[num_cells+1].u3<<'\t'<<'-'<<'\t'<<'-'<<'\t'<<'-'<<'\t'<<'-'<<'\t'<<'-'<<'\t'<<'-'<<'\t'<<pressure[num_cells+1]<<'\t'<<data[num_cells+1].u2/data[num_cells+1].u1<<'\t'<<'-'<<'\t'<<'-'<<'\t'<<phi[num_cells+1]<<'\t'<<'-'<<'\t'<<'-'<<'\t'<<'-'<<'\t'<<'-'<<'\t'<<'-'<<endl;
+
+        cout<<"Sucessfully written file "<<filename<<" with smoothed gravity = "<<rs_at_moment<<" at time "<<globalTime<<" and dt="<<dt<<endl;
     }
-    else cout << "Unable to open file.";
+    else cout << "Unable to open file" << filename << endl; 
     outfile.close();
 
-    cout<<"Sucessfully written file "<<filenamedummy.str()<<" with smoothed gravity = "<<rs_at_moment<<" at time "<<globalTime<<" and dt="<<dt<<endl;
 
 }
diff --git a/init_and_bounds.cpp b/init_and_bounds.cpp
@@ -1,4 +1,4 @@
-#include "main.h"
+#include "aiolos.h"
 
 //extern AOS* init_AOS(int num);
 
@@ -11,7 +11,9 @@ Initializer for a new simulation. Stores all parameters in the object, given to
     
         par_physics: List of physics parameters
 */
-hydro_run::hydro_run(string filename) {
+hydro_run::hydro_run(string filename, double debug_) {
+
+        debug = debug_ ;
 
         ////~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
         //
@@ -340,12 +342,9 @@ hydro_run::hydro_run(string filename) {
             //Conversion from shock tube parameters (given as dens, velocity, pressure) to conserved variables (dens, momentum, internal energy)
             BACKGROUND_U = AOS(u1, u1*u2, 0.5*u1*u2*u2 + u3/(gamma_adiabat-1.) );
 
-            if (boundaries_number == 1) {
-                boundary_left = BoundaryType::reflecting ;
-                boundary_right = BoundaryType::fixed ;
-                SHOCK_TUBE_UR = BACKGROUND_U ;
-            }
-            //mdot boundaries
+
+            /* mdot boundaries
+             * needs fixing... Should just be specified as a BoundaryType::fixed boundary
             if(boundaries_number == 3) {
                 mdot = read_parameter_from_file<double>(filename,"ACCRETION_RATE", debug).value; //Accretion rate in numerical units
     
@@ -357,13 +356,21 @@ hydro_run::hydro_run(string filename) {
                 boundary_left = BoundaryType::reflecting ;
                 boundary_right = BoundaryType::fixed ;
             }
+            */
 
 
             cout<<"Problem 2 pos3."<<endl;
 
             initialize_background(BACKGROUND_U);
             init_grav_pot();
 
+            //
+            // If we want to initialize a hydrostatic atmosphere, then we overwrite the so far given density/pressure data and set every velocity to 0
+            //
+            if(init_static_atmosphere == 1) {
+                initialize_hydrostatic_atmosphere_nonuniform();
+            }
+
             if(debug > 0) 
                 cout<<"Successfully initialized problem 2."<<endl;
         }
@@ -383,14 +390,8 @@ hydro_run::hydro_run(string filename) {
         speed           = np_zeros(num_cells+2);
         cs              = np_zeros(num_cells+2);
 
-        //
-        // If we want to initialize a hydrostatic atmosphere, then we overwrite the so far given density/pressure data and set every velocity to 0
-        //
-        if(init_static_atmosphere == 1) {
-            initialize_hydrostatic_atmosphere_nonuniform();
-        }
 
-        USE_WAVE = read_parameter_from_file<int>(filename,"WAVE_USE", debug).value; 
+        USE_WAVE = read_parameter_from_file<int>(filename,"WAVE_USE", debug, 0).value; 
         if(USE_WAVE==1) {
             WAVE_AMPLITUDE = read_parameter_from_file<double>(filename,"WAVE_AMPLITUDE", debug).value; 
             WAVE_PERIOD    = read_parameter_from_file<double>(filename,"WAVE_PERIOD", debug).value; 
@@ -523,7 +524,7 @@ void hydro_run::initialize_shock_tube_test(const AOS &leftval,const AOS &rightva
 
     for(int i=0; i<=num_cells+1; i++) {
 
-        if(x_i[i] < SHOCK_TUBE_MID) 
+        if(x_i12[i] < SHOCK_TUBE_MID) 
             u[i] = leftval;
         else
             u[i] = rightval;
@@ -541,119 +542,6 @@ void hydro_run::initialize_background(const AOS &background) {
         u[i] = background;
 }
 
-//
-// Non-changing boundary conditions. Give values as leftval and rightval (for example the initial values from SHOCK_TUBE_UL and UR)
-//                                   and they are assigned to the ghost zones
-//
-//                                  In the case a planet is involved, the planet is by default assumed to be sitting at x=0 and 
-//                                  then the left boundary is WALL and the right is CONST.
-void hydro_run::boundaries_const_both(AOS &left_ghost, const AOS &leftval, const AOS &rightval, AOS &right_ghost ){
-
-    if(debug >= 3)
-        cout<<"Inside boundaries_const_both. Assigning values left = "<<leftval.u1<<" "<<leftval.u2<<" "<<leftval.u3<<" and right = "<<leftval.u1<<" "<<rightval.u2<<" "<<rightval.u3<<endl;
-
-    //For the shock tube test we use the initial values as boundaries
-    if(problem_number == 1) {
-        left_ghost  = SHOCK_TUBE_UL;
-        right_ghost = SHOCK_TUBE_UR;
-    } 
-    else if(problem_number == 2) {
-
-        left_ghost = AOS( leftval.u1, -leftval.u2, leftval.u3); 
-        right_ghost= BACKGROUND_U;
-
-    } else {
-        left_ghost  = leftval;
-        right_ghost = rightval;
-    }
-}
-
-
-
-
-//
-// Open boundaries on both ends: Takes the last and second last values on each side and extrapolates them to generate the ghost values
-//
-//                                  In the case a planet is involved, the planet is by default assumed to be sitting at x=0 and 
-//                                  then the left boundary is WALL and the right is OPEN.
-void hydro_run::boundaries_open_both(AOS &left_ghost, const AOS &leftval, const AOS &leftval2, const AOS &rightval2, const AOS &rightval, AOS &right_ghost ){
-
-    if(debug >= 3)
-        cout<<"Inside boundaries_open_both. Assigning values left = "<<leftval.u1<<" "<<leftval.u2<<" "<<leftval.u3<<" and right = "<<leftval.u1<<" "<<rightval.u2<<" "<<rightval.u3<<endl;
-
-    double dl1 = (leftval2.u1 - leftval.u1); /// (x_i[1] - x_i[0])
-    double dl2 = (leftval2.u2 - leftval.u2); /// (x_i[1] - x_i[0])
-    double dl3 = (leftval2.u3 - leftval.u3); /// (x_i[1] - x_i[0])
-    double dr1 = (rightval.u1 - rightval2.u1); /// (x_i[cell_number] - x_i[cell_number-1])
-    double dr2 = (rightval.u2 - rightval2.u2); // (x_i[cell_number] - x_i[cell_number-1])
-    double dr3 = (rightval.u3 - rightval2.u3); // (x_i[cell_number] - x_i[cell_number-1])
-
-
-    if(problem_number == 2) {
-
-        //double small_momentum = 1e-10;
-        //double small_dens     = 1e-10;
-        double tmp_r_dens     = rightval.u1 + dr1;
-        //double tmp_r_momentum = rightval.u2 + dr2;
-
-        //Wall left
-        left_ghost = AOS( leftval.u1, -leftval.u2, leftval.u3); 
-
-        //Inflow/outflow boundary conditions, that do not allow sign change across the last cell boundary, otherwise horrible things happen
-
-        //Sign change detect
-        //if((rightval.u2 + dr2)/rightval.u2 < 0) {
-        //    tmp_r_momentum = rightval.u2;
-            //tmp_r_momentum = small_momentum;
-        //}
-        if(tmp_r_dens < 0.) {
-            tmp_r_dens = rightval.u1;
-
-        }
-        double press_boundary   = (gamma_adiabat-1.) * (rightval.u3 - 0.5 * rightval.u2 * rightval.u2 / rightval.u1 );
-        double hydrostat_energy = (press_boundary - rightval.u1 * (phi[num_cells+1] - phi[num_cells]) ) / (gamma_adiabat-1.);
-
-
-        //right_ghost= AOS(tmp_r_dens, tmp_r_momentum, rightval.u3 + dr3 );
-        //right_ghost = AOS(rightval.u1, rightval.u2, rightval.u3);
-        right_ghost = AOS(rightval.u1, rightval.u2, hydrostat_energy);
-    }
-    else {
-
-        left_ghost  = AOS( leftval.u1 - dl1, leftval.u2  - dl2,  leftval.u3 - dl3 ); 
-        right_ghost = AOS(rightval.u1 + dr1, rightval.u2 + dr2, rightval.u3 + dr3 );
-
-    }
-
-}
-
-void hydro_run::boundaries_planet_mdot(AOS &left_ghost, const AOS &leftval, const AOS &rightval, AOS &right_ghost ) {
-
-    if(debug >= 3)
-        cout<<"Inside boundaries_planet_mdot. Assigning values left = "<<leftval.u1<<" "<<leftval.u2<<" "<<leftval.u3<<" and right = "<<leftval.u1<<" "<<rightval.u2<<" "<<rightval.u3<<endl;
-
-    left_ghost  = AOS( leftval.u1, -leftval.u2, leftval.u3); 
-
-    double rho0 = BACKGROUND_U.u1;
-    double e0   = BACKGROUND_U.u3 - 0.5 * BACKGROUND_U.u2 * BACKGROUND_U.u2 / rho0;
-
-    right_ghost = AOS(rho0, -mdot, e0 + 0.5 * mdot * mdot / rho0 );
-    //we compute back and forth between the two momenta and E's because the initial conditions could be 0 momentum, then we would have mdot=0    
-
-}
-
-//
-// Wall boundary conditions on both. 
-//                                   
-
-void hydro_run::boundaries_wall_both(AOS &left_ghost, const AOS &leftval, const AOS &rightval, AOS &right_ghost ){
-
-    if(debug >= 3)
-        cout<<"Inside boundaries_wall_both. Assigning values left = "<<leftval.u1<<" "<<leftval.u2<<" "<<leftval.u3<<" and right = "<<leftval.u1<<" "<<rightval.u2<<" "<<rightval.u3<<endl;
-
-    left_ghost = AOS( leftval.u1, -leftval.u2, leftval.u3); 
-    right_ghost= AOS( rightval.u1, -rightval.u2, rightval.u3);
-}
 
 void hydro_run::apply_boundary_left() {
     double E_kinetic, pressure_active, pressure_bound ;

diff --git a/main.cpp b/main.cpp
@@ -1,4 +1,4 @@
-#include "main.h"
+#include "aiolos.h"
 
 
 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -73,10 +73,9 @@ int main(int argc, char** argv)
         cout<<"In main, construction of simulation is about to start."<<endl;
 
         //Main simulation class object, is initialized with the simulation parameters from a file
-        hydro_run simulation1(simulationname);
+        hydro_run simulation1(simulationname, debug);
 
         simulation1.set_suppress_warnings(suppress_warnings_global);
-        simulation1.set_debug(debug);
 
         cout<<"In main, execution is about to start."<<endl;
 

diff --git a/makefile b/makefile
@@ -9,15 +9,27 @@ SRC = $(wildcard *.cpp)
 OBJ = $(SRC:.cpp=.o)
 INC = $(SRC:.cpp=.h)
 
+TEST_OBJ = $(subst main.o, test_files/main.o, $(OBJ))
+
+
+
 #euler_hllc: $(OBJ)
 #	$(CXX) -o $@ $(OBJ) $(CXXFLAGS) $(LDFLAGS) $(BFLAGS)
 
 aiolos: $(OBJ) makefile
 	$(CXX) -o $@ $(OBJ) $(CXXFLAGS) $(LDFLAGS)
 	#./$@
 
-%.o: %.cpp %.h makefile
-	$(CXX) $(CPPFLAGS) $(CXXFLAGS) -c $< $(BFLAGS)
+%.o: %.cpp makefile aiolos.h
+	$(CXX) $(CPPFLAGS) $(CXXFLAGS) -c $< $(BFLAGS) -o $@ 
+
+.PHONY: tests
+
+tests: $(TEST_OBJ) makefile aiolos.h
+	$(CXX) -o  $@  $(TEST_OBJ) $(CXXFLAGS) $(LDFLAGS)
+	./tests > /dev/null
+	cd test_files ; python test_shock_tube.py
+	cd test_files ; python test_steady_state.py
 
 clean:
-	rm *.o
+	rm *.o test_files/*.o
diff --git a/source.cpp b/source.cpp
@@ -1,4 +1,4 @@
-#include "main.h"
+#include "aiolos.h"
Original file line number	Diff line number	Diff line change
		@@ -1,4 +1,4 @@
		#include "main.h"
		#include "aiolos.h"



Expand Down