Get rid of grid stretching in the overlap zone

[hzfmer]

For stability purpose, we would like to make sure the distortion of the fine grids terminates above the overlap zone, which means the grids inside and below the shallowest overlap zone should be Cartesian instead of curvilinear.

A possible solution is to add an 1D array to mark whether a certain layer is above or below the overlap zone.

[ooreilly]

@hzfmer since you assigned yourself to this task are you working on it? I have been working on it but I haven't solved it yet (been out for this week, but plan to return to it on Monday).

[hzfmer]

Hi Ossian, I created the issue because Te-yang and I saw some instabilities when the model includes multiple blocks and topography. It seems to be of high priority. I am not sure if you are able to work on this, so I assigned myself, but it would be very helpful if you could work on this as well. The curvilinear finite difference part looks very complicated and I had a hard time understanding it. I could also provide any help you need if you decide to work on this issue. Cheers, Zhifeng

…

On Nov 22, 2019, 10:01 -0800, ooreilly ***@***.***>, wrote: @hzfmer since you assigned yourself to this task are you working on it? I have been working on it but I haven't solved it yet (been out for this week, but plan to return to it on Monday). — You are receiving this because you were assigned. Reply to this email directly, view it on GitHub, or unsubscribe.

[ooreilly]

@hzfmer understood. Has it been confirmed that the instabilities originate from the overlap zone?

Yes, the kernels are quite complicated due to the mix of interpolation, differentiation, and metric coefficients that are needed at different grid locations. I started adding all the masks and I put together a test, and a simulation test that went unstable. I need to prototype more to understand how the mapping needs to be designed.

[hzfmer]

It is very likely due to the overlap zone.

We have done multiple different tests:
Stable:

• One block, homogeneous, without topography
• One block, homogeneous, with flat topography
• Two/three blocks, homogeneous, without topography
• Two/three blocks, heterogeneous, without topography

Unstable:

• Two blocks, homogeneous, with topography
• Two blocks, heterogeneous, with topography

We tried point source and multiple sources, the stability was not changed. I believe you tested one heterogeneous block with Gaussian topography before, which was stable.

We are testing two homogeneous blocks with flat topography. But up to now, we suspect that the overlap zone is the source of instability.

[ooreilly]

Very good. Please keep me posted.

I didn't test heterogeneous properties with a Gaussian. I used an incline plane for the topography function. But yes, it was just one block.

This reminds me, did we ever test acoustic material properties? @hzfmer if not, you could please create an issue for it so that we don't forgot to do so in the future? Thanks.

[hzfmer]

@ooreilly I haven't tested acoustic materials, but @deanrockit might have done some tests. I will check with him and create an issue later.

[hzfmer]

I did some more tests:

Stable:

• One block, homogeneous, without topography
• One block, homogeneous, with flat topography
• One block, homogeneous, with Gaussian topography
• Two/three blocks, homogeneous, without topography
• Two/three blocks, heterogeneous, without topography

Unstable:

• One block, homogeneous, with real topography
• Two blocks, homogeneous, with Gaussian/real topography
• Two blocks, heterogeneous, with Gaussian/real topography

It looks that even for one homogeneous block, it can be unstable when having real topography (i.e. not an analytical function). We should check if the padding caused instability. In my tests, the padding zone comes from real topography as well.

[ooreilly]

@hzfmer

1. What time step are you using? The maximum stable time step depends on the length scale associated with the smallest cell size.
2. Are you using frequency dependent Q? I have only done stability proofs and tests for a purely elastic test case.
3. By padding, do you mean sponge layer? we have seen instabilities coming from having the sponge layer in the past ( I saw that when testing using an incline plane)

[deanrockit]

@ooreilly Since you mentioned, I have a question about how grid spacing in the curvilinear grid system.
In case of positive topographic height (above sea level), the dh of the mapped grids should be larger or smaller than the unmapped dh?
The second question is how should we define the grid spacing in the curvilinear grid system? Should we measure the distance to the neighboring grid in all three directions or just along z?

[hzfmer]

I am out of town today, but I will be beck tonight. 1. I have checked that the maximum cmax * dt / dx = 0.22, which is smaller than the CFL. 2. It is feequemcy dependent q, i will test constant q later. 3. By padding i mean the extra 8 layers surrounding topography.

…

-Zhifeng Hu

On Nov 23, 2019, 11:59 -0800, ooreilly ***@***.***>, wrote: @hzfmer 1. What time step are you using? The maximum stable time step depends on the length scale associated with the smallest cell size. 2. Are you using frequency dependent Q? I have only done stability proofs and tests for a purely elastic test case. 3. By padding, do you mean sponge layer? we have seen instabilities coming from having the sponge layer in the past ( I saw that when testing using an incline plane) — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.

[ooreilly]

@deanrockit If the topography function is positive, the grid spacing in the vertical direction of each cell is stretched, and therefore this length scale increases compared to the flat case. In contrast, if the topography function is negative, the grid spacing in the vertical direction of each cell is squashed, and therefore the length scale decreases compared to the flat case. However, the grid spacing also changes in the other two directions as well.

You need to determine the non-uniform grid spacing in each direction of each cell using the mapping function. The shortest such length scale would be necessary for stability. However, note that I said "necessary" not sufficient. The scheme might still be unstable due to stiffness caused by discretization at the boundary. It is difficult to deduce how the time step might need to be altered because of this. Usually, I perform an eigenvalue analysis to study this behavior. In practice, that is quite infeasible to do.

@hzfmer I would be surprised if the instability is due to CFL if you are running at cmax *dt / dx = 0.22. Are you running with a sponge layer?

[deanrockit]

@ooreilly In practice, this will be very important and should be the first thing to address when setting up a simulation. Since you mentioned it's not feasible to do the eigenvalue analysis you normally do, would you recommend any other method that can assure stability of the simulation?

I tested a simulation with a topo field characterized by white noise with values ranging from 0 to 1000m. All the values are positive. With positive topography values, grids should be stretched and therefore dx increased. CFL should be fine in this case, but it went unstable after 100 time steps.
I'm still investigating the starting point of instability. Meanwhile, would you give us a heads on up on what else we should look out when performing the curvilinear calculations?

[ooreilly]

@deanrockit Yes, we should be able to estimate it from the metric coefficients (derivatives of the topography function). However, there is more to take into consideration.

For a reliable simulation, we need both stability, and accuracy. My recommendation would be to do perform some systematic experiments to establish how time step, and grid spacing depends on the frequency content of the topography function. But before doing any such experiments and discussing how to set them up, it would be good ensure that the problems you are reporting are not due to some bug, or incorrect configuration.

The metrics coefficients are computed discretely by using fourth order difference approximations. It is the approximation of the first derivative that is used to numerically compute the stretched grid spacings. Since white noise is not differentiable, the metric coefficients cannot be accurately computed. As a result, you will most likely be forced to take a very tiny time step to obtain a stable simulation. In any case, it will not matter much because the simulation will not be accurate. Much like that you need have a certain number of grid points per wavelength to resolve the wavefield on uniform grid, there is now an additional criterion that depends on the length scales of the topography.
Hence, some filtering might be necessary to make sure that the topography function is differentiable.

@hzfmer @deanrockit We might want an additional tool that shows you what the metric coefficients look like. This tool would output 2D "images" that correspond to taking derivatives in the x- and y-directions of the topography function. These derivatives are computed by topography/metrics/kernel.c (runs on the CPU) and tested by tests/metrics/test_metrics.c.

[deanrockit]

@ooreilly I'm pretty sure that the topography is configured correctly since I've tested flat topo of 1000m height, and it produced the same result as the non-topo case where the source was offset downward by also 1000m.
To see what happens when it's non-flat, I simply replaced the z values in the flat topo file by random numbers. Random numbers are severely fluctuating throughout the space, so what you just said about random field being non-differentiable makes a lot of sense.
However, there is something I couldn't understand. Since real topography should be considered non-differentiable as well, did you apply any sort of filter to the topography field in your San Jacinto calculations that you showed at SCEC. Topography should be less complicated than random field, which is much less demanding than the random case in terms of dt. Any thoughts?

[ooreilly]

@deanrockit For the San Jancinto test, I didn't any apply any filtering. But it is possible that this topography is much smoother than what you are testing with. That said, I didn't run this test with frequency dependent Q enabled. Are you using frequency dependent Q, or are you running in purely elastic mode? Also, are you using a sponge layer?

I believe that the gmt grd tool is setup to use cubic spline interpolation to query the elevation map model at the grid points of the finite difference grid. I don't know if the original model had been pre-processed in any way before this query step.

[deanrockit]

@ooreilly I ran the case with sponge layer and Q(f), since I would like to know what would happen in a realistic setting. Like you said, it's a great idea to conduct systematic experiments to investigate the stability and accuracy condition and as a function of topographic length scale.
Let me first try some Gaussian hill models first. Start small.

[ooreilly]

@deanrockit I understand that you would like to test the code in a realistic setting. However, we do not fully understand how all the features work together. Therefore, I would advocate a systematic procedure in which each feature is added in a step by step fashion. I'm glad to hear that you agree with me on this point. But before we start any extensive investigation that will provide insight into how to choose the time step, and frequency content (in space) of the topography, I would like to get a sense of what causes the instability in your current simulation (with realistic topography).

Here is how I would proceed:

• Disable the sponge layer, disable frequency dependent Q, run with a single block.
  • If unstable, reduce the time step by a factor of 2 to check that the instability is not caused by stiffness. Keep adjusting the time step if the instability persists. If you are unable to obtain a stable simulation after maybe 2 - 3 of these adjustments, then it is most likely not an instability caused by stiffness, and is then a sign that something is violating energy conservation.
  • If still unstable, we need to re-conduct an energy stability test. More details below.
  • If stable, proceed to next step.
• Enable frequency dependent Q Check stability, and adjust time step if unstable, as before.
  • If unstable. This could be an indication that the stability does not hold for frequency dependent Q. I then need to do a stability analysis for this modification and see if there is some way to fix it. More details below.
  • If stable, proceed to next step
• Enable the sponge layer, and repeat stability test.
  • If unstable, make sure that the topography is flat in the sponge layer region. Make sure to maintain the height of the topography in the inner region so that you don't cause an abrupt jump between the inner region and the sponge layer.

We don't have any stability proofs for topography with frequency dependent Q, nor sponge layer. Simulations with an incline plane have shown that the sponge layer and topography can induce instability. No stability issues are known for when frequency dependent Q has been enabled, but that doesn't mean that there cannot be any as far as I know. If you find issues, please notify me and I will attempt to see if I can do stability analysis for this part.

It has also been a while since I ran the energy stability test. This test is quite difficult to perform, and I didn't perform it after reimplementing topography into the updated frequency dependent Q kernels in AWP. Instead of focusing on ensuring that the curvilinear grid does not stretch within the overlap zone, I'm considering reintroducing this test as my top priority.

If you want to experiment with the Gaussian hill, I want to remind you that we did Gaussian hill tests earlier that were compared against Edge for a small test-case. We also did a large-scale test to study dispersion effects.
https://github.com/SCECcode/awp-benchmarks/tree/master/tests/topography/gaussian_hill
https://github.com/SCECcode/awp-benchmarks/tree/master/tests/topography/gaussian_hill_large
It's been a while since I ran these tests so it is possible that I need to go in update them.

One more thing, would it be possible for you to share your tests to the awp-benchmarks repository so that you, @hzfmer, and I can review it together?

[ooreilly]

I decided to investigate the CFL-type stability condition a bit more. It turns out that I was wrong. I apologize for the mistake.

The stability condition doesn't quite dependent on the smallest grid cell. Using some rather hand-wavy arguments, I find that the time step is inversely proportional to the Euclidean norm of the resultant contravariant basis vector. See the notes that I posted in the awp-topo repository: https://github.com/SCECcode/awp-topo, notes/stability-analysis (you need to be on the stability-analysis branch until the PR has been merged).

[hzfmer]

@ooreilly
It sounds great!

I have found that the sponge layer and frequency-dependent Q are probably not causing instabilities, but dt can be sensitive. Also, when I applied a low pass filter to the real topography, an unstable case can become stable.

In a word, we need look more closely at dt and differentiability of the topography, while sponge layers and Q are okay for now.
I will wrap up my tests when we are clear about how to determine dt and filter the topography.

[ooreilly]

@hzfmer I'm pretty sure that the sponge layer can cause instabilities that are unrelated to time step size. At least that is what I observed when I tested with an incline plane. I was unable to get a stable result by adjusting the time step. I had to disable the sponge layer to resolve the instability.

Now, it might be that for the actual topographies that you use that the sponge layer Q work. If so, that's great! The time step restrictions you are observing sound consistent with the estimate I derived. If you look at the notes you can see that if f_x and f_y (derivatives of the topography function) are large then the time step needs to decrease.

[deanrockit]

@ooreilly It's great that you can derive the stability condition from theoretical point of view instead of empirically guessing. If I understand correctly, in case of highly variable topography where f_x and f_y are large we may want to apply some sort of smoothing operator to avoid using very small dt. Would you suggest this approach?

[ooreilly]

@deanrockit Yes, smoothing should help managing the time step. Smoothing should also help with delivering accurate simulations that are not contaminated with high frequency garbage due to having unresolved features on the grid. I don't understand this resolution criterion yet. Before investigating that further, I would like to put together test and experiments for the time step estimate as previously discussed.