Draft: Add magma2 SPH scheme

I just added the MAGMA2 scheme in the branch MAGMA2 and are there some advices for the future implementation?

I added the two extra profiles:

• hydro_iact_cosmo.h and hydro_iact_sphenix.h. One is added the cosmological term and another is using the sphenix's diffusion parameter.
• There is a parameter called h_crit in the runner_iact_force and runner_iact_nonsym_force , which is the inverse of the resolution eta. And everytime I just change the parameter manually
• Also for the matrix inversion, I just paste the function invert_dimension_by_dimension_matrixin swift, and it might need clear a bit and call the function instead
• And I just directly use constant viscosity and diffusion parameters into the function without define it in the hydro_parameter.h. It also might need to change a bit
• Last thing is MAGMA2 use the quadratically mid point reconstruction and it is also worth to test the linearly reconstruction. Since it may increases the computing effciency.

requested review from @bvandenbroucke

changed the description

added SPH enhancement labels

@jborrow will also be interested.

marked this merge request as draft

changed title from Draft: add magma2 scheme to Draft: Add magma2 SPH scheme

added 1 commit

• 5279c958 - Also add the new files to the outside build

Compare with previous version

Also I added the MAGMA2 scheme based on the SPHENIX scheme, and I didnt delete the SPHENIX part in all the profiles

added 1 commit

• 189be751 - Applied formatting script.

Compare with previous version

Please use 1.f instead of 1 in these expressions.

Could any of these denominators be 0?

changed this line in version 6 of the diff

Change the value into float.

All these denominators can be zero, so that I write a few lines to prevent the issue.

if (A_i == 0.f && A_j == 0.f){
	  A_ij = 1.f;
	  A_ji = 1.f;
  } else if ((A_i == 0.f && A_j != 0.f) || (A_j == 0.f && A_i != 0.f)){
	  A_ij = 0.f;
	  A_ji = 0.f;
  } else{
	  A_ij = A_i / A_j;
	  A_ji = A_j / A_i;
  }

  if (Av_i == 0.f && Av_j == 0.f){
          Av_ij = 1.f;
          Av_ji = 1.f;
  } else if ((Av_i == 0.f && Av_j != 0.f) || (Av_j == 0.f && Av_i != 0.f)){
          Av_ij = 0.f;
          Av_ji = 0.f;
  } else{
          Av_ij = Av_i / Av_j;
          Av_ji = Av_j / Av_i;
  }

In what scenario can these actually be zero? Is that preventable in a different way?

if (A_i == 0.f && A_j == 0.f) { /* For smooth internal energy field, we turn off diffusion term*/
    A_ij = 1.f;
    A_ji = 1.f;
  } else if ((A_i == 0.f && A_j != 0.f) || (A_j == 0.f && A_i != 0.f)) { /* For extreme values, we add diffusion term*/
    A_ij = 0.f;
    A_ji = 0.f;
  } else {
    A_ij = A_i / A_j;
    A_ji = A_j / A_i;
  }

When A_i and A_j are both zero, it means two particles in the smoothed field. We will set A_ij and A_ji to 1, and it will turn off the diffusion term, same as the viscosity term. When one of A_i/A_j is zero, it encounters discontinuity or some extreme values. We set A_ij/A_ji to zero, so that it will add diffusion term to fix the discontinuity.

This looks very good! I have left some comments for you. Since there is a lot of duplicate code in hydro_iact, some comments apply to multiple lines and I have not always put them everywhere.

Can you explain to me again what the difference is between hydro_iact, hydro_iact_sphenix and hydro_iact_cosmo and why we have all 3 of them? I think only hydro_iact is actually used? A quick diff shows that these files are almost the same. There are probably better ways to have different versions without having to duplicate all the code.