mass loss rate in wind is not equal to the value specified in the input file

[danieljprice]

Using SETUP=wind, there appears to be a factor of two offset between the actual mass loss rate (mass injected in a simulation over a period of time; e.g. by measuring decrease in sink particle mass) and the desired wind_mass_rate in the input file.

From experiments with Taissa Danilovich the two numbers appear to be discrepant by exactly a factor of 2, as in actual Mdot is a factor of two higher than the desired Mdot.

Would be good to get to the bottom of this.

[lsiess]

I fixed it. I will soon make a pull request

[danieljprice]

What was the cause?

[lsiess]

On Thu, 06 Apr 2023 03:52:07 -0700 Daniel Price @.***> wrote:

What was the cause?

mass was removed twice. In inject_particles we had the following sequence:

do i=inner_boundary_sphere,outer_sphere,-1

.....

   ! update the sink particle mass
   if (nptmass > 0 .and. wind_emitting_sink <= nptmass) then
      xyzmh_ptmass(4,wind_emitting_sink) = xyzmh_ptmass(4,wind_emitting_sink) - mass_of_spheres
   endif
   print '(" ##### eject sphere ",4(i4),i7,9(1x,es12.5))',i,inner_sphere,nboundaries,&
        outer_sphere,npart,time,local_time,r/xyzmh_ptmass(iReff,1),v,u,rho
endif
!cs2max = max(cs2max,gamma*(gamma-1)*u)

enddo ! update sink particle properties if (nptmass > 0 .and. wind_emitting_sink <= nptmass) then mass_lost = mass_of_spheres (inner_sphere-outer_sphere+1) xyzmh_ptmass(4,wind_emitting_sink) = xyzmh_ptmass(4,wind_emitting_sink) - mass_lost if (pulsating_wind) then inner_radius = wind_injection_radius + deltaR_oscsin(omega_osc*time) !v2 xyzmh_ptmass(5,wind_emitting_sink) = (inner_radius3-dr3)(1./3.) !accretion radius xyzmh_ptmass(5,wind_emitting_sink) = inner_radius endif endif

As you see, mass was removed twice from the wind emitting sink particle. I suspect this error was introduced in a merge. It is not the first time I am discovering unpleasant mistakes like this after a merge...

Cheers Lionel

[danieljprice]

Am closing the issue but what this says to me is we should have a test for this. Ideally we should add the tests from Siess+ to the testsuite and run them at very low resolution on each pull request, checking exactly this kind of thing cannot be broken