[aspect-devel] Questions about Integrating the GyPSuM Tomography Model into ASPECT

[Wolfgang Bangerth]

David,

>> When testing these things, it's ok to hard-code the inner and outer radius for 
>> a moment. You can later generalize this by reading parameters from some model 
>> at run-time, but first make sure that the results you get are actually correct 
>> before you start making the code more complicated.
> 
> We did hardcode the inner and outer radii for testing purposes, but I 
> figured I'd ask about that functionality as part of this thread since it 
> was something that was going to have to be resolved at some point. I'm 
> also generally curious about the proper way to call parameters that are 
> being read by other modules.

The way we typically do this is to query the geometry model. As a simple 
example, take a look at the 'box' initial temperature model in 
source/initial_temperature/box.cc -- it has this kind of code that 
specifically looks at properties of the box (namely, the get_origin() 
and get_extents() methods in the last line):

     initial_temperature (const Point<dim> &position) const
     {
       // this initial condition only makes sense if the geometry is a
       // Box. verify that it is indeed
       const GeometryModel::Box<dim> *geometry
         = dynamic_cast<const GeometryModel::Box<dim>*> 
(&this->get_geometry_model());
       AssertThrow (geometry != 0,
                    ExcMessage ("This initial condition can only be used 
if the geometry "
                                "is a box."));

       double perturbation = 1;
       for (unsigned int d=0; d<dim; ++d)
         perturbation *= 
std::sin(numbers::PI*(position[d]-geometry->get_origin()[d])/geometry->get_extents()[d]);


There are other models that are specific to the spherical shell geometry 
and query the inner and outer radii.


>> * When you run the code on one processor, does it run through in debug mode
>>  without errors?
> 
> When all of the calculations were in the template, the code would 
> compile in debug mode, but the coefficients were inaccurate and ASPECT 
> would crash shortly after the runs were started (the temperature 
> perturbations were being reported as inf).

Can you explain what the 'template' is? And did you figure out why the 
perturbations were reported as 'inf'?


> After moving the calculations 
> to the SphericalHarmonicsLookup class, the code compiles and the 
> coefficients are correct (any depth-dependent calculations were done at 
> a hardcoded depth rather than at the position in the mesh that ASPECT is 
> populating). However, I haven't been able to test whether ASPECT will 
> run with these calculations because I don't know how to pass the mesh 
> position to the class (so these coefficients are only being tested at 
> one point).

I have to admit that I don't know how the SphericalHarmonicsLookup class 
works. But how does the S40RTS or SAVANI models do this?


> I make it about halfway down the list. If I leave the code in the 
> template, I gain access to the mesh position, but the coefficients are 
> inaccurate.

That would be the place to leave it then and to debug it. If you don't 
understand why the results are wrong there, you're likely only going to 
make things worse by moving the code to another place that may be even 
more complicated to understand.

You may not need this advice, but since we're on a public mailing list 
where others can learn from these sorts of exchanges, I'll offer it 
anyway: When I identify something that looks like a bug, I typically try 
to debug it as soon as I can find it, rather than work around it or move 
the code to some other place. That's because a bug typically represents 
a situation that I thought I understand, but *apparently do not*. Unless 
I intend to completely rewrite the code, my lack of understanding of 
what is going on is not going to disappear and experience shows that 
sooner or later, the bug will reappear and I will have to sit down and 
deal with it anyway. So my strategy is to do it right away -- if I find 
that I don't understand why something doesn't work, I try to figure it 
out right there and then. When I've made it work, at least I have a 
baseline to compare against if I want to move the code somewhere else or 
rewrite it. (This baseline typically comes in the form of a testcase 
that goes into the testsuite.)

Best
  W.

-- 
------------------------------------------------------------------------
Wolfgang Bangerth          email:                 bangerth at colostate.edu
                            www: http://www.math.colostate.edu/~bangerth/