Negative Jacobian

[kuangdai]

I think there may be a bug in the geometric models. Marjolaine has been unable to get hers to run for several weeks now, and I am trying to run a simulation for Kristiina, and it does not work either. It's always a negative Jacobian error. I suspect it's something to do with coordinate transformations in regional meshes, but I can't say for certain. It is equally possible that I am doing something wrong. But there are two people depending on me to figure this out asap and I am stuck, so I'd be very grateful for your help!

I've been trying to figure out what is happening and I've attached two figures I've created in the process: The bathymetry_map shows the geometric model. The layer is at 2000m depth and the range is [0 12170]. The circle shows where the mesh should be (source at lat 37.5, lon -16.5, curved mesh with max colatitude 3.2 degrees). I marked data points in red where the bathymetry goes above 2000m, i.e. where the islands are - it should be a safe distance from the mesh. geometric_gone_wrong is a scatter plot of some data I put out in Quad::setupGLL. This is for the nr=0 slice (nr is set to constant 100, bound by inplane): x-axis: theta of the GLL points obtained via geodesy::sz2rtheta(getPointSZ(), false) z-axis: undulation (see below) + r of the GLL points obtained via geodesy::sz2rtheta(getPointSZ(), false) colormap: undulation of the GLL points at nr=0, obtained via mUndulation->getPointwise()ipnt This scatter plot looks very wrong. Do you have any idea what might be the issue?

Thanks a lot! Claudia

[kuangdai]

The way you made the second plot is not 100% right because you only used the zeroth order coefficient.

Please upload the complete input folder.

[tnissen]

And here is a description from Marjolaine on this topic, with input files and sketch attached. For her, it is really quite urgent as she has to submit her thesis in a few weeks. So either we get it to run in some fashion very soon, or not at all. One compromise could be to just undulate the topography and leave deeper layers flat, but it seems that even that doesnt work right now.

Model space:

Model space covers everything south of -60S. Base model consists of PREM, with lower, mid and upper crust above (from crust 1.0) and then an ice layer above that. I have five geometric models to undulate the top five interfaces, which are: atm_ice, ice_uc, uc_mc, mc_lc, lc_mant. The undulation data for the top two of these interfaces are interpolated from ‘BedMachine Antarctica v2’; the bottom three interfaces from crust 1.0.

It’s worth noting that these interfaces are sometimes quite close together (~1km) and not only do the ranges of the geometric models overlap, but a lot of the interfaces also lie in the range of neighbouring geometric models. A sketch of the situation is attached (this just shows the base model interfaces and the range of undulations for each geometric model).

Errors:

I’ve been trying to run simulations with this full-complexity model but have had a number of issues. When I activated all 5 geometric models, I was getting a ‘negative jacobian or distorted element occurred’ error. Claudia explained that if an interface lies in the range of two geometric models, it’s final relabelling will be a sum of the effects of both, so the issue was likely that I hadn’t removed the effects of neighbouring geometric models (and so the final undulations were no longer within the given range). However, I also get the negative jacobian error when I apply only one geometric model at a time (with nr = 1000 but not when nr ~ 5). This was the case for all of the interfaces, except uc_mc, which started running without an error but blew up at ~25% completion. This interface also happens to be the only one whose geometric model range doesn’t overlap with either of the neighbouring interfaces – surely this is not a coincidence?

Claudia also suggested that it could be an issue arising from the spatial sampling. I tried with a 2s mesh (instead of the 10s one) and still got the negative jacobian error. The BedMachine data has a spacing of 1km at -60S (but is on a regular lat lon grid so resolution increases towards the south pole) but crust 1.0 is clearly lower resolution, so I'm not sure how a spacing issue could be affecting the geometric models of those deeper interfaces as well.

[chaindl]

input.zip Here you go. Thanks Kuangdai!

[kuangdai]

1) In bathy_geo_resampled.nc, we can see that the shape of dz (the undulation data) is 1201x1051, so the data matrix is arranged by lon x lat not lat x lon, and so the input variable data_rank should [1, 0] instead of [0, 1]. The messy plot you made is likely to be caused by this dimension error. Actually if you run the programme in debug mode, Eigen will abort for size incompatibility. There should not be a bug for calculating spherical coordinates -- this part of code is modularised and used everywhere. If there were an error, the code wouldn't work at all.

2) After changing 1), we still get negative Jacobian, because there's another problem. You are using depth_below_solid_surface=false, which means all depth are in reference to the water surface instead of ocean floor. Given that, min_max: [-1000., 12170.] does not make too much sense, as it works mostly like [0, 12170] (-1000 is in the air). Because interface=2000., the maximum possible undulation should be 2000, otherwise that interface will rise above the sea. However, the programme shows that the data ranges between [-4388.67, 3689.65]. Reducing the undulation by changing factor can make it run. 0.7 works.

Please check your data and input.

[chaindl]

Thank you for the suggestion with the data ranks - I forgot about that. But depth_below_solid_surface=false is the intended setting. There is an island in the model but it is outside the modeled domain, as shown in the first figure: All areas with dr>=2000 are highlighted in red and the circle shows the region in range of the mesh. So, the problematic areas should not be sampled.

[kuangdai]

Because the simulation runs with a slightly smaller undulation, we can infer that one or two elements are distorted because of some large/sharp local undulations. You can try

1) smooth the data (e.g., with 2D Gaussian filter) a little bit, because raw data sampled from ETOPO models are very detailed. 2) use a higher frequency mesh and a larger Nr -- the sparser the sampling, the more likely some elements are distorted

This could not be a code problem. Let us close, please.

[kuangdai]

Changing undulation range to min_max: [-1000., 50000.] will solve the issue. Some downside undulation is too large and we can make them vary more smoothly by extend the lower bound a little bit.

[kuangdai]

The attached figures show the model as seen by AxiSEM3D. Nothing seems wrong. I don't know why you had the vertical bands in your plot.

The negative Jacobian is caused by the fact that, while you wanted to use 2000 m as the interface, the first element layer (surface layer) has a thickness of 2350 m. In other words, the interface cut across the lower part of all surface elements. For a downside undulation, the upper part of these elements will be stretched (in great amount) and the lower part squeezed. Such a sharp change caused an overshooting when the element computes the gradient of undulation and led to negative Jacobian.

The solution is simple: I changed the interface location to 2350 m. If you want to use 2000 m, make a factious discontinuity there.

[axelwang]

Hey @kuangdai ! I tried to plot the GLL points like the first figure you show above, but couldn't find where they are. I looked into the exodus files and found they only contain the geometrical nodes. I am able to plot all the elements since I know all the node coordinates but I couldn't find where are the coordinates stored for the GLL points?