Wind test including radiation force + ieos = 17 + bug fixes #542
Conversation
…a dusty (bowen) wind with radiative force included
|
this all looks good, remaining failure is just that the tolerance needs to be a bit higher for the internal energy: |
src/tests/test_wind.f90
Outdated
| call checkval(npart,12180,0,nfailed(3),'number of ejected particles') | ||
| call checkval(xyzmh_ptmass(15,1),1.591640703559762E-06,epsilon(0.),nfailed(4),'wind mass loss rate') | ||
| if (testcyl) then | ||
| call checkval(eint,3.360686893182378E+03,eps_sum,nfailed(5),'total internal energy') |
There was a problem hiding this comment.
It would really help here to know where these numbers come from, e.g. if we improve things how do we know if it is better? Are these just initial values, or are they from some semi-analytic solution?
There was a problem hiding this comment.
There is no analytical solution and these numbers are associated to a non converged simulation. To reach convergence on the analytical solution would require very high resolution and very long integration time.
This is just a check that the code consistently gives the same numbers.
If the timestepping or the discretisation of the equations (e..g 2nd or fourth order) change the numbers will be completely different (relative error can exceed several percents).
The tests on Eint and Ekin were not included in the previous version of test_wind.
src/tests/test_wind.f90
Outdated
| call checkval(eint,3.367417540822784E+03,eps_sum,nfailed(5),'total internal energy') | ||
| call checkval(ekin,5.524867074648306E+01,eps_sum,nfailed(6),'total kinetic energy') | ||
| else!if (test) | ||
| call checkval(eint,3.179016341424608E+03,eps_sum,nfailed(5),'total internal energy') |
There was a problem hiding this comment.
same comment, can these be calculated from first principles or are they just regression tests?
There was a problem hiding this comment.
they are just regression tests
| eint = sum(vxyzu(4,1:npart)) | ||
| ekin = sqrt(sum(vxyzu(1,1:npart)**2+vxyzu(2,1:npart)**2+vxyzu(3,1:npart)**2)) | ||
| if (vb) print '(5(1x,es22.15),i8)',xyzmh_ptmass(4,1),xyzmh_ptmass(7,1),xyzmh_ptmass(15,1),eint,ekin,npart | ||
| call checkval(xyzmh_ptmass(4,1),1.199987815414834E+00,epsilon(0.),nfailed(1),'sink particle mass') |
There was a problem hiding this comment.
should be m_ini-m_inj
| ! check 1D wind profile | ||
| i = size(trvurho_1D(1,:)) | ||
| !print '((5(1x,es22.15)))',trvurho_1D(:,i),massoftype(igas) | ||
| call checkval(trvurho_1D(2,i),7.058624412798283E+13,epsilon(0.),nfailed(2),'outer wind radius') |
There was a problem hiding this comment.
same comment with values here, can they be computed or are these just regression tests also?
| call checkval(xyzmh_ptmass(4,1),1.199987894518367E+00,epsilon(0.),nfailed(6),'sink particle mass') | ||
| call checkval(xyzmh_ptmass(7,1),0.,epsilon(0.),nfailed(7),'mass accreted') | ||
| call checkval(npart,12180,0,nfailed(8),'number of ejected particles') | ||
| call checkval(xyzmh_ptmass(15,1),1.591640703559762E-06,epsilon(0.),nfailed(9),'wind mass loss rate') |
There was a problem hiding this comment.
is this prescribed somewhere?
There was a problem hiding this comment.
xyzmh_ptmass(15) is the mass loss rate. It is defined in the .in file. Don't know why I checked this variable. It is useless.
Type of PR:
Bug fix / modification to existing code (fixes #540)
Description:
Include a new series of tests in the test_wind module :
Implement ieos = 17 to keep gamma constant but allow the mean molecular weight to change according to H2 formation
Bug fixes
Testing:
run simple wind models with ieos=5,17 to check that mu and gamma where not held constant
Did you run the bots? no
Did you update relevant documentation in the docs directory? no