From fc6d389cbd5222f435d5bcf16aa8e0bcf1c0f99f Mon Sep 17 00:00:00 2001 From: Richard Booth Date: Mon, 3 Jul 2023 09:57:10 +0100 Subject: [PATCH 01/10] fixed boundaries for shock tube problems --- h | 5 ----- init_and_bounds.cpp | 53 +++++---------------------------------------- 2 files changed, 6 insertions(+), 52 deletions(-) delete mode 100644 h diff --git a/h b/h deleted file mode 100644 index 254bf2b..0000000 --- a/h +++ /dev/null @@ -1,5 +0,0 @@ -Binary file aiolos matches -Binary file aiolos_chem01 matches -chemistry.cpp:double thermo_get_r(double T, double r1, double delta_stoch); -chemistry.cpp:double thermo_get_r(double T, double r1, double delta_stoch) { -Binary file chemistry.o matches diff --git a/init_and_bounds.cpp b/init_and_bounds.cpp index 08b4931..58cb056 100644 --- a/init_and_bounds.cpp +++ b/init_and_bounds.cpp @@ -1782,40 +1782,6 @@ void c_Species::apply_boundary_left(std::vector& u) { switch(boundary_left) { case BoundaryType::user: user_boundary_left(u); - /*double dm = 0; - - mass_reservoir = 0.; - t_evap = 50.; - latent_heat = 3.34e5*1e7/1e3; //3.34e5 J/kg - p_sat = 1e-3 * std::exp(latent_heat/(Rgas*T*(1.-T_evap/prim[i].temperature))); - - for (int i=num_ghosts; i > 0; i--) { - - dp = prim[2].press - p_sat; - - if(dp < 0) { // Condensation - dm = -1; - } - else{ //Evaporation - - if(mass_reservoir < 0.) { - dp = 0.; - dm = 0; - } - - dm = 1; - } - - mass_reservoir += base->dt * dm; - - AOS_prim prim ; - eos->compute_primitive(&u[i],&prim, 1) ; - prim.pres = prim[2] - dp; - prim.pres = std::max( prim.pres, 0.0) ; - eos->compute_conserved(&prim, &u[i-1], 1) ; - } - - */ break; case BoundaryType::open: for (int i=num_ghosts; i > 0; i--) { @@ -1841,28 +1807,21 @@ void c_Species::apply_boundary_left(std::vector& u) { case BoundaryType::fixed: for (int i=0; i < num_ghosts; i++) { - /*double dens_wall; - if(base->problem_number == 1) - dens_wall = SHOCK_TUBE_UL.u1 * mass_amu; - else { - dens_wall = BACKGROUND_U.u1 * mass_amu; - }*/ - double dens_wall; - AOS_prim prim; if(base->problem_number == 1) - dens_wall = SHOCK_TUBE_UL.u1; + u[i] = SHOCK_TUBE_UL ; else { + AOS_prim prim; prim.density = BACKGROUND_U.u1; if(this->prim[2].speed > 0) prim.speed = this->prim[2].speed; else prim.speed = 0.; prim.temperature = const_T_space; + + eos->update_eint_from_T(&(prim), 1); + eos->update_p_from_eint(&(prim), 1); + eos->compute_conserved(&(prim), &(u[i]), 1) ; //Ghost cell u is fixed after init, only need to update p in prim } - - eos->update_eint_from_T(&(prim), 1); - eos->update_p_from_eint(&(prim), 1); - eos->compute_conserved(&(prim), &(u[i]), 1) ; //Ghost cell u is fixed after init, only need to update p in prim } break; From 65d080706fa7ddd87a490604dfc489b6010f8561 Mon Sep 17 00:00:00 2001 From: Richard Booth Date: Tue, 4 Jul 2023 12:04:14 +0100 Subject: [PATCH 02/10] Fixed 2 species friction --- source.cpp | 210 ++++++++++++++------------------ test_files/friction_2spc.par | 1 + test_files/friction_2spc_an.par | 1 + test_files/load_aiolos.py | 14 +-- test_files/test_drag.py | 2 +- test_files/test_steady_state.py | 2 +- 6 files changed, 104 insertions(+), 126 deletions(-) diff --git a/source.cpp b/source.cpp index cc815bb..41416b9 100644 --- a/source.cpp +++ b/source.cpp @@ -260,41 +260,16 @@ void c_Sim::compute_friction_analytical() { for(int j=0; j <= num_cells+1; j++){ + fill_alpha_basis_arrays(j); + compute_alpha_matrix(j); + alpha = friction_coefficients(0, 1) ; + v1b = species[0].prim[j].speed; v2b = species[1].prim[j].speed; - if(collision_model == 'C') { - alpha = alpha_collision; - } else { - // Physical model - m0 = species[0].mass_amu*amu; - m1 = species[1].mass_amu*amu; - mumass = m0 * m1 / (m0 + m1); - meanT = (m1 * species[0].prim[j].temperature + m0 * species[1].prim[j].temperature) / (m0 + m1); - if (species[0].is_dust_like || species[1].is_dust_like) { - // Use Epstein drag law (Hard sphere model) - double RHO_DUST = 1; - double s0=0, s1=0 ; - if (species[0].is_dust_like) - s0 = std::pow(3*m0/(4*M_PI*RHO_DUST), 1/3.) ; - if (species[0].is_dust_like) - s1 = std::pow(3*m1/(4*M_PI*RHO_DUST), 1/3.) ; - - double A = M_PI*(s0+s1)*(s0+s1) ; - double v_th = std::sqrt(8*kb*meanT/(M_PI * mumass)) ; - - alpha = 4/3. * species[1].prim[j].density * v_th * A/(m0+m1) ; - } else { // Gas-gas collisions - coll_b = 5.0e17 * std::pow(meanT, 0.75) ; // from Zahnle & Kasting 1986 Tab. 1 - alpha = kb * meanT * species[1].prim[j].number_density / (m0 * coll_b) ; // From Burgers book, or Schunk & Nagy. - alpha *= alpha_collision; - } - } - - alphas_sample(j) = alpha; - friction_sample(j) = alphas_sample(j) * (v2b - v1b); eps = species[0].u[j].u1 / species[1].u[j].u1; + f1 = (1. + dt*alpha)/(1. + dt*alpha*(1.+eps)) ; f2 = dt*alpha*eps / (1. + dt * alpha); @@ -537,109 +512,110 @@ void c_Sim::compute_alpha_matrix(int j) { //Called in compute_friction() and com if(steps< 4 && j ==100) cout< sj) alpha_local *= dens_vector(sj) / dens_vector(si) ; - } else { - if(si==sj) { - alpha_local = 0.; - } else { // Physical drag laws - //mtot = mass_vector(si) + mass_vector(sj); - mumass = mass_vector(si) * mass_vector(sj) * inv_totmasses(si,sj); /// (mass_vector(si) + mass_vector(sj)); - mumass_amu = mumass/amu; - meanT = (mass_vector(sj)*temperature_vector(si) + mass_vector(si)*temperature_vector(sj)) * inv_totmasses(si,sj); // / (mass_vector(si) + mass_vector(sj)); //Mean collisional mu and T from Schunk 1980 + for(int si=0; si sj) alpha_local *= dens_vector(sj) / dens_vector(si) ; + } else { + if(si==sj) { + alpha_local = 0.; + } else { // Physical drag laws + //mtot = mass_vector(si) + mass_vector(sj); + mumass = mass_vector(si) * mass_vector(sj) * inv_totmasses(si,sj); /// (mass_vector(si) + mass_vector(sj)); + mumass_amu = mumass/amu; + meanT = (mass_vector(sj)*temperature_vector(si) + mass_vector(si)*temperature_vector(sj)) * inv_totmasses(si,sj); // / (mass_vector(si) + mass_vector(sj)); //Mean collisional mu and T from Schunk 1980 + + if (species[si].is_dust_like || species[sj].is_dust_like) { //One of the collision partners is dust + // Use Epstein drag law (Hard sphere model) + double RHO_DUST = 1; + double s0=0, s1=0 ; + if (species[si].is_dust_like) + s0 = std::pow(3*mass_vector(si)/(4*M_PI*RHO_DUST), 1/3.) ; + if (species[sj].is_dust_like) + s1 = std::pow(3*mass_vector(sj)/(4*M_PI*RHO_DUST), 1/3.) ; - if (species[si].is_dust_like || species[sj].is_dust_like) { //One of the collision partners is dust - // Use Epstein drag law (Hard sphere model) - double RHO_DUST = 1; - double s0=0, s1=0 ; - if (species[si].is_dust_like) - s0 = std::pow(3*mass_vector(si)/(4*M_PI*RHO_DUST), 1/3.) ; - if (species[sj].is_dust_like) - s1 = std::pow(3*mass_vector(sj)/(4*M_PI*RHO_DUST), 1/3.) ; - - double A = M_PI*(s0+s1)*(s0+s1) ; - double v_th = std::sqrt(8*kb*meanT/(M_PI * mumass)) ; + double A = M_PI*(s0+s1)*(s0+s1) ; + double v_th = std::sqrt(8*kb*meanT/(M_PI * mumass)) ; - alpha_local = 4/3. * dens_vector(sj) * v_th * A * inv_totmasses(si,sj) ; - } else { // Gas-gas collisions. Used binary diffusion coefficient + alpha_local = 4/3. * dens_vector(sj) * v_th * A * inv_totmasses(si,sj) ; + } else { // Gas-gas collisions. Used binary diffusion coefficient - int ci = (int)(species[si].static_charge*1.01); - int cj = (int)(species[sj].static_charge*1.01); - double qn = (std::fabs(species[si].static_charge) + std::fabs(species[sj].static_charge)); - //string ccase = ""; //Debug string to make sure the cases are picked right + int ci = (int)(species[si].static_charge*1.01); + int cj = (int)(species[sj].static_charge*1.01); + double qn = (std::fabs(species[si].static_charge) + std::fabs(species[sj].static_charge)); + //string ccase = ""; //Debug string to make sure the cases are picked right + + if(std::abs(ci) + std::abs(cj) == 0) { //n-n collision + coll_b = 5.0e17 * std::sqrt(std::sqrt(meanT*meanT*meanT)) ; //std::pow(meanT, 0.75) // from Zahnle & Kasting 1986 Tab. 1 + alpha_local = kb * meanT * numdens_vector(sj) * species[si].inv_mass / coll_b ; // From Burgers book, or Schunk & Nagy. + //ccase = " n-n "; + } + else if(std::abs(ci) == 0 || std::abs(cj) == 0) {//i-n collision - if(std::abs(ci) + std::abs(cj) == 0) { //n-n collision - coll_b = 5.0e17 * std::sqrt(std::sqrt(meanT*meanT*meanT)) ; //std::pow(meanT, 0.75) // from Zahnle & Kasting 1986 Tab. 1 - alpha_local = kb * meanT * numdens_vector(sj) * species[si].inv_mass / coll_b ; // From Burgers book, or Schunk & Nagy. - //ccase = " n-n "; - } - else if(std::abs(ci) == 0 || std::abs(cj) == 0) {//i-n collision + if(species[si].mass_amu < 1.1 && species[sj].mass_amu < 1.1) { //Resonant H,H+ collisions + alpha_local = 2.65e-10 * numdens_vector(sj) * std::sqrt(meanT) * std::pow((1. - 0.083 * std::log10(meanT)), 2.); - if(species[si].mass_amu < 1.1 && species[sj].mass_amu < 1.1) { //Resonant H,H+ collisions - alpha_local = 2.65e-10 * numdens_vector(sj) * std::sqrt(meanT) * std::pow((1. - 0.083 * std::log10(meanT)), 2.); - - //ccase = "resonant H-p+ "; - - } else { - - alpha_local = 1*2.21 * 3.141592 * numdens_vector(sj) * mass_vector(sj) * inv_totmasses(si,sj); // /(mass_vector(sj)+mass_vector(si)); - alpha_local *= std::sqrt(0.66 / mumass ) * 1e-12 * qn * elm_charge; //0.66 is the polarizability of neutral atomic H - alpha_local *= resonant_pair_matrix(si,sj); //Crude way to emulate resonant collision cross-sections: pairs are multiplied by 10, everything else by 1 - - //ccase = "resonant or nonresonant i-n "; - } - } - else { //i-i collision - alpha_local = alpha_collision_ions * 1.27 * species[si].static_charge * species[si].static_charge * species[sj].static_charge * species[sj].static_charge * std::sqrt(mumass_amu) * amu * species[si].inv_mass; // / mass_vector(si); - alpha_local *= numdens_vector(sj) / std::sqrt(meanT*meanT*meanT); + //ccase = "resonant H-p+ "; - //ccase = " i-i "; + } else { + + alpha_local = 1*2.21 * 3.141592 * numdens_vector(sj) * mass_vector(sj) * inv_totmasses(si,sj); // /(mass_vector(sj)+mass_vector(si)); + alpha_local *= std::sqrt(0.66 / mumass ) * 1e-12 * qn * elm_charge; //0.66 is the polarizability of neutral atomic H + alpha_local *= resonant_pair_matrix(si,sj); //Crude way to emulate resonant collision cross-sections: pairs are multiplied by 10, everything else by 1 + + //ccase = "resonant or nonresonant i-n "; } + } + else { //i-i collision + alpha_local = alpha_collision_ions * 1.27 * species[si].static_charge * species[si].static_charge * species[sj].static_charge * species[sj].static_charge * std::sqrt(mumass_amu) * amu * species[si].inv_mass; // / mass_vector(si); + alpha_local *= numdens_vector(sj) / std::sqrt(meanT*meanT*meanT); - if(si!=sj && debug>3) - cout<<"cell "<3) + cout<<"cell "< Date: Tue, 4 Jul 2023 13:45:56 +0100 Subject: [PATCH 03/10] got coll heating test working --- advection.cpp | 940 ++++++++++++++++--------------- init_and_bounds.cpp | 2 +- test_files/collheat_2spc.par | 1 + test_files/collheat_2spc_rad.par | 1 + test_files/friction_6spc.par | 1 + test_files/load_aiolos.py | 6 +- test_files/test_conservation.py | 2 +- test_files/test_irradiation.py | 2 +- 8 files changed, 510 insertions(+), 445 deletions(-) diff --git a/advection.cpp b/advection.cpp index c1a2191..8c75787 100644 --- a/advection.cpp +++ b/advection.cpp @@ -1,7 +1,7 @@ /** * advection.cpp * - * This file contains the base routines to start and run a simulation, + * This file contains the base routines to start and run a simulation, * as well as hydrodynamics, i.e. the HLLC Riemann solver. */ @@ -9,390 +9,433 @@ /** * Main simulation loop - * + * * After successful simulation init, execute the main simulation from globalTime = [0.,t_max] - * Every loop iteration starts with an estimator for all species' mid-cell pressures, + * Every loop iteration starts with an estimator for all species' mid-cell pressures, * to be able to determine the sound-speed and hence the CFL factor. * Following are output checks and then the main code modules are executed in order. */ -void c_Sim::execute() { - +void c_Sim::execute() +{ + steps = 0; double output_counter = 0; - double monitor_counter= 0; + double monitor_counter = 0; const double dt_initial = dt_min_init; - double next_print_time = dt_initial; + double next_print_time = dt_initial; const signed long long int maxsteps = 1e12; - + int crashed_T = 0, crashed_J = 0; - int crash_T_imin = num_cells+2, crash_T_imax = 0, crash_T_numcells = 0; - int crash_J_imin = num_cells+2, crash_J_imax = 0, crash_J_numcells = 0; + int crash_T_imin = num_cells + 2, crash_T_imax = 0, crash_T_numcells = 0; + int crash_J_imin = num_cells + 2, crash_J_imax = 0, crash_J_numcells = 0; double crashtime, crashed_temperature, crashed_meanintensity; int crashtime_already_assigned = 0; - - cout<<" VERSION 0.2"< 0. && globalTime > start_hydro_time) { //Comment in if a radiative equilibrium phase is desired before starting hydro - do_hydrodynamics = 1; - } - - if(steps==0) { - for(int s = 0; s < num_species; s++) + for (globalTime = 0; (globalTime < t_max) && (steps < maxsteps);) + { + + if (start_hydro_time > 0. && globalTime > start_hydro_time) + { // Comment in if a radiative equilibrium phase is desired before starting hydro + do_hydrodynamics = 1; + } + + if (steps == 0) + { + for (int s = 0; s < num_species; s++) species[s].compute_pressure(species[s].u); compute_total_pressure(); } - + dt = get_cfl_timestep(); - dt = std::min(dt, timestep_rad2) ; - dt = std::min(dt, t_max - globalTime) ; - if(steps == 0) + dt = std::min(dt, timestep_rad2); + dt = std::min(dt, t_max - globalTime); + if (steps == 0) dt = std::min(dt, dt_initial); - - if( globalTime > next_print_time) { - cout<<" Beginning step "< next_print_time) + { + cout << " Beginning step " << steps << " @ globalTime " << globalTime << " dt " << dt; + cout << ", CFL " << cfl_step << ", energy dt " << timestep_rad2 << "\n"; next_print_time *= 10.; } - + // // Save internal energy before we update it // - for(int s = 0; s < num_species; s++) - for(int i=num_cells; i>=0; i--) { + for (int s = 0; s < num_species; s++) + for (int i = num_cells; i >= 0; i--) + { species[s].primlast[i].internal_energy = species[s].prim[i].internal_energy; } - + // // Step 0: Update gravity. Important for self-gravity and time-dependent smoothing length // - if(debug >= 2) - cout<<"Beginning timestep "<= 2) + cout << "Beginning timestep " << steps << endl; + update_mass_and_pot(); - - for(int s=0; s= 2) - cout<<"Before fluxes... "; - + + for (int s = 0; s < num_species; s++) + species[s].update_kzz_and_gravpot(s); // Recomputes the homopause boundary and adjusts species-specific potentials + + if (debug >= 2) + cout << "Before fluxes... "; + // - // Output data, when required. Keep the output here, between the flux and conserved variable update, so that the + // Output data, when required. Keep the output here, between the flux and conserved variable update, so that the // zeroth output has the initialized conserved values, but already the first fluxes. // - - if(steps==0 || steps==99e99) { - for(int s=0; s output_counter*output_time + output_time_offset){ - if(debug >= 1) - cout<<" Globaltime is "< monitor_counter*monitor_time) { + + monitor_counter += 1.; + output_counter += 1.; + } + if (globalTime > output_counter * output_time + output_time_offset) + { + if (debug >= 1) + cout << " Globaltime is " << globalTime << " and comparevalue is " << output_counter << " " << output_time << endl; + + print_diagnostic_file((int)output_counter); + for (int s = 0; s < num_species; s++) + species[s].print_AOS_component_tofile((int)output_counter); + + output_counter += 1.; + } + if (globalTime > monitor_counter * monitor_time) + { print_monitor(steps); - - monitor_counter += 1.; + + monitor_counter += 1.; } - + ////~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//// // // Proper start // Do all explicit and implicit operations for one timestep on the entire grid for all species // ////~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//// - + // // Step 0: Hydrodynamics, if so desired // - if (do_hydrodynamics == 1) { - - for(int s = 0; s < num_species; s++) + if (do_hydrodynamics == 1) + { + + for (int s = 0; s < num_species; s++) species[s].execute(species[s].u, species[s].dudt[0]); - - for(int s=0; s < num_species; s++) { - species[s].u0 = species[s].u ; - for(int j=0; j < num_cells+2; j++) - species[s].u[j] = species[s].u[j] + species[s].dudt[0][j]*dt ; + + for (int s = 0; s < num_species; s++) + { + species[s].u0 = species[s].u; + for (int j = 0; j < num_cells + 2; j++) + species[s].u[j] = species[s].u[j] + species[s].dudt[0][j] * dt; } } - - if (order == IntegrationType::first_order) { + + if (order == IntegrationType::first_order) + { globalTime += dt; } - else if (order == IntegrationType::second_order) { + else if (order == IntegrationType::second_order) + { // 2nd step evaluated at t+dt globalTime += dt; - if(use_self_gravity==1) + if (use_self_gravity == 1) update_mass_and_pot(); - if (do_hydrodynamics == 1) { + if (do_hydrodynamics == 1) + { - if (use_drag_predictor_step) { - for(int s = 0; s < num_species; s++) + if (use_drag_predictor_step) + { + for (int s = 0; s < num_species; s++) species[s].compute_pressure(species[s].u); compute_total_pressure(); - - compute_drag_update() ; + + compute_drag_update(); if (use_collisional_heating) - compute_collisional_heat_exchange() ; //Disable if radiation is used? + compute_collisional_heat_exchange(); // Disable if radiation is used? } - //debug = 2; - for(int s = 0; s < num_species; s++) + // debug = 2; + for (int s = 0; s < num_species; s++) species[s].execute(species[s].u, species[s].dudt[1]); - for(int s = 0; s < num_species; s++){ - for(int j = 0; j < num_cells+2; j++) { + for (int s = 0; s < num_species; s++) + { + for (int j = 0; j < num_cells + 2; j++) + { // Need to re-compute u[j] if it was updated in the drag-predictor if (use_drag_predictor_step) - species[s].u[j] = species[s].u0[j] + species[s].dudt[0][j]*dt ; - - species[s].u[j] += (species[s].dudt[1][j] - species[s].dudt[0][j])*dt / 2 ; + species[s].u[j] = species[s].u0[j] + species[s].dudt[0][j] * dt; + + species[s].u[j] += (species[s].dudt[1][j] - species[s].dudt[0][j]) * dt / 2; } } - } else { - for(int s = 0; s < num_species; s++) { - species[s].apply_boundary_left(species[s].u) ; - species[s].apply_boundary_right(species[s].u) ; + } + else + { + for (int s = 0; s < num_species; s++) + { + species[s].apply_boundary_left(species[s].u); + species[s].apply_boundary_right(species[s].u); } } } - - for(int s = 0; s < num_species; s++) + + for (int s = 0; s < num_species; s++) species[s].compute_pressure(species[s].u); - - //Computes the velocity drag update after the new hydrodynamic state is known for each species - if (do_hydrodynamics == 1) - compute_drag_update() ; - - if (do_hydrodynamics == 0 && friction_solver > 0) - compute_drag_update() ; - + + // Computes the velocity drag update after the new hydrodynamic state is known for each species + if (do_hydrodynamics == 1) + compute_drag_update(); + + if (do_hydrodynamics == 0 && friction_solver > 0) + compute_drag_update(); + // If either switch is set we need to think more carefully about what should be done - if( (photochemistry_level + use_rad_fluxes ) > 0 ) { - + if ((photochemistry_level + use_rad_fluxes) > 0) + { + update_opacities(); - if(photochemistry_level == 0 && use_rad_fluxes > 0) + if (photochemistry_level == 0 && use_rad_fluxes > 0) reset_dS(); - + // Compute high-energy dS and ionization - if(photochemistry_level == 1) { //C2Ray scheme + if (photochemistry_level == 1) + { // C2Ray scheme reset_dS(); - - if(false) { - cout<<"Pos 1 dS_UV = "< 1) + cout << "timestep in execute()=" << dt << " stepnum " << steps << " totaltime" << globalTime << endl; + // // Detection of negative J and T and soft exit // - for(int s = 0; s < num_species; s++) { - for(int i=num_cells-1; i>=0; i--) { - if(species[s].prim[i].temperature < 0 || std::isnan(species[s].prim[i].temperature) ) { - - if(crashtime_already_assigned == 0) { - crashtime = globalTime; - crashtime_already_assigned = 1; - } - - crashed_temperature = species[s].prim[i].temperature; - crashed_T = s+1; - crash_T_imin = (icrash_T_imax)? i : crash_T_imax; - crash_T_numcells++; - - globalTime = 1.1*t_max; //This secures that the program exits smoothly after the crash and produces a -1 output file of the crashed state + for (int s = 0; s < num_species; s++) + { + for (int i = num_cells - 1; i >= 0; i--) + { + if (species[s].prim[i].temperature < 0 || std::isnan(species[s].prim[i].temperature)) + { + + if (crashtime_already_assigned == 0) + { + crashtime = globalTime; + crashtime_already_assigned = 1; } + + crashed_temperature = species[s].prim[i].temperature; + crashed_T = s + 1; + crash_T_imin = (i < crash_T_imin) ? i : crash_T_imin; + crash_T_imax = (i > crash_T_imax) ? i : crash_T_imax; + crash_T_numcells++; + + globalTime = 1.1 * t_max; // This secures that the program exits smoothly after the crash and produces a -1 output file of the crashed state + } + } + + if (crashed_T > 0) + { + cout << endl + << ">>> CRASH <<< DUE TO NEGATIVE TEMPERATURES, crash_imin/imax = " << crash_T_imin << "/" << crash_T_imax << " sample T = " << crashed_temperature << " num of crashed cells/total cells = " << crash_T_numcells << "/" << num_cells << " crashed species name = " << species[crashed_T - 1].speciesname << endl; + cout << " @t/dt = " << crashtime << "/" << dt << " stepnum " << steps << endl; + cout << "Writing crash dump into last output and exiting program." << endl; } - - if(crashed_T > 0) { - cout<>> CRASH <<< DUE TO NEGATIVE TEMPERATURES, crash_imin/imax = "<=0; i--) { - if(Jrad_FLD(i,b) < 0) { - - if(Jrad_FLD(i,b) > -1e-10) { //Ignore extremely small flux overshoots - Jrad_FLD(i,b) *= -1.; + if (couple_J_into_T) + { + for (int b = 0; b < num_bands_out; b++) + { + for (int i = num_cells - 1; i >= 0; i--) + { + if (Jrad_FLD(i, b) < 0) + { + + if (Jrad_FLD(i, b) > -1e-10) + { // Ignore extremely small flux overshoots + Jrad_FLD(i, b) *= -1.; } - else { - + else + { + // If we're computing the planet's surface temperature // allow J < 0 in the boundary. if (i < num_cells && use_planetary_temperature && - Jrad_FLD(num_cells,b) > 0) - continue ; - - if(crashtime_already_assigned == 0) { + Jrad_FLD(num_cells, b) > 0) + continue; + + if (crashtime_already_assigned == 0) + { crashtime = globalTime; crashtime_already_assigned = 1; } - - crashed_meanintensity = Jrad_FLD(i,b); - crashed_J = b+1; - crash_J_imin = (icrash_J_imax)? i : crash_J_imax; + + crashed_meanintensity = Jrad_FLD(i, b); + crashed_J = b + 1; + crash_J_imin = (i < crash_J_imin) ? i : crash_J_imin; + crash_J_imax = (i > crash_J_imax) ? i : crash_J_imax; crash_J_numcells++; - - globalTime = 1.1*t_max; + + globalTime = 1.1 * t_max; } } - } - - if(crashed_J > 0 && couple_J_into_T) { - cout<>> CRASH <<< DUE TO NEGATIVE J, crash_imin/imax = "< 0 && couple_J_into_T) + { + cout << endl + << ">>> CRASH <<< DUE TO NEGATIVE J, crash_imin/imax = " << crash_J_imin << "/" << crash_J_imax << " sample J = " << crashed_meanintensity << " num of crashed cells/total cells = " << crash_J_numcells << "/" << num_cells << " crashed band number = " << crashed_J - 1 << endl; + cout << " @t/dt = " << crashtime << "/" << dt << " stepnum " << steps << endl; + cout << "Writing crash dump into last output and exiting program." << endl; + } } } - } - } - cout<=0; i--) { - checksum = (species[s].de_e[i] > checksum )?species[s].de_e[i]:checksum ; + for (int s = 0; s < num_species; s++) + { + for (int i = num_cells; i >= 0; i--) + { + checksum = (species[s].de_e[i] > checksum) ? species[s].de_e[i] : checksum; } } - - cout<<"Finished at time="<& u_in, std::vector& dudt) { - - - // - // Step 1: Boundary values - // - // - - apply_boundary_left(u_in) ; - apply_boundary_right(u_in) ; - - if(USE_WAVE==1) - add_wave(u_in, base->globalTime); - - if(debug >= 2) - cout<<"Done."<debug > 0) { //Can put in custom debug info here. We left some example code in just for referencing the code usage. - - if(prim[21].temperature < 0. ) { - cout<<" T<0 in execute, steps ="<steps<< " beforeP, T = "<= 2) - cout<<"Done. Starting fluxes."<= 2) - cout<<"Done. Starting sources."<source_pressure_prefactor_left[j] * prim_l[j].pres - - base->source_pressure_prefactor_right[j] * prim_r[j].pres) ,0); +void c_Species::execute(std::vector &u_in, std::vector &dudt) +{ + + // + // Step 1: Boundary values + // + // + + apply_boundary_left(u_in); + apply_boundary_right(u_in); + + if (USE_WAVE == 1) + add_wave(u_in, base->globalTime); + + if (debug >= 2) + cout << "Done." << endl + << " Before compute pressure in species " << speciesname << "... "; + + if (base->debug > 0) + { // Can put in custom debug info here. We left some example code in just for referencing the code usage. + + if (prim[21].temperature < 0.) + { + cout << " T<0 in execute, steps =" << base->steps << " beforeP, T = " << prim[20].temperature << " " << prim[21].temperature << " " << prim[22].temperature << " p_nominal = " << (u[21].u3 - 0.5 * std::pow(u[21].u2, 2.) / u[21].u1) << endl; } - - // - // Step 3: Add it all up to update the conserved variables - // - - for(int j=1; j<=num_cells; j++) { - dudt[j] = (flux[j-1] * base->surf[j-1] - flux[j] * base->surf[j]) / base->vol[j] + (source[j] + source_pressure[j]) ; - - if( debug > 3) { //Or put in your own conditions - char alpha; - cout<<"Debuggin fluxes in cell i= "<steps<= 2) + cout << "Done. Starting fluxes." << endl; + + // + // Step 2: Compute fluxes and sources + // + if (not is_dust_like) + { + for (int j = 0; j <= num_cells; j++) + flux[j] = hllc_flux(j); + } + else + { + for (int j = 0; j <= num_cells; j++) + flux[j] = dust_flux(j); + } + + if (debug >= 2) + cout << "Done. Starting sources." << endl; + + for (int j = 1; j <= num_cells; j++) + { + source[j] = source_grav(u_in[j], j); + source_pressure[j] = AOS(0, -(base->source_pressure_prefactor_left[j] * prim_l[j].pres - base->source_pressure_prefactor_right[j] * prim_r[j].pres), 0); + } + + // + // Step 3: Add it all up to update the conserved variables + // + + for (int j = 1; j <= num_cells; j++) + { + dudt[j] = (flux[j - 1] * base->surf[j - 1] - flux[j] * base->surf[j]) / base->vol[j] + (source[j] + source_pressure[j]); + + if (debug > 3) + { // Or put in your own conditions + char alpha; + cout << "Debuggin fluxes in cell i= " << j << " for species " << speciesname << " at time " << base->steps << endl; + cout << " fl.u1 = " << flux[j - 1].u1 << ": fr.u1 = " << flux[j].u1 << endl; + cout << " fl.u2 = " << flux[j - 1].u2 << ": fr.u2 = " << flux[j].u2 << endl; + cout << " fl.u3 = " << flux[j - 1].u3 << ": fr.u3 = " << flux[j].u3 << endl; + cout << " Cartesian fluxes: Fl-Fr+s = " << ((flux[j - 1].u2 - flux[j].u2) / base->dx[j] + source[j].u2) << endl; + cout << " D_surface/volume=" << (0.5 * (base->surf[j] - base->surf[j - 1]) / base->vol[j]) << " vs. 1/dx=" << (1. / base->dx[j]) << " AL = " << base->surf[j] << " AR = " << base->surf[j - 1] << " VOL = " << base->vol[j] << endl; + cout << endl; + cout << " Reminder, the following three lines must be +- equal for well-balancing:" << endl; + cout << " s = " << source[j].u1 << "/" << source[j].u2 << "/" << source[j].u3 << endl; + cout << " Al*Fl - Ar*Fr + sP = " << ((flux[j - 1].u2 * base->surf[j - 1] - flux[j].u2 * base->surf[j]) / base->vol[j] + source_pressure[j].u2) << endl; + cout << " dP/dr = " << ((prim_l[j].pres - prim_r[j].pres) / base->dx[j]) << endl; + cout << " " << endl; + cout << " sP = " + << "/" << source_pressure[j].u1 << "/" << source_pressure[j].u2 << "/" << source_pressure[j].u3 << endl; + cout << " Al*Fl - Ar*Fr = " << ((flux[j - 1].u2 * base->surf[j - 1] - flux[j].u2 * base->surf[j]) / base->vol[j]) << endl; + cout << " dP/dr + S = " << ((prim_l[j].pres - prim_r[j].pres) / base->dx[j] + source[j].u2) << endl; + cout << endl; + cout << " Al*Fl - Ar*Fr + s = " << ((flux[j - 1].u2 * base->surf[j - 1] - flux[j].u2 * base->surf[j]) / base->vol[j] + (source[j].u2)) << endl; + cout << " u1 : Al*Fl - Ar*Fr + s + sP = " << ((flux[j - 1].u1 * base->surf[j - 1] - flux[j].u1 * base->surf[j]) / base->vol[j] + (source[j].u1 + source_pressure[j].u1)) << endl; + cout << " u2 : Al*Fl - Ar*Fr + s + sP = " << ((flux[j - 1].u2 * base->surf[j - 1] - flux[j].u2 * base->surf[j]) / base->vol[j] + (source[j].u2 + source_pressure[j].u2)) << endl; + cout << " u3 : Al*Fl - Ar*Fr + s + sP = " << ((flux[j - 1].u3 * base->surf[j - 1] - flux[j].u3 * base->surf[j]) / base->vol[j] + (source[j].u3 + source_pressure[j].u3)) << endl; + cout << " " << endl; + cout << " u1 : dudt[1] = " << (dudt[j].u1) << endl; + cout << " u2 : dudt[2] = " << (dudt[j].u2) << endl; + cout << " u3 : dudt[3] = " << (dudt[j].u3) << endl; + + cin >> alpha; + } + } +} /** * The HLLC Riemann solver with pressure-based wave-speed estimate, according to Toro(2007) and references therein. - * - * @param[in] j cell interface number at which to compute the flux. + * + * @param[in] j cell interface number at which to compute the flux. * @return flux at cell interface j */ -AOS c_Species::hllc_flux(int j) +AOS c_Species::hllc_flux(int j) { - int jleft = j, jright = j+1; + int jleft = j, jright = j + 1; AOS flux; int option = 0; - - //Speed of gas - double ul = prim_r[jleft].speed; - double ur = prim_l[jright].speed; - - double pl = prim_r[jleft].pres; + + // Speed of gas + double ul = prim_r[jleft].speed; + double ur = prim_l[jright].speed; + + double pl = prim_r[jleft].pres; double pr = prim_l[jright].pres; double pl_e = pl; double pr_e = pr; - if(base->use_total_pressure) { - pl = base->total_press_r[jleft]; - pr = base->total_press_l[jright]; + if (base->use_total_pressure) + { + pl = base->total_press_r[jleft]; + pr = base->total_press_l[jright]; } - - double dl = prim_r[jleft].density; + + double dl = prim_r[jleft].density; double dr = prim_l[jright].density; - double mom_l = dl*ul ; - double mom_r = dr*ur ; + double mom_l = dl * ul; + double mom_r = dr * ur; - double El = dl*prim_r[jleft].internal_energy + 0.5*mom_l*ul ; - double Er = dr*prim_l[jright].internal_energy + 0.5*mom_r*ur ; + double El = dl * prim_r[jleft].internal_energy + 0.5 * mom_l * ul; + double Er = dr * prim_l[jright].internal_energy + 0.5 * mom_r * ur; - if( (debug > 2) && (j>0 && j<13) ) { - cout<<" IN HLLC, j = "< 2) && (j > 0 && j < 13)) + { + cout << " IN HLLC, j = " << j << " dl/dr = " << dl << "/" << dr << " ul/ul = " << ul << "/" << ur << " pl/pr = " << pl << "/" << pr << " El/Er = " << El << "/" << Er; } - - //Speed of shocks - double SL = ul - prim_r[jleft].sound_speed ; - double SR = ur + prim_l[jright].sound_speed ; - - //Intermediate values in the star region, equations 10.30 -10.39 in Toro - double SS = ( pr-pl+mom_l*(SL - ul)-mom_r*(SR-ur) )/(dl*(SL - ul)-dr*(SR-ur) ); - - if ((SL <= 0) && (SS >= 0)) { - AOS FL = AOS (mom_l, mom_l * ul + pl, ul * (El + pl_e) ); - double comp3_L = El/dl + (SS-ul)*(SS + pl_e/(dl*(SL-ul))); - AOS US_L = AOS(1,SS, comp3_L) * dl * (SL - ul)/(SL-SS); - AOS FS_L = FL + (US_L - AOS(dl, mom_l, El)) * SL; - + + // Speed of shocks + double SL = ul - prim_r[jleft].sound_speed; + double SR = ur + prim_l[jright].sound_speed; + + // Intermediate values in the star region, equations 10.30 -10.39 in Toro + double SS = (pr - pl + mom_l * (SL - ul) - mom_r * (SR - ur)) / (dl * (SL - ul) - dr * (SR - ur)); + + if ((SL <= 0) && (SS >= 0)) + { + AOS FL = AOS(mom_l, mom_l * ul + pl, ul * (El + pl_e)); + double comp3_L = El / dl + (SS - ul) * (SS + pl_e / (dl * (SL - ul))); + AOS US_L = AOS(1, SS, comp3_L) * dl * (SL - ul) / (SL - SS); + AOS FS_L = FL + (US_L - AOS(dl, mom_l, El)) * SL; + flux = FS_L; option = 1; } - else if ((SS <= 0) && (SR >= 0)) { - - AOS FR = AOS (mom_r, mom_r * ur + pr, ur * (Er + pr_e) ); - double comp3_R = Er/dr + (SS-ur)*(SS + pr_e/(dr*(SR-ur))); - AOS US_R = AOS(1,SS, comp3_R) * dr * (SR - ur)/(SR-SS); - AOS FS_R = FR + (US_R - AOS(dr, mom_r, Er)) * SR; - + else if ((SS <= 0) && (SR >= 0)) + { + + AOS FR = AOS(mom_r, mom_r * ur + pr, ur * (Er + pr_e)); + double comp3_R = Er / dr + (SS - ur) * (SS + pr_e / (dr * (SR - ur))); + AOS US_R = AOS(1, SS, comp3_R) * dr * (SR - ur) / (SR - SS); + AOS FS_R = FR + (US_R - AOS(dr, mom_r, Er)) * SR; + flux = FS_R; option = 2; } - else if (SL >= 0) { - AOS FL = AOS (mom_l, mom_l * ul + pl, ul * (El + pl_e) ); + else if (SL >= 0) + { + AOS FL = AOS(mom_l, mom_l * ul + pl, ul * (El + pl_e)); flux = FL; option = 3; } - else if (SR <= 0) { - AOS FR = AOS(mom_r, mom_r * ur + pr, ur * (Er + pr_e) ); + else if (SR <= 0) + { + AOS FR = AOS(mom_r, mom_r * ur + pr, ur * (Er + pr_e)); flux = FR; - option= 4 ; + option = 4; } - if((debug > 2) && (j>0 && j<13)) - cout<<" opt = "< Date: Tue, 4 Jul 2023 16:05:10 +0100 Subject: [PATCH 04/10] got collisional heating with radiation working --- radiation.cpp | 21 +++------------------ test_files/collheat_2spc.par | 1 - test_files/collheat_2spc.spc | 4 ++-- test_files/collheat_2spc_new.spc | 7 ------- test_files/test_conservation.py | 2 +- 5 files changed, 6 insertions(+), 29 deletions(-) delete mode 100644 test_files/collheat_2spc_new.spc diff --git a/radiation.cpp b/radiation.cpp index 57a5582..4e0a6d6 100644 --- a/radiation.cpp +++ b/radiation.cpp @@ -289,7 +289,6 @@ void c_Sim::update_dS_jb(int j, int b) { // for(int s=0; s 1) { @@ -562,7 +559,7 @@ void c_Sim::update_fluxes_FLD() { d[idx_s ] = 1 / dt ; r[idx_rs] = Ts / dt ; r[idx_rs] += (species[s].dS(j) + species[s].dG(j)) / species[s].u[j].u1 / species[s].cv; //Misc heating terms that are not directly related to self-consistent temperature are just added to the rhs - + for(int b=0; b 1) { - - double tau = total_opacity(j,0) * (x_i12[j+1]-x_i12[j]); - + //Heuristic fix for boundary heat bug at large optical depth fill_alpha_basis_arrays(j); compute_collisional_heat_exchange_matrix(j); @@ -597,20 +592,10 @@ void c_Sim::update_fluxes_FLD() { for (int si=0; si < num_species; si++) { int idx = j*stride + (si + num_bands_out) * num_vars + num_bands_out; for (int sj=0; sj < num_species; sj++) { - - if (tau < 1e3) { - d[idx+sj] -= friction_coefficients(si,sj) ; - - } else - d[idx+sj] -= friction_coefficients(si,sj) * 1.e-4 ; - //cout<<" si/sj = "< Date: Tue, 4 Jul 2023 16:49:45 +0100 Subject: [PATCH 05/10] initial attempt at github actions --- .github/ISSUE_TEMPLATE/bug_report.md | 27 ++++++++++++++++ .github/ISSUE_TEMPLATE/feature_request.md | 20 ++++++++++++ .github/workflows/tests.yml | 39 +++++++++++++++++++++++ makefile | 4 ++- test_files/test_coll_heating.py | 13 ++++++-- test_files/test_conservation.py | 15 ++++++--- test_files/test_drag.py | 10 +++++- test_files/test_dustyshock.py | 13 ++++++-- test_files/test_dustywave.py | 20 ++++++++---- test_files/test_irradiation.py | 17 ++++++++-- test_files/test_shock_tube.py | 26 ++++++++++----- test_files/test_steady_state.py | 15 ++++++--- 12 files changed, 187 insertions(+), 32 deletions(-) create mode 100644 .github/ISSUE_TEMPLATE/bug_report.md create mode 100644 .github/ISSUE_TEMPLATE/feature_request.md create mode 100644 .github/workflows/tests.yml diff --git a/.github/ISSUE_TEMPLATE/bug_report.md b/.github/ISSUE_TEMPLATE/bug_report.md new file mode 100644 index 0000000..0857b1e --- /dev/null +++ b/.github/ISSUE_TEMPLATE/bug_report.md @@ -0,0 +1,27 @@ +--- +name: Bug report +about: Create a report to help us improve +title: '' +labels: bug +assignees: '' + +--- + +**Describe the bug -** +1. A clear and concise description of what the bug is. +2. Include the full Traceback of the error. + +**To reproduce -** +1. Include a simple, standalone code snippet that demonstrates the issue. +2. If relevant, provide the input parameters for the frank fit (or attach the .json parameter file). + +**Expected behavior -** +A clear and concise description of what you expected to happen. + +**Your setup -** + - Code version + - Compiler + - Operating system + +**Additional context -** +Add any other context about the problem (e.g., any peculiarity of the data or of the problem you're trying to solve?). diff --git a/.github/ISSUE_TEMPLATE/feature_request.md b/.github/ISSUE_TEMPLATE/feature_request.md new file mode 100644 index 0000000..8170920 --- /dev/null +++ b/.github/ISSUE_TEMPLATE/feature_request.md @@ -0,0 +1,20 @@ +--- +name: Feature request +about: Suggest an idea for frank +title: '' +labels: enhancement +assignees: '' + +--- + +**Is your feature request related to a problem? Please describe -** +A clear and concise description of what the problem is. + +**Describe the solution you'd like -** +A description of what you want to happen. + +**Describe alternatives you've considered -** +A description of any alternative solutions or features you've considered. + +**Additional context -** +Add any other context or screenshots about the feature request. diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml new file mode 100644 index 0000000..83cc9e8 --- /dev/null +++ b/.github/workflows/tests.yml @@ -0,0 +1,39 @@ +name: tests + +# Run this workflow for every commit +on: push + +jobs: + test: + name: Run tests + runs-on: ubuntu-latest + + steps: + - name: Checkout code + uses: actions/checkout@v3 + + - name: Set up Python 3.11 + uses: actions/setup-python@v2 + with: + python-version: 3.11 + + - name: Install test dependencies + run: | + pip install numpy, scipy, matplotlib, pytest + + - name: Install code dependencies + run: | + sudo apt-get install -y libgsl-dev libeigen3-dev + + - name: Build tests + run: make tests + + - name: Run the code tests + run: | + ./tests > /dev/null + + - name: Check the results + run: | + cd test_files + mkdir -p test-reports + py.test -v --junitxml=test-reports/junit.xml . diff --git a/makefile b/makefile index 32b98ee..c31ece3 100644 --- a/makefile +++ b/makefile @@ -41,10 +41,12 @@ aiolos: $(OBJ) makefile %.o: %.cpp makefile aiolos.h $(CXX) $(CPPFLAGS) $(CXXFLAGS) -c $< $(BFLAGS) -o $@ -.PHONY: tests +.PHONY: run_tests tests: $(TEST_OBJ) makefile aiolos.h $(CXX) -o $@ $(TEST_OBJ) $(CXXFLAGS) $(LDFLAGS) + +run_tests: tests ./tests > /dev/null cd test_files ; python3 test_shock_tube.py -p cd test_files ; python3 test_steady_state.py diff --git a/test_files/test_coll_heating.py b/test_files/test_coll_heating.py index 3c47907..b1ea3b9 100644 --- a/test_files/test_coll_heating.py +++ b/test_files/test_coll_heating.py @@ -21,7 +21,7 @@ def _solution_2sp(aij_dt, cv1, cv2, T1, T2): return T1, T2 -def check_collisional_heating_2spc(problem, L1_target=None): +def check_collisional_heating_2spc(problem, L1_target=None, raise_error=False): param_file = problem + '.par' params = load_aiolos_params(param_file) @@ -66,12 +66,19 @@ def check_collisional_heating_2spc(problem, L1_target=None): else: print('Test {} L1 checked failed:'.format(problem)) print('\tL1={}, target={}'.format(L1, L1_target)) + + np.testing.assert_array_less(L1, L1_target) + else: print('Test {} L1 values:'.format(problem)) print('\tL1={}'.format(L1)) +def test_coll_heating(): + check_collisional_heating_2spc('collheat_2spc', [1.5e-3, 2.5e-3]) if __name__ == "__main__": - check_collisional_heating_2spc('collheat_2spc', [1.5e-3, 2.5e-3]) - + try: + test_coll_heating() + except AssertionError: + pass diff --git a/test_files/test_conservation.py b/test_files/test_conservation.py index c6e5948..374ae60 100644 --- a/test_files/test_conservation.py +++ b/test_files/test_conservation.py @@ -70,15 +70,15 @@ def check_conservation(problem, species='hydrogen.spc', L1_target=None): else: print('Conservation test {} L1 checked failed:'.format(problem)) print('\tL1={}, target={}'.format(L1, L1_target)) + + np.testing.assert_array_less(L1, L1_target) + else: print('Conservation test {} L1 values:'.format(problem)) print('\tL1={}'.format(L1)) -if __name__ == "__main__": - - - +def test_conservation(): check_conservation("soundwave_128", L1_target=[5e-16, 2e-11,3e-16]) check_conservation("dustywave_stiff", L1_target=[5e-16, 5e-8, 2e-14]) @@ -97,3 +97,10 @@ def check_conservation(problem, species='hydrogen.spc', L1_target=None): check_conservation("collheat_2spc_rad", "collheat_2spc.spc", L1_target=[4e-16, 1e-18, 2e-5]) + + +if __name__ == "__main__": + try: + test_conservation() + except AssertionError: + pass \ No newline at end of file diff --git a/test_files/test_drag.py b/test_files/test_drag.py index 1e4521f..4e876f8 100644 --- a/test_files/test_drag.py +++ b/test_files/test_drag.py @@ -29,13 +29,21 @@ def check_drag(problem, spc, L1_target=None): else: print('Test {} L1 checked failed:'.format(problem)) print('\tL1={}, target={}'.format(L1, L1_target)) + + np.testing.assert_array_less(L1, L1_target) + else: print('Test {} L1 values:'.format(problem)) print('\tL1={}'.format(L1)) -if __name__ == "__main__": +def test_drag(): check_drag("friction_2spc", 'H0', [3e-16, 3e-14, 2e-14]) check_drag("friction_2spc_phys", 'H0', [2e-30, 2e-29, 2e-20]) +if __name__ == "__main__": + try: + test_drag() + except AssertionError: + pass diff --git a/test_files/test_dustyshock.py b/test_files/test_dustyshock.py index 2058174..f36273e 100644 --- a/test_files/test_dustyshock.py +++ b/test_files/test_dustyshock.py @@ -116,10 +116,16 @@ def check_dusty_shock(L1s=None): else: print('Test {} L1 checked failed:'.format(name)) print('\tL1={}, target={}'.format(L1, L1s)) + + np.testing.assert_array_less(L1, L1s) + else: print('Test {} L1 values:'.format(name)) print('\tL1={}'.format(L1)) +def test_dusty_shock(): + check_dusty_shock([0.541, 0.540]) + if __name__ == "__main__": import argparse @@ -130,8 +136,11 @@ def check_dusty_shock(L1s=None): help="Make plots of the results") args = parser.parse_args() - check_dusty_shock([0.541, 0.540]) - + try: + test_dusty_shock() + except AssertionError: + pass + if args.make_plots: f = plot_dusty_shock() #f.show() diff --git a/test_files/test_dustywave.py b/test_files/test_dustywave.py index 3ed1cdb..9076a2c 100644 --- a/test_files/test_dustywave.py +++ b/test_files/test_dustywave.py @@ -175,6 +175,9 @@ def check_dustywave_problem(name, L1s=None, make_plots=False): else: print('Test {} L1 checked failed:'.format(name)) print('\tL1={}, target={}'.format(L1, L1s)) + + np.testing.assert_array_less(L1, L1s) + else: print('Test {} L1 values:'.format(name)) print('\tL1={}'.format(L1)) @@ -182,6 +185,13 @@ def check_dustywave_problem(name, L1s=None, make_plots=False): if make_plots: plot_dustywave(name) +def test_dustywave(make_plots=False): + + check_dustywave_problem("dustywave_nonstiff", [1.25e-8, 1.34e-8], + make_plots=make_plots) + check_dustywave_problem("dustywave_stiff", [4.78e-8, 4.78e-8], + make_plots=make_plots) + if __name__ == "__main__": import argparse @@ -192,9 +202,7 @@ def check_dustywave_problem(name, L1s=None, make_plots=False): help="Make plots of the results") args = parser.parse_args() - - - check_dustywave_problem("dustywave_nonstiff", [1.25e-8, 1.34e-8], - make_plots=args.make_plots) - check_dustywave_problem("dustywave_stiff", [4.78e-8, 4.78e-8], - make_plots=args.make_plots) + try: + test_dustywave (args.make_plots) + except AssertionError: + pass diff --git a/test_files/test_irradiation.py b/test_files/test_irradiation.py index 6060473..7b338e6 100644 --- a/test_files/test_irradiation.py +++ b/test_files/test_irradiation.py @@ -20,7 +20,7 @@ def Guillot_2band(tau, Tint4, Tirr40, gamma0, Tirr41, gamma1, return (Tint4*term_int + Tirr40*term0 + Tirr41*term1)**0.25 -def test_structure(sim, Tint, L1_target, make_plots): +def check_structure(sim, Tint, L1_target, make_plots): filename = 'diagnostic_' + sim + '_t-1.dat' data = load_aiolos_diag(filename) @@ -44,6 +44,8 @@ def test_structure(sim, Tint, L1_target, make_plots): else: print("Irradiation test L1 check failed. " + "L1={}, expected={}".format(L1,L1_target)) + + if make_plots: import matplotlib.pyplot as plt @@ -63,6 +65,13 @@ def test_structure(sim, Tint, L1_target, make_plots): os.makedirs('plots', exist_ok=True) plt.savefig('plots/Irradiation.png') + np.testing.assert_array_less(L1, L1_target) + + +def test_structure(make_plots=False): + check_structure('irradiation', 175.0, 1.5e-3, make_plots) + + if __name__ == "__main__": import argparse @@ -72,5 +81,7 @@ def test_structure(sim, Tint, L1_target, make_plots): help="Make plots of the results") args = parser.parse_args() - test_structure('irradiation', 175.0, 1.5e-3, args.make_plots) - + try: + test_structure(args.make_plots) + except AssertionError: + pass diff --git a/test_files/test_shock_tube.py b/test_files/test_shock_tube.py index 13e6271..47037c4 100644 --- a/test_files/test_shock_tube.py +++ b/test_files/test_shock_tube.py @@ -265,7 +265,8 @@ def plot_riemann_solution(par_num): hspace=0.05, wspace=0.235) return f -def check_riemann_problem(par_num, L1_target=0, make_plots=False): +def check_riemann_problem(par_num, L1_target=0, + make_plots=False, raise_error=False): param_file = 'shock_tube{}.par'.format(par_num) snap_file = 'output_shock_tube{}_H2_t-1.dat'.format(par_num) @@ -280,6 +281,11 @@ def check_riemann_problem(par_num, L1_target=0, make_plots=False): L1 = np.abs(data['density'] - exact.density(data['x']-x0, t)).mean() + + if raise_error: + np.testing.assert_array_less(L1, L1_target) + + if L1 <= L1_target: print('Riemann test {} L1 check passed'.format(par_num)) else: @@ -291,7 +297,13 @@ def check_riemann_problem(par_num, L1_target=0, make_plots=False): fig_name = os.path.join('plots', 'shock_tube{}.png'.format(par_num)) plot_riemann_solution(par_num).savefig(fig_name) - + +def test_shock_tube(make_plots=False): + L1_errors = [np.nan, 0.016, 0.033, 0.25, + 1.06, 0.078, 0.0040, 0.0035] + for i in range(1, 8): + check_riemann_problem(i, L1_errors[i], make_plots, not make_plots) + if __name__ == "__main__": import argparse @@ -301,9 +313,7 @@ def check_riemann_problem(par_num, L1_target=0, make_plots=False): help="Make plots of the results") args = parser.parse_args() - - L1_errors = [np.nan, 0.016, 0.033, 0.25, - 1.06, 0.078, 0.0040, 0.0035] - for i in range(1, 8): - check_riemann_problem(i, L1_errors[i], args.make_plots) - + try: + test_shock_tube(args.make_plots) + except AssertionError: + pass diff --git a/test_files/test_steady_state.py b/test_files/test_steady_state.py index 8b9203a..99f6a1f 100644 --- a/test_files/test_steady_state.py +++ b/test_files/test_steady_state.py @@ -29,14 +29,16 @@ def check_steady_state(problem, L1_target=None): else: print('Test {} L1 checked failed:'.format(problem)) print('\tL1={}, target={}'.format(L1, L1_target)) + + np.testing.assert_array_less(L1, L1_target) + else: print('Test {} L1 values:'.format(problem)) print('\tL1={}'.format(L1)) - -if __name__ == "__main__": - check_steady_state("shock_tube6_H2", [0,0,0]) - check_steady_state("shock_tube7_H2", [4.9e-4, 4.9e-4, 0]) +def test_steady_states(): + check_steady_state("shock_tube6_H2", [1e-30,1e-30,1e-30]) + check_steady_state("shock_tube7_H2", [4.9e-4, 4.9e-4, 1e-30]) check_steady_state("soundwave_32_H2", [8.72e-8,]*3) check_steady_state("soundwave_64_H2", [2.95e-8,]*3) @@ -47,6 +49,11 @@ def check_steady_state(problem, L1_target=None): check_steady_state("planet_cartesian_H2", [4e-13, 4e-10, 0.2]) check_steady_state("planet_spherical_H2", [1e-11, 2e-10, 2.5]) +if __name__ == "__main__": + try: + test_steady_states() + except AssertionError: + pass From 566f9f3c52fc30e151bde3d757871b92b36b7628 Mon Sep 17 00:00:00 2001 From: Richard Booth Date: Tue, 4 Jul 2023 16:50:58 +0100 Subject: [PATCH 06/10] typos --- .github/workflows/tests.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 83cc9e8..736c7dd 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -19,7 +19,7 @@ jobs: - name: Install test dependencies run: | - pip install numpy, scipy, matplotlib, pytest + pip install numpy scipy matplotlib pytest - name: Install code dependencies run: | From 29b95b5b0d1dcbc37be091c4bcdd7269d05af37c Mon Sep 17 00:00:00 2001 From: Richard Booth Date: Tue, 4 Jul 2023 16:52:21 +0100 Subject: [PATCH 07/10] parallelized test build --- .github/workflows/tests.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 736c7dd..3181dc0 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -26,7 +26,7 @@ jobs: sudo apt-get install -y libgsl-dev libeigen3-dev - name: Build tests - run: make tests + run: make -j tests - name: Run the code tests run: | From ac1aa1bdac4279e09f87baf5ffe2548bf32957dc Mon Sep 17 00:00:00 2001 From: Richard Booth Date: Wed, 5 Jul 2023 10:52:03 +0100 Subject: [PATCH 08/10] removed unnecessary libraries --- .github/workflows/tests.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 3181dc0..982dc9c 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -19,7 +19,7 @@ jobs: - name: Install test dependencies run: | - pip install numpy scipy matplotlib pytest + pip install numpy pytest - name: Install code dependencies run: | From a65382b853dadbcaed4454e1a74e7bf08d81bb62 Mon Sep 17 00:00:00 2001 From: Richard Booth Date: Wed, 5 Jul 2023 13:09:48 +0100 Subject: [PATCH 09/10] Added working radiation tests --- test_files/irradiation1.par | 85 +++++++++++++++++++++ test_files/irradiation2.par | 85 +++++++++++++++++++++ test_files/load_aiolos.py | 76 +++++++++--------- test_files/main.cpp | 56 +++++++------- test_files/static5_euv1e3_C_1s2b_rerun2.par | 85 +++++++++++++++++++++ test_files/test_conservation.py | 4 +- test_files/test_irradiation.py | 59 ++++++++++---- 7 files changed, 369 insertions(+), 81 deletions(-) create mode 100644 test_files/irradiation1.par create mode 100644 test_files/irradiation2.par create mode 100644 test_files/static5_euv1e3_C_1s2b_rerun2.par diff --git a/test_files/irradiation1.par b/test_files/irradiation1.par new file mode 100644 index 0000000..7d087a6 --- /dev/null +++ b/test_files/irradiation1.par @@ -0,0 +1,85 @@ +# +# Parameterfile for atmosphere simulations +# Do not use tabs, only white spaces as value separators +# +SPECIES_FILE irradiation.spc +DENSITY_FLOOR 1e-20 +PARI_TIME_TMAX 1.01e11 +CFL_BREAK_TIME 1e3 +PARI_TIME_OUTPUT 2.e10 +PARI_CFLFACTOR 0.1 +PARI_TIME_DT 1e10 +PARI_DPHI_FACTOR 1.5e-2 +ENERGY_EPSILON 0.01 +MAX_TIMESTEP_CHANGE 1.01 +XI_RAD 2. +# Control parameters +# +# Problem numbers explained: shock_tube, no grav == 1, planet with grav == 2# +# Friction solver: analytic == 0, numerical == 1 +# Collision model: constant == C, physical == P +PARI_PROBLEM_NUMBER 2 +PARI_NUM_SPECIES 1 +PARI_NUM_BANDS 2 +NUM_BANDS_OUT 1 +DO_HYDRO 0 +#CONST_TEMP 1430. +PARI_CONST_TEMP 2.7 +INIT_TEMPERATURE_MODEL P +INIT_J_FACTOR 1e0 +INIT_T_TEMP 2.7 +USE_PLANET_TEMPERATURE 1 +PARI_TPLANET 350. +PARI_CORE_CV 1e9 + +# Solar radiation temps at 0.1 AU: 645K for bond albedo=0.75, 913K for bond albedo = 0 +# print 5777*(1./20./2.)**0.5*(1-0.75)**0.25 + +FRICTION_SOLVER 1 +PARI_COLL_MODEL P +PARI_ALPHA_COLL 1. +PARI_ORDER 2 +PARI_USE_RADIATION 1 +PARI_OPACITY_MODEL C +PRESSURE_BROADENING 0. +CONSTOPA_SOLAR_FACTOR 0.35 +PARI_CONST_OPAC 1.0 + +PARI_SELF_GRAV_SWITCH 0 +PARI_LINEAR_GRAV 0 +PARI_DEBUGLEVEL 1 +PARI_INIT_STATIC 1 +USE_TIDES 1 +#Physical parameters + +PARI_PLANET_MASS 224. +PARI_INIT_DATA_SWITCH 1 +PARI_INIT_DATA_U1 1e-5 +PARI_INIT_DATA_U2 1.0 +PARI_INIT_DATA_U3 1e-68 + +#Radiation parameters +PARI_LAM_MIN 9.0e-2 +PARI_LAM_MAX 9.0e-2 +PARI_LAM_PER_DECADE 0.1 +PARI_TSTAR 6070. +PARI_RSTAR 0.1227 +PARI_UVSTAR 1.e1 + +PARI_PLANET_DIST 0.05 + +#Numerical parameters +PARI_BOUND_TYPE 2 +PARI_BOUND_TYPE_LEFT 2 +PARI_BOUND_TYPE_RIGHT 1 +PARI_DOMAIN_DX 1e-3 +PARI_GRID_TYPE 1 +PARI_GEOMETRY 2 +PARI_CELLS_PER_DECADE 1000 +#Domain min/max are 8 me core radius for earth density, and 100*hill radius +PARI_DOMAIN_MIN 9.200e9 +PARI_DOMAIN_MAX 9.2e11 + +#Scenarios +PARI_INIT_WIND 1 +REVERSE_HYDROSTAT_CONSTRUCTION 1 diff --git a/test_files/irradiation2.par b/test_files/irradiation2.par new file mode 100644 index 0000000..abdcc84 --- /dev/null +++ b/test_files/irradiation2.par @@ -0,0 +1,85 @@ +# +# Parameterfile for atmosphere simulations +# Do not use tabs, only white spaces as value separators +# +SPECIES_FILE mix_static_paperplot.spc +DENSITY_FLOOR 1e-20 +PARI_TIME_TMAX 1.01e13 +CFL_BREAK_TIME 1e3 +PARI_TIME_OUTPUT 2.e13 +PARI_CFLFACTOR 0.1 +PARI_TIME_DT 1e10 +PARI_DPHI_FACTOR 1.5e-2 +ENERGY_EPSILON 0.01 +MAX_TIMESTEP_CHANGE 1.01 +XI_RAD 2. +# Control parameters +# +# Problem numbers explained: shock_tube, no grav == 1, planet with grav == 2# +# Friction solver: analytic == 0, numerical == 1 +# Collision model: constant == C, physical == P +PARI_PROBLEM_NUMBER 2 +PARI_NUM_SPECIES 1 +PARI_NUM_BANDS 2 +NUM_BANDS_OUT 1 +DO_HYDRO 0 +#CONST_TEMP 1430. +PARI_CONST_TEMP 2.7 +INIT_TEMPERATURE_MODEL P +INIT_J_FACTOR 1e0 +INIT_T_TEMP 2.7 +USE_PLANET_TEMPERATURE 1 +PARI_TPLANET 350. +PARI_CORE_CV 1e9 + +# Solar radiation temps at 0.1 AU: 645K for bond albedo=0.75, 913K for bond albedo = 0 +# print 5777*(1./20./2.)**0.5*(1-0.75)**0.25 + +FRICTION_SOLVER 1 +PARI_COLL_MODEL P +PARI_ALPHA_COLL 1. +PARI_ORDER 2 +PARI_USE_RADIATION 1 +PARI_OPACITY_MODEL C +PRESSURE_BROADENING 0. +CONSTOPA_SOLAR_FACTOR 0.35 +PARI_CONST_OPAC 1.0 + +PARI_SELF_GRAV_SWITCH 0 +PARI_LINEAR_GRAV 0 +PARI_DEBUGLEVEL 1 +PARI_INIT_STATIC 1 +USE_TIDES 1 +#Physical parameters + +PARI_PLANET_MASS 224. +PARI_INIT_DATA_SWITCH 1 +PARI_INIT_DATA_U1 1e-5 +PARI_INIT_DATA_U2 1.0 +PARI_INIT_DATA_U3 1e-68 + +#Radiation parameters +PARI_LAM_MIN 9.0e-2 +PARI_LAM_MAX 9.0e-2 +PARI_LAM_PER_DECADE 0.1 +PARI_TSTAR 6070. +PARI_RSTAR 0. +PARI_UVSTAR 0.e1 + +PARI_PLANET_DIST 0.05 + +#Numerical parameters +PARI_BOUND_TYPE 2 +PARI_BOUND_TYPE_LEFT 2 +PARI_BOUND_TYPE_RIGHT 1 +PARI_DOMAIN_DX 1e-3 +PARI_GRID_TYPE 1 +PARI_GEOMETRY 2 +PARI_CELLS_PER_DECADE 1000 +#Domain min/max are 8 me core radius for earth density, and 100*hill radius +PARI_DOMAIN_MIN 9.200e9 +PARI_DOMAIN_MAX 9.2e11 + +#Scenarios +PARI_INIT_WIND 1 +REVERSE_HYDROSTAT_CONSTRUCTION 1 diff --git a/test_files/load_aiolos.py b/test_files/load_aiolos.py index d198b99..6117d99 100644 --- a/test_files/load_aiolos.py +++ b/test_files/load_aiolos.py @@ -8,14 +8,11 @@ def load_aiolos_snap(filename): ('pressure', 'f8'), ('velocity', 'f8'), ('temperature', 'f8'), - ('dt_cs', 'f8'), - ('dt_v', 'f8'), - ('soundspeed', 'f8'), - ('dt_E', 'f8'), - - - ] - cols=[0,1,2,3,10,11,12,13,14,15,16] + ('soundspeed', 'f8'), + ('potential', 'f8'), + ('cooling', 'f8'), + ('heating', 'f8')] + cols=[0,1,2,3,10,11,12,15,19,21,22] data = np.genfromtxt(filename, dtype=dtype, usecols=cols) return data @@ -34,14 +31,14 @@ def load_aiolos_params(filename): params[key.strip()] = val.strip() return params -def load_aiolos_diag(filename, bands=None, species=None): - if bands is None or species is None: +def load_aiolos_diag(filename, bands_in=None, bands_out=1, species=None): + if bands_in is None or species is None: # Work out the number of bands: with open(filename) as f: line = f.readline() cols = len(line.strip().split()) - def c(b, s): return 5 + 7*b + 3*b*s + s + def c(b, s): return 5 + (3 + s)*b + s + (4+2*s)*bands_out res = [] b = 1 while True: @@ -58,42 +55,45 @@ def c(b, s): return 5 + 7*b + 3*b*s + s raise ValueError("Error when finding bands and species") if len(res) > 1: for b, s in res: - if b == bands or s == species: - bands = b + if b == bands_in or s == species: + bands_in = b species = s break else: raise ValueError("Error determining the number of bands and" "species, help out by specifying them") else: - bands, species = res[0] + bands_in, species = res[0] # Construct the dtype: - def add_cols(name, dt='f8'): + def add_cols(name, bands=bands_out, dt='f8'): return [ (name + str(i), dt) for i in range(bands)] def add_kappa_cols(): cols = [] - for i in range(bands): + for i in range(bands_in): + for j in range(species): + cols += [('kappa_P(Tsun)' + str(i) + '_' + str(j), 'f8')] + for i in range(bands_out): for j in range(species): cols += [('kappa_R' + str(i) + '_' + str(j), 'f8')] + for j in range(species): cols += [('kappa_P' + str(i) + '_' + str(j), 'f8')] - cols += [('kappa_P(Tsun)' + str(i) + '_' + str(j), 'f8')] return cols dtype = [('x', 'f8'), ('Jtot', 'f8'), ('Stot', 'f8')] - dtype += add_cols('J') - dtype += add_cols('S') - dtype += add_cols('dS') - dtype += add_cols('rho_kappaR') - dtype += add_cols('tau_cell') - dtype += add_cols('tau_radial') + dtype += add_cols('J',bands_out) + dtype += add_cols('S',bands_in) + dtype += add_cols('dS', bands_in) + dtype += add_cols('rho_kappaR',bands_out) + dtype += add_cols('tau_cell',bands_out) + dtype += add_cols('tau_radial', bands_in) dtype += add_kappa_cols() dtype += [('T_rad', 'f8')] dtype += [ ('T_gas' + str(i), 'f8') for i in range(species)] dtype += [('pressure' , 'f8')] - dtype += add_cols('Flux') + dtype += add_cols('Flux',bands_out) # Load the data return np.genfromtxt(filename, dtype=dtype) @@ -125,9 +125,11 @@ def load_aiolos_species(filename): import sys import matplotlib.pyplot as plt - f, ax = plt.subplots(3,1) + f, ax = plt.subplots(4,1, sharex=True) - for snap in sys.argv[1:]: + snaps = [s for s in sys.argv[1:] if s[0] == 'o'] + [s for s in sys.argv[1:] if s[0] != 'o'] + + for snap in snaps: if snap.startswith('output'): data = load_aiolos_snap(snap) if 'Gas' in snap: @@ -141,26 +143,30 @@ def load_aiolos_species(filename): ax[1].semilogx(data['x'], data['temperature'], ls='', marker=m, label=snap) + #ax[1].semilogx(data['x'], data['pressure']/data['density']**1.29, ls='', + # marker=m, label=snap) - #ax[2].semilogx(data['x'], data['dt_cs'], ls='-', - # marker=m, label=snap) - #ax[2].semilogx(data['x'], data['dt_v'], ls='--', - # marker=m, label=snap) - #ax[2].semilogx(data['x'], data['dt_E'], ls=':', - # marker=m, label=snap) ax[2].semilogx(data['x'], data['velocity'], ls='', marker=m, label=snap) - ax[2].semilogx(data['x'], data['soundspeed'], ls='--', - marker=m) + + l,= ax[3].loglog(data['x'], data['heating'], ls='-', + label=snap) + ax[3].semilogx(data['x'], np.abs(data['cooling']), ls='--', + c=l.get_color()) + #ax[2].semilogx(data['x'], data['soundspeed'], ls='--', + # marker=m) elif snap.startswith('diagnostic'): - data = load_aiolos_diag(snap) + data = load_aiolos_diag(snap, species=2) ax[1].semilogx(data['x'], data['T_rad'], ls='', marker='^', label=snap) + #ax[3].semilogx(data['x'], np.abs(data['dS0']), marker='^', ls='') ax[-1].set_xlabel('r') ax[0].set_ylabel('density') ax[1].set_ylabel('temperature') + ax[1].set_ylim(0, 4000) ax[2].set_ylabel('velocity') + ax[3].set_ylabel('heating / cooling') ax[1].legend(ncol=3) diff --git a/test_files/main.cpp b/test_files/main.cpp index ddb6751..2873541 100644 --- a/test_files/main.cpp +++ b/test_files/main.cpp @@ -37,42 +37,42 @@ bool run_test(std::string test_name, std::string spc = "hydrogen.spc", int main() { - int count = 0 ; - count += !run_test("shock_tube1.par") ; - count += !run_test("shock_tube2.par") ; - count += !run_test("shock_tube3.par") ; - count += !run_test("shock_tube4.par") ; - count += !run_test("shock_tube5.par") ; - count += !run_test("shock_tube6.par") ; - count += !run_test("shock_tube7.par") ; + int fail_count = 0 ; + fail_count += !run_test("shock_tube1.par") ; + fail_count += !run_test("shock_tube2.par") ; + fail_count += !run_test("shock_tube3.par") ; + fail_count += !run_test("shock_tube4.par") ; + fail_count += !run_test("shock_tube5.par") ; + fail_count += !run_test("shock_tube6.par") ; + fail_count += !run_test("shock_tube7.par") ; - count += !run_test("soundwave_32.par") ; - count += !run_test("soundwave_64.par") ; - count += !run_test("soundwave_128.par") ; - count += !run_test("soundwave_256.par") ; - count += !run_test("soundwave_512.par") ; + fail_count += !run_test("soundwave_32.par") ; + fail_count += !run_test("soundwave_64.par") ; + fail_count += !run_test("soundwave_128.par") ; + fail_count += !run_test("soundwave_256.par") ; + fail_count += !run_test("soundwave_512.par") ; - count += !run_test("friction_2spc.par", "friction_2spc.spc") ; - count += !run_test("friction_2spc_an.par", "friction_2spc.spc") ; - count += !run_test("friction_2spc_phys.par", "friction_2spc.spc") ; - count += !run_test("friction_2spc_phys_an.par", "friction_2spc.spc") ; - count += !run_test("friction_6spc.par", "friction_6spc.spc") ; + fail_count += !run_test("friction_2spc.par", "friction_2spc.spc") ; + fail_count += !run_test("friction_2spc_an.par", "friction_2spc.spc") ; + fail_count += !run_test("friction_2spc_phys.par", "friction_2spc.spc") ; + fail_count += !run_test("friction_2spc_phys_an.par", "friction_2spc.spc") ; + fail_count += !run_test("friction_6spc.par", "friction_6spc.spc") ; - count += !run_test("collheat_2spc.par", "collheat_2spc.spc") ; - count += !run_test("collheat_2spc_rad.par", "collheat_2spc.spc") ; + fail_count += !run_test("collheat_2spc.par", "collheat_2spc.spc") ; + fail_count += !run_test("collheat_2spc_rad.par", "collheat_2spc.spc") ; - count += !run_test("dustywave_nonstiff.par") ; - count += !run_test("dustywave_stiff.par") ; + fail_count += !run_test("dustywave_nonstiff.par") ; + fail_count += !run_test("dustywave_stiff.par") ; - count += !run_test("dusty_shock.par") ; + fail_count += !run_test("dusty_shock.par") ; - count += !run_test("planet_cartesian.par") ; - count += !run_test("planet_spherical.par") ; - - count += !run_test("irradiation.par", "irradiation.spc") ; + fail_count += !run_test("planet_cartesian.par") ; + fail_count += !run_test("planet_spherical.par") ; + fail_count += !run_test("irradiation1.par", "irradiation.spc") ; + fail_count += !run_test("irradiation2.par", "irradiation.spc") ; - return count ; + return fail_count ; } //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ diff --git a/test_files/static5_euv1e3_C_1s2b_rerun2.par b/test_files/static5_euv1e3_C_1s2b_rerun2.par new file mode 100644 index 0000000..7d087a6 --- /dev/null +++ b/test_files/static5_euv1e3_C_1s2b_rerun2.par @@ -0,0 +1,85 @@ +# +# Parameterfile for atmosphere simulations +# Do not use tabs, only white spaces as value separators +# +SPECIES_FILE irradiation.spc +DENSITY_FLOOR 1e-20 +PARI_TIME_TMAX 1.01e11 +CFL_BREAK_TIME 1e3 +PARI_TIME_OUTPUT 2.e10 +PARI_CFLFACTOR 0.1 +PARI_TIME_DT 1e10 +PARI_DPHI_FACTOR 1.5e-2 +ENERGY_EPSILON 0.01 +MAX_TIMESTEP_CHANGE 1.01 +XI_RAD 2. +# Control parameters +# +# Problem numbers explained: shock_tube, no grav == 1, planet with grav == 2# +# Friction solver: analytic == 0, numerical == 1 +# Collision model: constant == C, physical == P +PARI_PROBLEM_NUMBER 2 +PARI_NUM_SPECIES 1 +PARI_NUM_BANDS 2 +NUM_BANDS_OUT 1 +DO_HYDRO 0 +#CONST_TEMP 1430. +PARI_CONST_TEMP 2.7 +INIT_TEMPERATURE_MODEL P +INIT_J_FACTOR 1e0 +INIT_T_TEMP 2.7 +USE_PLANET_TEMPERATURE 1 +PARI_TPLANET 350. +PARI_CORE_CV 1e9 + +# Solar radiation temps at 0.1 AU: 645K for bond albedo=0.75, 913K for bond albedo = 0 +# print 5777*(1./20./2.)**0.5*(1-0.75)**0.25 + +FRICTION_SOLVER 1 +PARI_COLL_MODEL P +PARI_ALPHA_COLL 1. +PARI_ORDER 2 +PARI_USE_RADIATION 1 +PARI_OPACITY_MODEL C +PRESSURE_BROADENING 0. +CONSTOPA_SOLAR_FACTOR 0.35 +PARI_CONST_OPAC 1.0 + +PARI_SELF_GRAV_SWITCH 0 +PARI_LINEAR_GRAV 0 +PARI_DEBUGLEVEL 1 +PARI_INIT_STATIC 1 +USE_TIDES 1 +#Physical parameters + +PARI_PLANET_MASS 224. +PARI_INIT_DATA_SWITCH 1 +PARI_INIT_DATA_U1 1e-5 +PARI_INIT_DATA_U2 1.0 +PARI_INIT_DATA_U3 1e-68 + +#Radiation parameters +PARI_LAM_MIN 9.0e-2 +PARI_LAM_MAX 9.0e-2 +PARI_LAM_PER_DECADE 0.1 +PARI_TSTAR 6070. +PARI_RSTAR 0.1227 +PARI_UVSTAR 1.e1 + +PARI_PLANET_DIST 0.05 + +#Numerical parameters +PARI_BOUND_TYPE 2 +PARI_BOUND_TYPE_LEFT 2 +PARI_BOUND_TYPE_RIGHT 1 +PARI_DOMAIN_DX 1e-3 +PARI_GRID_TYPE 1 +PARI_GEOMETRY 2 +PARI_CELLS_PER_DECADE 1000 +#Domain min/max are 8 me core radius for earth density, and 100*hill radius +PARI_DOMAIN_MIN 9.200e9 +PARI_DOMAIN_MAX 9.2e11 + +#Scenarios +PARI_INIT_WIND 1 +REVERSE_HYDROSTAT_CONSTRUCTION 1 diff --git a/test_files/test_conservation.py b/test_files/test_conservation.py index 374ae60..e51e95e 100644 --- a/test_files/test_conservation.py +++ b/test_files/test_conservation.py @@ -42,11 +42,11 @@ def check_conservation(problem, species='hydrogen.spc', L1_target=None): c4pi = 2.99792458e10/(4*np.pi) N = int(params.get('NUM_BANDS', 1)) diag = 'diagnostic_' + problem - d0 = load_aiolos_diag(diag + '_t0.dat') + d0 = load_aiolos_diag(diag + '_t0.dat', bands_in=N) for i in range(N): init[2] += d0['J{}'.format(i)][num_bound-1:-num_bound+1].sum()/c4pi - df = load_aiolos_diag(diag + '_t-1.dat') + df = load_aiolos_diag(diag + '_t-1.dat', bands_in=N) for i in range(N): final[2] += df['J{}'.format(i)][num_bound-1:-num_bound+1].sum()/c4pi diff --git a/test_files/test_irradiation.py b/test_files/test_irradiation.py index 7b338e6..155b81e 100644 --- a/test_files/test_irradiation.py +++ b/test_files/test_irradiation.py @@ -1,8 +1,11 @@ import os import numpy as np -from load_aiolos import load_aiolos_snap, load_aiolos_diag - +from load_aiolos import ( + load_aiolos_snap, + load_aiolos_diag, + load_aiolos_params +) sigma_rad = 5.67e-5 def Guillot_2band(tau, Tint4, Tirr40, gamma0, Tirr41, gamma1, @@ -20,29 +23,45 @@ def Guillot_2band(tau, Tint4, Tirr40, gamma0, Tirr41, gamma1, return (Tint4*term_int + Tirr40*term0 + Tirr41*term1)**0.25 -def check_structure(sim, Tint, L1_target, make_plots): +def check_structure(sim, L1_target, make_plots): filename = 'diagnostic_' + sim + '_t-1.dat' - data = load_aiolos_diag(filename) - + data = load_aiolos_diag(filename, species=1) + Tirr40 = 0.25*data['S0'][-1]/sigma_rad gamma0 = (data['kappa_PTsun0_0'] / data['kappa_P0_0'])[-1] - gamma1 = (data['kappa_PTsun1_0'] / data['kappa_P1_0'])[-1] + if 'S1' in data.dtype.names: + gamma1 = (data['kappa_PTsun1_0'] / data['kappa_P0_0'])[-1] + Tirr41 = 0.25*data['S1'][-1]/sigma_rad + else: + gamma1 = 1 + Tirr41 = 0 - Tirr40 = 0.25*data['S0'][-1]/sigma_rad - Tirr41 = 0.25*data['S1'][-1]/sigma_rad + tau = data['tau_radial0']/gamma0 - tau = data['tau_radial0'] - T_guillot = Guillot_2band(tau, Tint**4, Tirr40, gamma0, Tirr41, gamma1) + # Get the planet's internal temperature + pars = load_aiolos_params(f"{sim}.par") + Tint = 0 + if 'USE_PLANET_TEMPERATURE' in pars: + if 'PARI_TINT' in pars: + Tint = float(pars['PARI_TINT']) + else: + Tint = float(pars.get('PARI_TPLANET',0)) + + # Get the flux-limiter factor + fH = 1/float(pars.get('XI_RAD',1)) + + T_guillot = Guillot_2band(tau, Tint**4, Tirr40, gamma0, Tirr41, gamma1, + fH=fH) idx = (tau < 100) & (tau > 1e-10) L1 = np.abs(1 - data['T_gas0'][idx]/T_guillot[idx]).mean() if L1 <= L1_target: - print('Irradiation test L1 check passed') + print(f'Irradiation test {sim[-1]} L1 check passed') else: - print("Irradiation test L1 check failed. " + + print(f"Irradiation test {sim[-1]} L1 check failed. " + "L1={}, expected={}".format(L1,L1_target)) @@ -50,7 +69,7 @@ def check_structure(sim, Tint, L1_target, make_plots): if make_plots: import matplotlib.pyplot as plt - plt.semilogy(data['T_gas0'], tau, label='sim') + plt.semilogy(data['T_gas0'], tau, label=f'{sim}') plt.loglog(T_guillot, tau, 'k--', label='Guillot (2010), $f_H=1$') ylim = plt.ylim() @@ -63,13 +82,21 @@ def check_structure(sim, Tint, L1_target, make_plots): plt.legend(loc='best') os.makedirs('plots', exist_ok=True) - plt.savefig('plots/Irradiation.png') - + plt.savefig(f'plots/{sim}.png') + plt.clf() + np.testing.assert_array_less(L1, L1_target) def test_structure(make_plots=False): - check_structure('irradiation', 175.0, 1.5e-3, make_plots) + par_files = [ + ("irradiation1", 0.02), + ("irradiation2", 0.06) + ] + + for f, tol in par_files: + check_structure(f, tol, make_plots) + #check_structure('irradiation', 175.0, 1.5e-3, make_plots) if __name__ == "__main__": From 94a3f3f5edaccaa2448953ec966556bdebdebcee Mon Sep 17 00:00:00 2001 From: Richard Booth Date: Wed, 5 Jul 2023 16:47:58 +0100 Subject: [PATCH 10/10] removed old test file --- test_files/irradiation.par | 78 -------------------------------------- 1 file changed, 78 deletions(-) delete mode 100644 test_files/irradiation.par diff --git a/test_files/irradiation.par b/test_files/irradiation.par deleted file mode 100644 index 4c9c71e..0000000 --- a/test_files/irradiation.par +++ /dev/null @@ -1,78 +0,0 @@ -# -# Parameterfile for atmosphere simulations -# Do not use tabs, only white spaces as value separators -# -DENSITY_FLOOR 1e-20 -PARI_TIME_TMAX 1.01e8 -CFL_BREAK_TIME 1e100 -PARI_TIME_OUTPUT 2.e7 -PARI_CFLFACTOR 1e5 -PARI_TIME_DT 1e10 -PARI_DPHI_FACTOR 1.0 -ENERGY_EPSILON 0.01 -# Control parameters -# -# Problem numbers explained: shock_tube, no grav == 1, planet with grav == 2# -# Friction solver: analytic == 0, numerical == 1 -# Collision model: constant == C, physical == P -PARI_PROBLEM_NUMBER 2 -PARI_NUM_SPECIES 1 -PARI_NUM_BANDS 2 -DO_HYDRO 0 -#CONST_TEMP 1430. -PARI_CONST_TEMP 730. -INIT_TEMPERATURE_MODEL C -INIT_J_FACTOR 1e0 -INIT_T_TEMP 730 - -# Solar radiation temps at 0.1 AU: 645K for bond albedo=0.75, 913K for bond albedo = 0 -# print 5777*(1./20./2.)**0.5*(1-0.75)**0.25 - -FRICTION_SOLVER 1 -PARI_COLL_MODEL P -PARI_ALPHA_COLL 1. -PARI_USE_COLL_HEAT 0 -PARI_ORDER 2 -PARI_USE_RADIATION 1 -PARI_OPACITY_MODEL C -PRESSURE_BROADENING 0. -CONSTOPA_SOLAR_FACTOR 1e0 -PARI_CONST_OPAC 1.e0 - -PARI_SELF_GRAV_SWITCH 0 -PARI_LINEAR_GRAV 0 -PARI_DEBUGLEVEL 1 -PARI_INIT_STATIC 1 -USE_TIDES 1 -#Physical parameters - -PARI_PLANET_MASS 224. -PARI_INIT_DATA_SWITCH 1 -PARI_INIT_DATA_U1 1e-68 -PARI_INIT_DATA_U2 1.0 -PARI_INIT_DATA_U3 1e-68 - -#Radiation parameters -PARI_LAM_MIN 9.0e-2 -PARI_LAM_MAX 9.0e-2 -PARI_LAM_PER_DECADE 0.1 -PARI_TSTAR 6070. -PARI_RSTAR 0.1 -PARI_UVSTAR 1.e31 - -PARI_PLANET_DIST 0.05 - -#Numerical parameters -PARI_BOUND_TYPE 2 -PARI_BOUND_TYPE_LEFT 2 -PARI_BOUND_TYPE_RIGHT 1 -PARI_DOMAIN_DX 1e-3 -PARI_GRID_TYPE 1 -PARI_GEOMETRY 2 -PARI_CELLS_PER_DECADE 2000 -#Domain min/max are 8 me core radius for earth density, and 100*hill radius -PARI_DOMAIN_MIN 9.200e9 -PARI_DOMAIN_MAX 1.0e11 - -#Scenarios -PARI_INIT_WIND 1