Adiabatic density profile #1274
Conversation
gassmoeller
changed the title
gassmoeller
force-pushed
the
adiadatic_density_profile
branch
from
443b618
to
9b239e2
Compare
| template <int dim> | ||
| double Function<dim>::density_derivative (const Point<dim> &p) const | ||
| { | ||
| // TODO: better eps or make it a user input |
There was a problem hiding this comment.
I am still not sure if we should require the user to specify this derivative. Thoughts?
| if (z < delta_z * std::numeric_limits<double>::epsilon()) | ||
| return temperatures.front(); | ||
|
|
||
| const unsigned int i = static_cast<unsigned int>((z/delta_z) * (1. - 2. * std::numeric_limits<double>::epsilon())); |
There was a problem hiding this comment.
We are now repeating this logic many times in this file. How about moving this into a helper function?
| template <int dim> | ||
| double Function<dim>::density_derivative (const Point<dim> &p) const | ||
| { | ||
| // TODO: better eps or make it a user input |
There was a problem hiding this comment.
I am still not sure if we should require the user to specify this derivative. Thoughts?
| @@ -187,23 +201,80 @@ namespace aspect | |||
| { | |||
| const double z = this->get_geometry_model().depth(p); | |||
|
|
|||
| Assert (z >= 0, ExcInternalError()); | |||
There was a problem hiding this comment.
Isn't this too strict? Maybe we should allow being eps below 0. You catch this case below anyways (by returning temperatures.front).
| @@ -300,14 +321,19 @@ namespace aspect | |||
| "The format in which the output shall be produced. The " | |||
| "format in which the output is generated also determines " | |||
| "the extension of the file into which data is written."); | |||
| const std::string variables = | |||
There was a problem hiding this comment.
this is a great idea, I would have done the same. ;-)
gassmoeller
force-pushed
the
adiadatic_density_profile
branch
2 times, most recently
from
e902f6c
to
5f62db4
Compare
| in.temperature[0] = temperatures[i-1]; | ||
| in.pressure[0] = pressures[i-1]; | ||
| in.temperature[0] = temperatures[i]; | ||
| in.pressure[0] = pressures[i]; | ||
| in.velocity[0] = Tensor <1,dim> (); |
There was a problem hiding this comment.
can we move this and the strainrate, etc. outside the loop?
|
|
||
| const double floating_index = z/delta_z; | ||
| const unsigned int i = static_cast<unsigned int>(floating_index); | ||
| if (std::abs(floating_index-std::round(floating_index)) < 1e-6) |
There was a problem hiding this comment.
can you add a comment here please?
| if (std::abs(floating_index-std::round(floating_index)) < 1e-6) | ||
| return property[i]; | ||
|
|
||
| Assert (i+1 < temperatures.size(), ExcInternalError()); |
| // now do the linear interpolation | ||
| const double d=1.0+i-z/delta_z; | ||
| Assert ((d>=0) && (d<=1), ExcInternalError()); | ||
| const unsigned int i = static_cast<unsigned int>((z/delta_z) * (1. - 2. * std::numeric_limits<double>::epsilon())); |
There was a problem hiding this comment.
can you add a comment that says something like "if z/delta_z is within [k-eps, k+eps] of a whole number k, round it down to k-1"
| @@ -42,7 +42,7 @@ namespace aspect | |||
| // note: the pressure case for adiabatic conditions initialized and | |||
| // the pressure case for non-initialized adiabatic conditions | |||
| // are the same, so the pressure case is included in the if statement for adiabatic conditions. | |||
| if (this->get_adiabatic_conditions().is_initialized() || use_depth ==false) | |||
| if (use_depth ==false) | |||
There was a problem hiding this comment.
Can you fix the spacing around == or use !use_depth?
gassmoeller
changed the title
|
I addressed the comments. For now I would leave the TODO about the better epsilon for the function derivative to make progress with the other parts. |
gassmoeller
force-pushed
the
adiadatic_density_profile
branch
2 times, most recently
from
f101078
to
a9820f2
Compare
In order to support simplified density approximations like BA, TALA, ALA we need a reference density profile. This is a first step introducing the profile, which will be used by another PR by @tjhei later on. I expect many test failures, because this PR also solves a small bug in
AdiabaticConditions::InitialProfile, in which we used a shifted index for the density to compute the next step. The effect is small (10^-4 on temperature and pressure). I will update the PR, when this is ready for review.