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

[Ian Rose]

That is almost what I am describing, but I am also questioning whether it
is necessary to do the binning in radius. A general multipole expansion is
already written in terms of an integral over volume. The traditional
approach of splitting that up into radial slices is nice for a spherical
harmonic transform, but I do not think it is strictly required. That is to
say, the formula Wolfgang and Shangxin wrote have pulled the radius out of
the integral, and we could just keep it in there.

The downside of this (as Rene pointed out) is that extra work would be
required to get the values of the potential at points interior to the
domain. The radially binned approach has a very natural (if expensive) way
to do that. But if we are only interested in the potential at the surface
and the CMB, then this would be enough.

That being said, if we really wanted the internal potential, it would not
be too expensive to do a basic finite-element Poisson solve for the
interior, since we have the appropriate boundary conditions after the
expansion.

On Sun, Oct 25, 2015 at 9:01 AM, Wolfgang Bangerth <bangerth at tamu.edu>
wrote:

>
> Shangxin & all,
> I am also a bit late to the discussion, but let me give me $0.02 anyway:
>
> I think you started this by finding an algorithm that is well suited to
> one kind of numerical scheme (namely uniform meshes in spherical
> coordinates) and then trying to fit that into a code that implements a
> different method. This sometimes works, but it often leads to algorithms
> that are unnecessarily complicated and inefficient.
>
> Instead, let's start from what finite element codes do really well on
> general unstructured meshes (such as the ones used by ASPECT): computing
> volume integrals. In your case, if you subdivide the Earth into layers,
> then expanding the density into spherical harmonics requires you to compute
> integrals of the form
>
>   r_{lmi} =
>     \int_{r_i}^{r_{i+1}} \int \int  drho(r,phi,theta) Y_{lm} dtheta dphi dr
>
> which you can quite easily implement as follows:
> {
> for (cell=...)
>   if (cell->is_locally_owned())
>     for (q_point=...)
>        {
>           Point<dim> x = fe_values.quadrature_point(q_point);
>           double r = x.norm();
>           int i = ...find index of layer for radius r...;
>
>           rho = ...get rho at x;
>           drho = rho - rho_reference(r);
>
>           for (l=...)
>             for (m=...)
>               r[l][m][i] += drho * Y(l,m,x) * fe_values.JxW(q_point);
>
> Y = Utilities::MPI::sum(Y, mpi_communicator);
> }
>
>
> I believe this is in essence what Ian is describing.
>
> Now the question is how you get your layer subdivision r_i ... r_{i+1}.
> You tried to derive it from the mesh itself, but that's a brittle approach
> that, as you noticed, is not well suites to parallel computing or adaptive
> meshes. My suggestion would be to simply specify a subdivision yourself,
> either equidistant or otherwise, with a number of layers provided by the
> user.
>
> I hope this makes sense! Best
>  Wolfgang
>
>
> --
> ------------------------------------------------------------------------
> 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
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.geodynamics.org/pipermail/aspect-devel/attachments/20151025/80a13d72/attachment.html>