Artifacts in dynamic topography postprocessor

[ian-r-rose]

In connection with developing a geoid postprocessor, I have been looking at the dynamic topography postprocessor. When I put in very simple initial conditions,a degree-two temperature perturbation in the mid-mantle, I would want to see a corresponding degree-two dynamic topography pattern. Visually, the velocity field seems reasonable, but the dynamic topography (which is related to the gradients in the velocity field) does not. Instead I see large artifacts which look related to the initial dodecahedron used to generate the starting mesh.

Any ideas as to what could be causing this (especially from @jaustermann )? This example is with a fairly coarse mesh, but I feel it should be more than enough to give reasonable results for such a long-wavelength feature.

Here is a full prm which reproduces this:

set Dimension                              = 3
set Use years in output instead of seconds = true
set End time                               = 0.
set Output directory                       = output
set Use direct solver for Stokes system    = false
set Linear solver tolerance                = 1.e-9

subsection Material model
  set Model name = simple
  subsection Simple model
    set Thermal expansion coefficient = 1e-4
    set Viscosity                     = 1e22
  end
end

subsection Geometry model
  set Model name = spherical shell

  subsection Spherical shell
    set Inner radius  = 3.5e6
    set Outer radius  = 6.371e6
  end
end

subsection Model settings
  set Zero velocity boundary indicators       = 
  set Tangential velocity boundary indicators = outer, inner
  set Prescribed velocity boundary indicators =
  set Remove nullspace                        = angular momentum
  set Include shear heating                   = false
end

subsection Boundary temperature model
  set Model name = spherical constant
  subsection Spherical constant
    set Inner temperature = 100.
    set Outer temperature = 100.
  end
end

subsection Initial conditions
  set Model name = harmonic perturbation
  subsection Harmonic perturbation
    set Magnitude = 10.
    set Reference temperature = 100.
    set Vertical wave number      = 1
    set Lateral wave number one = 2
    set Lateral wave number two = 2
  end
end

subsection Gravity model
  set Model name = radial constant
  subsection Radial constant
    set Magnitude = 9.8
  end
end

subsection Mesh refinement
  set Initial global refinement          = 2
  set Initial adaptive refinement        = 0
  set Strategy                           = temperature
  set Time steps between mesh refinement = 0
end

subsection Postprocess
  set List of postprocessors = visualization, velocity statistics, temperature statistics, dynamic topography
  subsection Dynamic Topography
    set Subtract mean of dynamic topography = true
  end
  subsection Visualization
    set List of output variables = dynamic topography
    set Output format                 = vtu
    set Time between graphical output = 0e6
    set Number of grouped files       = 1
    subsection Dynamic Topography
      set Subtract mean of dynamic topography = true
    end
  end
end

[gassmoeller]

I see the same pattern in the output of the geoid postprocessor in #648 (reasonable, since it uses the same calculation for the topography) for coarse grids. It gets smaller with higher resolution (from 3-4 it is smaller than the real perturbation for our testcase). I thought it was a flaw in the geoid postprocessor implementation, but should have taken a look at the dynamic topography, I guess.

Can this be related to some kind of distortion of the cells at the edges of the initial mesh? Since the topography postprocessor uses a midpoint quadrature rule to evaluate the solution and then uses the dynamic pressure at this position to calculate the topography, a distorted cell would have a different midpoint depth and therefore a different pressure.

[jaustermann]

I've seen this pattern in the dynamic topography and discussed it with Wolfgang and Timo and also Shangxin at the hackathon this year. We came to the conclusion that it has to do with the way the sphere is partitioned. Most gridpoints have eight neighboring elements, however some few only have six neighbouring elements. In the screenshot you attached Ian you can see that most points have 4 neighboring elements (on the surface) but 8 specific points have only 3 neighboring elements. This asymmetry in the mesh introduces differences in the numerical accuracy and therefore differences in the solution. This problem gets better when using a mesh with manifolds (as you might already do) and gets smaller with higher resolution. I have a feeling it also depends on the perturbation you impose but haven't really come up with a way to quantify that.

To your point Rene, regarding the depth at which the solution is evaluated: I don't think that is the problem because if you look at the velocity field within the sphere you do see this structure as well (at least if I remember that correctly from when I did those runs).

Thanks for looking into this and let me know if this explanation made sense or you have ideas on how to improve on that.

Cheers, Jacky

On Oct 23, 2015, at 6:08 PM, Rene Gassmöller notifications@github.com<mailto:notifications@github.com> wrote:

I see the same pattern in the output of the geoid postprocessor in #648https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_geodynamics_aspect_pull_648&d=CwMCaQ&c=WO-RGvefibhHBZq3fL85hQ&r=j5AqZvMsoErn2L-vpbdTErPRtyT4BhzQUwKzsenbbTc&m=kb2ZrfJVoRiEZUv0PpFT2vtc_jenw9ai1Kg9EEVjA9I&s=Pv_yavuJKp8uqESm0RXrUM-b_CQEvbAUtpmuLeiE_qQ&e= (reasonable, since it uses the same calculation for the topography) for coarse grids. It gets smaller with higher resolution (from 3-4 it is smaller than the real perturbation for our testcase). I thought it was a flaw in the geoid postprocessor implementation, but should have taken a look at the dynamic topography, I guess.

Can this be related to some kind of distortion of the cells at the edges of the initial mesh? Since the topography postprocessor uses a midpoint quadrature rule to evaluate the solution and then uses the dynamic pressure at this position to calculate the topography, a distorted cell would have a different midpoint depth and therefore a different pressure.

— Reply to this email directly or view it on GitHubhttps://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_geodynamics_aspect_issues_651-23issuecomment-2D150704209&d=CwMCaQ&c=WO-RGvefibhHBZq3fL85hQ&r=j5AqZvMsoErn2L-vpbdTErPRtyT4BhzQUwKzsenbbTc&m=kb2ZrfJVoRiEZUv0PpFT2vtc_jenw9ai1Kg9EEVjA9I&s=NwRBty5gxtJOdgjMTeK5us4RogOaCdPXafinSlOptIE&e=.

[bangerth]

Yes, it boils down to the fact that the symmetry of the mesh is broken at some vertices. In finite element theory, you can often show that in places where four quadrilaterals come together in a regular way, the error happens to be one order better than what one would expect, whereas at other places it simply obeys the expected convergence order. You observe the equivalent effect in 3d: the error just happens to be large at the locations where only 3 hexes come together (and at the vertices below where 6 hexes come together). The effect disappears/becomes less noticeable/the error becomes smaller than the solution as you refine the mesh.

I think that there's really nothing very much you can do about other than use a finer mesh :-(

[ian-r-rose]

More thoughts on this problem:

I have been doing a bit of experimentation with the quadrature used to average cell stress at the surface, and it makes a big difference. The following are three images of dynamic topography due to a Y20 perturbation at depth. They are, respectively, using midpoint, trapezoidal, and Gauss quadrature rules. The current implementation uses midpoint.

Gauss quadrature for the averaging seems to be significantly better to my eyes.

[jaustermann]

Hi Ian,

I looked into this issue of where to evaluate the stress a while ago when we did the plume in a box comparison runs between stress based, free surface and sticky air. I compared the following three approaches: evaluating the stress at the midpoint, at the surface and doing a linear extrapolation to the surface from two points within the cell. For the plume example I found that resolution is actually a stronger control than the method of where to evaluate the stress. That means at low resolution the three methods yielded very different results, at medium resolution they yielded very similar results but the topography still changed significantly with increasing resolution. Therefore I figured that if your resolution is good enough it shouldn't make too much of a difference what method you use. Now this might be different for the artifacts on the spherical shell, which I didn't look at. The midpoint, trapezoidal and gauss quadrature rules are different ways to average over the surface cell, would it make sense to use an extrapolation scheme (to the surface face) rather than a whole cell averaging scheme?

Cheers, Jacky

On Nov 3, 2015, at 2:19 PM, Ian Rose wrote:

More thoughts on this problem:

I have been doing a bit of experimentation with the quadrature used to average cell stress at the surface, and it makes a big difference. The following are three images of dynamic topography due to a Y20 perturbation at depth. They are, respectively, using midpoint, trapezoidal, and Gauss quadrature rules. The current implementation uses midpoint.

[midpoint_quadrature]https://urldefense.proofpoint.com/v2/url?u=https-3A__cloud.githubusercontent.com_assets_5728311_10918850_a4a98e20-2D821c-2D11e5-2D8267-2Daadd3505f168.png&d=CwMCaQ&c=WO-RGvefibhHBZq3fL85hQ&r=j5AqZvMsoErn2L-vpbdTErPRtyT4BhzQUwKzsenbbTc&m=rKAjTS75Dsv7_mta6OQuqKuGSoTPZiV-Dc0Upd_jTVk&s=HPdT6KglmXix1HOGQipSi3o_e6bpPU_1fMmYi42w8Dw&e= [trapezoidal_quadrature]https://urldefense.proofpoint.com/v2/url?u=https-3A__cloud.githubusercontent.com_assets_5728311_10918853_a8d5ccd4-2D821c-2D11e5-2D8112-2Df0e45a10c009.png&d=CwMCaQ&c=WO-RGvefibhHBZq3fL85hQ&r=j5AqZvMsoErn2L-vpbdTErPRtyT4BhzQUwKzsenbbTc&m=rKAjTS75Dsv7_mta6OQuqKuGSoTPZiV-Dc0Upd_jTVk&s=K3YNcmIfjQ8xCm52HGqOnTuZSbi5O1oFgXn7SWlCUWM&e= [gauss_quadrature]https://urldefense.proofpoint.com/v2/url?u=https-3A__cloud.githubusercontent.com_assets_5728311_10918856_abf87a6a-2D821c-2D11e5-2D9adb-2D33046a66d51c.png&d=CwMCaQ&c=WO-RGvefibhHBZq3fL85hQ&r=j5AqZvMsoErn2L-vpbdTErPRtyT4BhzQUwKzsenbbTc&m=rKAjTS75Dsv7_mta6OQuqKuGSoTPZiV-Dc0Upd_jTVk&s=BYYyHOFSzvbBh9aAPWMCp80L1Cmjlx2-CotYMMJEFn8&e=

Gauss quadrature for the averaging seems to be significantly better to my eyes.

— Reply to this email directly or view it on GitHubhttps://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_geodynamics_aspect_issues_651-23issuecomment-2D153460012&d=CwMCaQ&c=WO-RGvefibhHBZq3fL85hQ&r=j5AqZvMsoErn2L-vpbdTErPRtyT4BhzQUwKzsenbbTc&m=rKAjTS75Dsv7_mta6OQuqKuGSoTPZiV-Dc0Upd_jTVk&s=rsHYDydY10gnpl3m3Y-BHDFiY1j6OxbYVizw-zO_xLA&e=.

[ian-r-rose]

Thanks for the info, Jacky. One thing that I have been looking at for the geoid has been to do a sort of hybrid cell surface/volume approach. The idea is to

1. Do averages of the pressure part on the surface. Pressure is a primary variable and hopefully the loss of accuracy near the boundary is not too bad, since we do not need to differentiate it. This also avoids the need to subtract off a mean topography.
2. Do averages of the viscous stress part over the volume of the surface cells. This would attempt to avoid problems with overshoots of the polynomial shape functions that you and Wolfgang have dealt with.

Thoughts about this approach?

[ian-r-rose]

More proof-by-pretty-pictures, Midpoint quadrature vs Gaussian quadrature for evaluating cell averages for surface stress:

This is with four levels of uniform refinement and a Y22 density field. All parameters are the same except for the quadrature rule. I am starting to convince myself that this is a worthwhile thing to do.

[bangerth]

These results certainly suggest that the midpoint rule is not a good choice :-(

From a mathematical perspective, one would of course want to evaluate the stress at the surface, not in the cell. But that leads to poor results with quadratic finite elements on coarse meshes because of an undershoot in the displacement field at the boundary. So one has to average the stress field on a cell, as @jaustermann mentioned. The midpoint rule was the simplest way to do this, but you may well have discovered that it pays to try being a bit more accurate.

[jaustermann]

Yep, that does look pretty convincing, thanks for looking into this Ian! I also realized that when I mentioned to Rene earlier that 'I don't think that the averaging is the problem because if you look at the velocity field within the sphere you do see this structure as well (at least if I remember that correctly from when I did those runs).' it referred to runs before I did the sphere with the manifold mesh. At that point the asymmetry in the mesh made these artifacts very pronounced and not only in the dynamic topography but also the velocity field itself. This certainly got much better (vanished?) with the manifold meshing.

Concerning your approach: That sounds good with the only caveat of course that you're using the pressure and stress from different depths and I don't have a good intuition for whether that's fine or not. I see your point about using the pressure from the surface which circumvents the issue of subtracting the mean topography. Doing a few test runs with cell averages vs. surface averages for the pressure at different resolution would probably give some insight (you might have already done that) but my guess would be that whether you evaluate the pressure as a cell average or the surface should not matter much once you get to a decent resolution (> 4).

In any case, great catch about the gaussian quadrature. We should adapt that in the DT postprocessor(s). Let me know if you want me to look at that at any point or if you just want to copy that over from the geoid postprocessor once you've decided what exactly is the best approach.

Cheers, Jacky

On Nov 4, 2015, at 7:49 AM, Wolfgang Bangerth notifications@github.com<mailto:notifications@github.com> wrote:

These results certainly suggest that the midpoint rule is not a good choice :-(

From a mathematical perspective, one would of course want to evaluate the stress at the surface, not in the cell. But that leads to poor results with quadratic finite elements on coarse meshes because of an undershoot in the displacement field at the boundary. So one has to average the stress field on a cell, as @jaustermannhttps://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_jaustermann&d=CwMCaQ&c=WO-RGvefibhHBZq3fL85hQ&r=j5AqZvMsoErn2L-vpbdTErPRtyT4BhzQUwKzsenbbTc&m=DdDsklltsyOpzT0GuakhQ5tN8g4rGGGy-A5dZElxFBU&s=PK0Zc6WXT3WR_0EkFmWtW9dgmJoueBMjbwlS2U5Rbb0&e= mentioned. The midpoint rule was the simplest way to do this, but you may well have discovered that it pays to try being a bit more accurate.

— Reply to this email directly or view it on GitHubhttps://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_geodynamics_aspect_issues_651-23issuecomment-2D153711654&d=CwMCaQ&c=WO-RGvefibhHBZq3fL85hQ&r=j5AqZvMsoErn2L-vpbdTErPRtyT4BhzQUwKzsenbbTc&m=DdDsklltsyOpzT0GuakhQ5tN8g4rGGGy-A5dZElxFBU&s=GNm_sVIV2IvCEbmHHTrvCZ2ea7cCEhEZDGZtk6TWopw&e=.

[gassmoeller]

Should be closed by #662.

[ian-r-rose]

I actually have one more idea I want to investigate this week, so if you don't mind I'll keep this open for a few more days.

[gassmoeller]

It has been more than a few days since the last comment here :smile:, also I can not longer reproduce the problems with the above parameters, so I guess it has been resolved successfully. Feel free to reopen if anyone encounters the problem again.