Bugfixes for init conds

Author: FredJenningsAstro

src/dust/dust.hpp

@@ -1328,11 +1328,7 @@ if(!(gamma_star > -3.001 && gamma_star < -2.999)){ // Guard against gamma_star =
 Real norm_n = stellar_mass_cent * gthr;
 Real norm_d = 4. * Kokkos::numbers::pi * (Kokkos::pow(stellar_density_profile_r_up, gthr) - Kokkos::pow(stellar_density_profile_r_low, gthr));
 norm = norm_n/norm_d;
-
-
-
 } else {
-printf("Doing log routine! \n");
  norm = stellar_mass_cent * 1./ (4. * Kokkos::numbers::pi * Kokkos::log(stellar_density_profile_r_up/stellar_density_profile_r_low));
 }
 // printf("r = %e norm = %e gamma_star = %e stellar_density_profile_r_low =%e stellar_density_profile_r_up=%e stellar_mass_cent=%e\n", r, norm, gamma_star, stellar_density_profile_r_low, stellar_density_profile_r_up, stellar_mass_cent);
@@ -1866,10 +1862,6 @@ for(int gc_i = 0; gc_i < num_grain_compositions; gc_i ++ ){
         total_mass_over_whole_dist_and_comps +=  (added_silicate_mass *  agb_normalised_silicate_mass_distribution_array[gs_i]);
       }
 
-
-
-
-
   }
 }
 

src/pgen/cluster.cpp

@@ -1065,17 +1065,15 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData<Real> *md) {
                 Real total_mass_C = 0; // Total Mass, not a density
                 Real total_mass_S = 0; // Total Mass, not a density
                 Real stellar_mass_this_cell = 0;
-                Real total_mass_over_whole_dist_and_comps = 0;
-                DustCalculateAGBWindContribution(total_mass_C,total_mass_S, total_mass_over_whole_dist_and_comps, stellar_mass_this_cell, b, k, j, i, 
+                DustCalculateAGBWindContribution(total_mass_C,total_mass_S, total_dust_mass, stellar_mass_this_cell, b, k, j, i, 
                   cons_pack, DustDevObj, init_run_stellar_injection_time);
                   // printf("A) total_mass_C=%e total_mass_S=%e \n",total_mass_C,total_mass_S);
                   // insert the amount of dust added by stellar injection, or the minimum amount allowed by the floor:
 
 
-                total_dust_mass = total_mass_C + total_mass_S; 
 
-                if(total_mass_over_whole_dist_and_comps > 1e-50){
-                Real renorm_fractor_for_dtgfloor = std::max(total_mass_over_whole_dist_and_comps, dtgfloor*u(IDN, k, j, i)*volume ) / total_mass_over_whole_dist_and_comps;
+                if(total_dust_mass > 1e-50){
+                Real renorm_fractor_for_dtgfloor = std::max(total_dust_mass, dtgfloor*u(IDN, k, j, i)*volume ) / total_dust_mass;
                 total_mass_C *= renorm_fractor_for_dtgfloor;
                 total_mass_S *= renorm_fractor_for_dtgfloor;
                 //if total_mass_C == 0 or total_mass_S == 0 we are outside range of AGB winds, so need to be clever how we get the S/C ratio
@@ -1087,7 +1085,7 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData<Real> *md) {
                     total_dust_mass_this_gc_i = total_mass_S;
                         }
                   } else {
-                    total_dust_mass_this_gc_i = total_dust_mass; // Only one of these will be non-zero in this case if we only have 1 composition
+                    total_dust_mass_this_gc_i = total_dust_mass *renorm_fractor_for_dtgfloor; // Only one of these will be non-zero in this case if we only have 1 composition
                   } 
                 } else {