[aspect-devel] How to grab the radius of every interior sphere made up from cell faces in the domain and some MPI questions

[Shangxin Liu]

Hi Rene,

I tested my code in a simple degree 3 order 2 harmonic perturbation case. I
think you also do that in your code. In my code there is a main
execute template performing the computation of geoid from spherical
harmonic expansions. There are other three templates each as a function to
return the spherical harmonic coefficients of surface topo, CMB topo, and
density anomaly. There is also another template as a function to perform
the data transformation to spherical harmonic coefficients and return the
values. Maybe this is not the best way, but it now works in the single
processor for the test.

As for grabbing different layers' radii, you're right, my approach for now
is exactly

vector R_sphere initialized with inner radius and outer radius, (r_CMB,
r_surface);
Then loop through all cells
  if (the radius of the midpoint of the upper face is already in the vector)
    skip cell;
  else
    push_back radius into R_sphere;
reorder the elements in R_sphere from the smallest r_CMB to largest
r_surface.
Then perform the numerical integral of rho_l_m(r) through each layer.

This approach now works fine in a single processor (I've took care of the
floating point math in my code), but doesn't work in parallel since each
processor will only cover a sub-domain and in this case each processor will
hold a vector R_sphere and each of them may not include all the radius
values throughout the mantle. At the beginning, I once considered directly
searching for the radii with equal division, 'minimal_radius + (n /
number_of_layers) * model_depth' with
n being the number of the layer I'm in. However, I found that in 0 initial
global refinement degree, the output three radii (3480,4501,6371) are not
in equal interval. This puzzled me a bit and made me resort to the above
"awkward" way. So how can I solve this problem in the parallel code?

As for the adaptive refinement problem, yes, from the beginning I've
realized this "aligning layers with the cell's boundary" approach will
break down in AMR. But for now I still don't know the detail of the
adaptive mesh refinement in ASPECT so I've never used this parameter in my
global model. Maybe We should use a better approach to also accommodate AMR.

About the second problem, let me describe in detail what I'm doing:

In my previous template(function) to compute spherical harmonic
coefficients, to deal with the term sin(theta)*delta_theta*delta_phi (theta
is polar and phi is the azimuth), I directly regard each cell's upper face
on the sphere as a small region for the sphere integral, and use the area
of each cell's upper face over the square of the current sphere to compute,
i.e., sin(theta)*delta_theta*delta_phi = delta_area/(r_sphere^2). However,
this approach works fine in higher refinement degree but generates a
considerable error in the coarse refinement. So now I make a new function
to generate a new 2D sphere grid for theta and phi to compute sphere
integral (theta = 0:size:180; phi = 0:size:360). The size of the grid is
controlled by an input parameter. In this template, I need to make some
kind of 2D sphere interpolation to get the value at the new grid point from
the output value at each cell upper face's midpoint. For now I just use the
simplest "nearest interpolation" to find the closest cell upper face
midpoint to the new grid point and then set the value at the new grid point
to that of the closest cell upper face midpoint. In this way, the result in
the lower refinement degree indeed becomes more accurate. However, this
induces another similar MPI problem as above. Each processor will only
cover part of the total cells so I need to further compare between
different processor to find the closest one among all the local_closest
results and also find the cell upper face's midpoint of that "closest"
processor.

Hope this further clarify the problems and my workflow,

Best,
Shangxin

>
>
> Message: 1
> Date: Thu, 22 Oct 2015 11:31:44 -0500
> From: Rene Gassmoeller <rengas at gfz-potsdam.de>
> To: aspect-devel at geodynamics.org
> Subject: Re: [aspect-devel] How to grab the radius of every interior
>         sphere made up from cell faces in the domain and some MPI questions
> Message-ID: <56290F70.80402 at gfz-potsdam.de>
> Content-Type: text/plain; charset=utf-8
>
> Hi Shangxin,
> I think I can help you with some of your questions, but first:
> I think I had mentioned at the Hackathon and also at one of the video
> user meetings earlier this year that I was working together with some
> others on a geoid postprocessor already. At the hackathon and afterwards
> Jacky, Ian, Anthony Osei-Tutu (a PhD student with us in Potsdam) and
> myself put something together that is already close to working. Since
> the project that gave us the incentive to write the postprocessor
> developed in a slightly different way it was never brought into a shape
> where we opened a pull request for the main version (it was a branch on
> my github repository). I did not know you where working on something
> similar, otherwise I would have opened one. But I did that now, so you
> can have a look at:
>
> https://github.com/geodynamics/aspect/pull/648
>
> The code is not exactly finished, we never benchmarked it properly and
> it gives some curious patterns for low resolutions, but it seems to work
> ok for higher resolutions. Maybe you have more experience with testing
> geoid output or can ask Scott for some testcases?
>
> >From what I read in your mail, we might be able to merge the best of the
> two versions into something that actually makes it into the main code.
> Our approach was mostly to take the algorithm used in CitcomS and
> transfer it to ASPECT with as little content changes as possible, while
> still putting it into our usual code style. Instead of a lengthy
> description of our approach, let me just try to answer your questions to
> point out what we do to circumvent your problems.
>
> > After getting the
> > topography value of each cell along the boundary, there is a another
> > template in my code to transfer them to the form of spherical harmonic
> > coefficients.
>
> I am not sure what you mean by template here. In most places we use only
> one template parameter for most of our classes, which is the dimension.
> Do you mean function?
>
> >
> > As for the density anomaly, The formula needs to first get the density
> > anomaly value per cell of each depth layer (in my code, the number of
> total
> > sphere layers is determined by initial refinement degree, i.e, there are
> 2
> > layers in 0 refinement degree, each sphere layer has 48 density values,
> one
> > value per cell).
>
> I think this approach is nice for a uniform mesh, but it might break
> down for AMR. Even if you start with only global refinement, your mesh
> might get coarsened over time. In our plugin we therefore follow the
> approach of the depth average postprocessor, in which we make the number
> of layers an input parameter for the user. Most likely the user will
> know best, what the 'average' resolution of the model will be, and what
> resolution might be reasonable for the number of layers. This of course
> poses the problem that some cells are partly in one layer and partly in
> another. For now, we decided to simply attribute the cell to the layer,
> which contains its midpoint. The same problem will occur in your
> approach when you switch on AMR.
>
> Then I transfer the density anomaly values for each layers
> > into spherical harmonic coefficients, rho_l_m(r), and integral the
> > coefficients from r_CMB to r_surface (in the 0 refinement degree, the
> > numerical integral contains two layers, 3480 to 4501km and 4501 to 6371
> > km). So in this integral, I need to grab the radius value of each
> internal
> > sphere, i.e, the radius of the cell's upper face globally (in simplest
> > refinement degree 0 case there are three radius 3480,4501,and 6371). My
> > current approach is to loop over all the cells and get the location of
> each
> > cell's upper face midpoint, then compute the distance between earth
> center
> > and these midpoint, remove the same distance value of the cell at the
> same
> > depth, finally push_back them to a "R_sphere vector" ((3480,4501,6371)
> in 0
> > refinement degree). But I think this approach is not an elegant way and
> > also will not work in parallel since each processor only holds a
> > sub-domain. So is there a fast way to directly grab the radius of each
> > internal sphere determined by the cell face in deal.ii/ASPECT?
>
> This is indeed more complicated in parallel and also you need to be
> careful about floating-point math here. As far as I understand you ask:
> for all cells
>   if (the radius of the midpoint of the upper face is already in the
> vector)
>     skip cell
>   else
>     add radius
>
> However, the above mentioned approach will break down with an adaptive
> mesh, because there might be more refined cells only in a part of the
> domain. I guess the problem you are encountering mostly stems from the
> fact that you want to align your layers with the cell boundaries
> (something we did not enforce, although it might be more accurate).
> But another approach could be the following: The cells divide the total
> model depth into equal parts. Therefore, the vector you are searching is
> essentially 'minimal_radius + (n / number_of_layers) * model_depth' with
> n being the number of the layer you are in. No need for looping over
> cells, or comparing radii of the face midpoints. This is the approach we
> used.
>
> >
> > My second problem is the comparison between the results of different
> > processor. To be concise, there is a loop through all the cells to find
> > which cell's upper face midpoint is the closest to a fixed point.
>
> What is this fixed point? And what do you want to do there? If you just
> want to evaluate the solution at that point, consider using the
> FEFieldFunction class of deal.II to get the solution at an arbitrary
> point in the domain. If you indeed need to find the cell that contains
> this point (and the closest face midpoint will always be the one of this
> cell -- or one of the neighbouring cells, if the neighbours are more
> refined than the containing cell) then consider using the
> GridTools::find_active_cell_around_point function.
>
> Hope that helps. Best regards,
> Rene
>
> >
> > On Wed, Oct 21, 2015 at 5:34 PM, Wolfgang Bangerth <bangerth at tamu.edu>
> > wrote:
> >
> >>
> >> Shangxin,
> >> I have difficulty following your questions, and I believe it is because
> >> you have already picked a specific algorithm and you are asking how to
> >> implement it. But we know neither the algorithm, nor what the problem is
> >> you are trying to solve.
> >>
> >> Can you give us an overview *in formulas* of what you are actually
> trying
> >> to compute, along with an overview of *how* you try to compute it? Maybe
> >> that already helps us give you suggestions how to avoid some of the
> >> problems you see.
> >>
> >> Best
> >>  W.
> >>
> >>
> >> On 10/21/2015 03:23 PM, Shangxin Liu wrote:
> >>
> >>> Hi;
> >>>
> >>> I'm now making a new plugin into the postprocess directory to compute
> >>> geoid anomaly. The code has been finished but it needs some further
> >>> modification to run in parallel. To realize this, I run into these two
> >>> problems:
> >>>
> >>> 1. The geoid computation needs the integral of the density anomaly from
> >>> CMB to surface. For now, I do the integral based on cells. So to get
> all
> >>> the radii of the internal spheres, the code loops over all of the cells
> >>> and get the radius of each cell's upper face and push_back the radius
> >>> into a vector after removing the same values from the cell faces at the
> >>> same depth. But I suppose this solution is not so elegant and also will
> >>> not work in the multiple MPI processors since each processor will only
> >>> cover part of the domain. So is there a fast way in deal.ii/ASPECT to
> >>> directly grab all the radii of the spheres determined by the cell faces
> >>> within the domain?
> >>>
> >>> 2. In the spherical harmonic coefficients computation, I need to
> >>> evaluate the values at each grid points I created prepared for sphere
> >>> integral from the output values at the center of each cell face. My
> >>> current solution is to loop over all the sphere cells and find the
> >>> closest cell face midpoint to the grid point (nearest point
> >>> interpolation). But since each processor will only cover part of the
> >>> domain so I also need to compare among all the processors to find the
> >>> closest point and the corresponding density anomaly/topo values at that
> >>> point. I notice that there is a Utilities::MPI::max/min command can
> >>> compare among different processors. But my purpose is not only to get
> >>> the closet distance, but also need to know which processor provide this
> >>> closest cell face midpoint and the corresponding location of this
> closet
> >>> cell midpoint hence to get the density anomaly/topo there. Are there
> any
> >>> MPI commands in deal.ii/ASPECT to realize this?
> >>>
> >>> Best,
> >>> Shangxin
> >>>
> >>
> >>
> >> --
> >> ------------------------------------------------------------------------
> >> Wolfgang Bangerth               email:
> bangerth at math.tamu.edu
> >>                                 www:
> http://www.math.tamu.edu/~bangerth/
> >>
> >>
> >
> >
> >
> > _______________________________________________
> > Aspect-devel mailing list
> > Aspect-devel at geodynamics.org
> > http://lists.geodynamics.org/cgi-bin/mailman/listinfo/aspect-devel
> >
>
>
> ------------------------------
>
> Message: 2
> Date: Thu, 22 Oct 2015 12:32:29 -0400
> From: Shangxin Liu <sxliu at vt.edu>
> To: Wolfgang Bangerth <bangerth at tamu.edu>
> Cc: aspect-devel <aspect-devel at geodynamics.org>
> Subject: Re: [aspect-devel] How to grab the radius of every interior
>         sphere made up from cell faces in the domain and some MPI questions
> Message-ID:
>         <
> CACg+8PRfGDny2B43n2o82t9G2YtbRWF9XJAsTE7vA8yGrPjG_w at mail.gmail.com>
> Content-Type: text/plain; charset="utf-8"
>
> Hi;
>
> To make the problems more explicit, here I attach the formulas. After
> transferring the surface topo, CMB topo, and each depth's density anomaly
> into spherical harmonic coefficients, the geoid anomaly is computed in this
> way:
>
> [image: Inline image 1]
> G is the gravitational constant, g is the gravity, n,m are spherical
> harmonic order and degree. The first item is the integral from core mantle
> boundary r_c to surface r_s of the density anomaly spherical harmonic
> coefficient. The second term is the contribution from surface dynamic topo,
> delta_rho_s is the density contrast between mantle and sea water. The third
> term is the contribution from CMB dynamic topo, delta_rho_c is the density
> contrast throughout the core mantle boundary.
>
> For the numerical integral of the first term, my question is how to fast
> grab the radius of each spherical layer determined by the cells' upper face
> in ASPECT/deal.ii (In 0 initial refinement degree, there are three values
> (3480,4501,6371); in 1 initial refinement degree, there are 5 values...). I
> think my current approach to loop over every cell and remove the same upper
> face radius from the cells at the same depth then put all the radius values
> into a vector is not efficient enough and will also not work in parallel
> computing.
>
> Best,
>
> Shangxin
>
> On Wed, Oct 21, 2015 at 10:41 PM, Shangxin Liu <sxliu at vt.edu> wrote:
>
> > Hi Wolfgang,
> >
> > Let me briefly describe my workflow because computing geoid in spherical
> > harmonic is a bit complicated in detail. There are three parts, surface
> > dynamic topo, CMB dynamic topo, and density anomaly throughout the
> > spherical domain that contribute to geoid. The surface and CMB dynamic
> > topography are easier in programming. I just use the results of surface
> > dynamic topography and then change the if loop condition to search for
> the
> > core mantle boundary cells to compute CMB topo. After getting the
> > topography value of each cell along the boundary, there is a another
> > template in my code to transfer them to the form of spherical harmonic
> > coefficients.
> >
> > As for the density anomaly, The formula needs to first get the density
> > anomaly value per cell of each depth layer (in my code, the number of
> total
> > sphere layers is determined by initial refinement degree, i.e, there are
> 2
> > layers in 0 refinement degree, each sphere layer has 48 density values,
> one
> > value per cell). Then I transfer the density anomaly values for each
> layers
> > into spherical harmonic coefficients, rho_l_m(r), and integral the
> > coefficients from r_CMB to r_surface (in the 0 refinement degree, the
> > numerical integral contains two layers, 3480 to 4501km and 4501 to 6371
> > km). So in this integral, I need to grab the radius value of each
> internal
> > sphere, i.e, the radius of the cell's upper face globally (in simplest
> > refinement degree 0 case there are three radius 3480,4501,and 6371). My
> > current approach is to loop over all the cells and get the location of
> each
> > cell's upper face midpoint, then compute the distance between earth
> center
> > and these midpoint, remove the same distance value of the cell at the
> same
> > depth, finally push_back them to a "R_sphere vector" ((3480,4501,6371)
> in 0
> > refinement degree). But I think this approach is not an elegant way and
> > also will not work in parallel since each processor only holds a
> > sub-domain. So is there a fast way to directly grab the radius of each
> > internal sphere determined by the cell face in deal.ii/ASPECT?
> >
> > My second problem is the comparison between the results of different
> > processor. To be concise, there is a loop through all the cells to find
> > which cell's upper face midpoint is the closest to a fixed point. So in
> > parallel running, each processor returns a local_min_distance value and
> the
> > location of cell face's midpoint(x,y,z) corresponding to this value. I
> need
> > to find the min distance among them and also the location of the cell
> > face's midpoint corresponding to this min_distance. I know in ASPECT the
> > namespace Utilities::MPI::min can return the min value among the results
> of
> > all the active processors, but I also need to grab the location of the
> cell
> > face's midpoint correspond to this min_distance contained by one of the
> > processors. I don't know how to do this. Is there any ways to loop
> through
> > all the MPI processors to find both the min_distance and the
> corresponding
> > cell face's midpoint? I wonder whether deal.ii/ASPECT has some commands
> > like "variable_name(MPI_id)" to return the result of each processor.
> >
> > Best,
> > Shangxin
> >
> > On Wed, Oct 21, 2015 at 5:34 PM, Wolfgang Bangerth <bangerth at tamu.edu>
> > wrote:
> >
> >>
> >> Shangxin,
> >> I have difficulty following your questions, and I believe it is because
> >> you have already picked a specific algorithm and you are asking how to
> >> implement it. But we know neither the algorithm, nor what the problem is
> >> you are trying to solve.
> >>
> >> Can you give us an overview *in formulas* of what you are actually
> trying
> >> to compute, along with an overview of *how* you try to compute it? Maybe
> >> that already helps us give you suggestions how to avoid some of the
> >> problems you see.
> >>
> >> Best
> >>  W.
> >>
> >>
> >> On 10/21/2015 03:23 PM, Shangxin Liu wrote:
> >>
> >>> Hi;
> >>>
> >>> I'm now making a new plugin into the postprocess directory to compute
> >>> geoid anomaly. The code has been finished but it needs some further
> >>> modification to run in parallel. To realize this, I run into these two
> >>> problems:
> >>>
> >>> 1. The geoid computation needs the integral of the density anomaly from
> >>> CMB to surface. For now, I do the integral based on cells. So to get
> all
> >>> the radii of the internal spheres, the code loops over all of the cells
> >>> and get the radius of each cell's upper face and push_back the radius
> >>> into a vector after removing the same values from the cell faces at the
> >>> same depth. But I suppose this solution is not so elegant and also will
> >>> not work in the multiple MPI processors since each processor will only
> >>> cover part of the domain. So is there a fast way in deal.ii/ASPECT to
> >>> directly grab all the radii of the spheres determined by the cell faces
> >>> within the domain?
> >>>
> >>> 2. In the spherical harmonic coefficients computation, I need to
> >>> evaluate the values at each grid points I created prepared for sphere
> >>> integral from the output values at the center of each cell face. My
> >>> current solution is to loop over all the sphere cells and find the
> >>> closest cell face midpoint to the grid point (nearest point
> >>> interpolation). But since each processor will only cover part of the
> >>> domain so I also need to compare among all the processors to find the
> >>> closest point and the corresponding density anomaly/topo values at that
> >>> point. I notice that there is a Utilities::MPI::max/min command can
> >>> compare among different processors. But my purpose is not only to get
> >>> the closet distance, but also need to know which processor provide this
> >>> closest cell face midpoint and the corresponding location of this
> closet
> >>> cell midpoint hence to get the density anomaly/topo there. Are there
> any
> >>> MPI commands in deal.ii/ASPECT to realize this?
> >>>
> >>> Best,
> >>> Shangxin
> >>>
> >>
> >>
> >> --
> >> ------------------------------------------------------------------------
> >> Wolfgang Bangerth               email:
> bangerth at math.tamu.edu
> >>                                 www:
> http://www.math.tamu.edu/~bangerth/
> >>
> >>
> >
> -------------- next part --------------
> An HTML attachment was scrubbed...
> URL: <
> http://lists.geodynamics.org/pipermail/aspect-devel/attachments/20151022/21fce3d3/attachment.html
> >
> -------------- next part --------------
> A non-text attachment was scrubbed...
> Name: Screen Shot 2015-10-22 at 12.12.18 PM.png
> Type: image/png
> Size: 16938 bytes
> Desc: not available
> URL: <
> http://lists.geodynamics.org/pipermail/aspect-devel/attachments/20151022/21fce3d3/attachment.png
> >
>
> ------------------------------
>
> Subject: Digest Footer
>
> _______________________________________________
> Aspect-devel mailing list
> Aspect-devel at geodynamics.org
> http://lists.geodynamics.org/cgi-bin/mailman/listinfo/aspect-devel
>
> ------------------------------
>
> End of Aspect-devel Digest, Vol 47, Issue 3
> *******************************************
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.geodynamics.org/pipermail/aspect-devel/attachments/20151022/778d1629/attachment-0001.html>