Adding SemUCB model to SpecFEM3D_globe #1718
Replies: 33 comments 32 replies
-
I have tried to isolate the effect of the crust by setting |
Beta Was this translation helpful? Give feedback.
-
agree, these flags are quite confusing and the naming is not very intuitive - and their functionality depends on a number of other settings... in general, the idea of the mesher is to take a 1D reference model (defined for inner core, outer core, mantle & crust, like PREM) and overimpose 3D mantle variations (usually given as perturbations wrt the 1D model value), and if possible 3D crustal variations. so,
the other flags are shortly explained in the code file
the use of these flags for crustal models depends if you're using a 1D crust from a 1D reference model, or if you're overimposing a 3D crustal model on top:
therefore, for a 3D crustal model, you should set |
Beta Was this translation helpful? Give feedback.
-
Beta Was this translation helpful? Give feedback.
-
that's correct, the and yes, there are differences related to how the geodetic (lat/lon) coordinates are converted to geocentric (theta/phi) coordinates and vice versa between version 7.0 and 8.1. also, the source positioning has slightly changed. |
Beta Was this translation helpful? Give feedback.
-
Hi Daniel,
and equivalent ones in specfem 7.0 ( it is not quite the same, because at this stage specfem 7.0 still works in cartesian coordinates while specfem 8.1 has already switched to spherical coordinates). When Is there any other place where the crustal model or the 3D model could make a difference? |
Beta Was this translation helpful? Give feedback.
-
Hi @danielpeter ,
8.1:
It seems thus that reading the model is not a problem, and yet the seismograms are substantially different (see earlier plot, I have only added and removed print statements so it shouldn't have changed). |
Beta Was this translation helpful? Give feedback.
-
Hi @danielpeter , |
Beta Was this translation helpful? Give feedback.
-
hi Dorian, sorry for this delay. there are many changes between version 7.0 and 8.1, so i think the mesh layout will have changed slightly from version 7.0 to 8.1 as well. comparing these version's mesh point-by-point, i would expect that they are at slightly different positions. |
Beta Was this translation helpful? Give feedback.
-
Hi Daniel, |
Beta Was this translation helpful? Give feedback.
-
Hi Daniel and Dorian,
Perhaps this could help direct the investigations.
I had first checked Dorian's version 7 against the one we have been running
for years at NERSC, and they match for 3D model SEMUCB with all flags
turned on in Pari_file (attenuation, ellipticity etc..).
I enclose 4 plots for a specific example source-station with a low pass
corner and cut-off at 53 and 40 s, respectively, so rather long period.
Model is SEMUCB with its 3D crust. The source is a deep earthquake (> 600
km) to contrast with the test that Dorian has been doing with a shallow
event.
In all cases, the attenuation model is our 1D model (QL6).
version 7 = "old"
version 8 = "new"
1) compares version 7 against version 8 (the one Dorian is working on) with
all flags turned on in Par_file (attenuation, ellipticity etc..). The
amplitudes in version 8 are much too large, especially for body waves that
include a reflection (much smaller differences for the direct P and S
waves). This might indicate some issue with the 3D crust implementation??.
I ruled out the effect of the flag ATTENUATION_1D_WITH_3D_STORAGE, as it
does not change anything whether it is turned on or off.
2) compares version 8 with and without attenuation (all other flags turned
on). They practically match in amplitude, there is "only" a variable phase
shift.
3) comparing version 7 with and without attenuation.
4) For reference I also include the comparison of 8 and 7 runs without
attenuation, which unfortunately do not match.
Note that surface wave Airy phases match more or less in amplitude, better
than the reflected SS, SP,SSS...
It looks like it might be an issue with the crust implementation, but
haven't found yet where that could happen.
I hope this example might be useful - happy to clarify as needed -
sorry the plots are not all exactly at the same y scale, but I didn't want
to put too many traces on the same plot.
regards
Barbara
…On Wed, Oct 2, 2024 at 12:27 PM Dorian Soergel ***@***.***> wrote:
I just uploaded it to this repo (in the main branch):
https://github.com/dosoe/Specfem3d_globe_7.0_Berkeley/tree/main
—
Reply to this email directly, view it on GitHub
<#1718 (reply in thread)>,
or unsubscribe
<https://github.com/notifications/unsubscribe-auth/BJEVZDOUXVUUHEZLX2R53GTZZRCIXAVCNFSM6AAAAABKSRHCXGVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTAOBSGQ3DEMQ>
.
You are receiving this because you were mentioned.Message ID:
***@***.***>
|
Beta Was this translation helpful? Give feedback.
-
no 3D attenuation in the model
…On Thu, Oct 3, 2024 at 4:30 PM Dorian Soergel ***@***.***> wrote:
So the attenuation works in the 1D case, but not in the 3d case, right? Do
we have 3d attenuation, i e is there a scaling between attenuation and 3d
perturbations? I don't believe so, but you will know better than me.
—
Reply to this email directly, view it on GitHub
<#1718 (reply in thread)>,
or unsubscribe
<https://github.com/notifications/unsubscribe-auth/BJEVZDOLKTRKL3FI7QJ3SGLZZXHRZAVCNFSM6AAAAABKSRHCXGVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTAOBTHAZDKNQ>
.
You are receiving this because you were mentioned.Message ID:
***@***.***>
|
Beta Was this translation helpful? Give feedback.
-
hi Dorian & Barbara, thanks for providing more details! i'm currently comparing this repo: so far, the discrepancies that I see all come from the crustal mesh stretching, source positioning and the geographic to geocentric conversions of the old version versus the new one. I haven't been able though to run it with the full SEMUCB model, as it takes ages to finish the meshing part. There might also be some issue with the It's probably easiest if I start adding your model to the "official" devel version of SPECFEM3D_GLOBE together with a conversion flag to make these comparisons easier for you. for that purpose, I would add your preferred model |
Beta Was this translation helpful? Give feedback.
-
Hi Daniel,
Thanks a lot for your attention.
I have checked that the 3D model that I used in the comparison figures that
I sent is the same for the version 7 and version 8 runs.
I am surprised you say the mesher takes forever to run, because my
experience is similar to Dorian's, it takes no more than 15 mn on Anvil
with NEX = 128 and NPROC = 8.
ABout the ellipticity question: In short, our model (both crust and mantle)
is referred to a spherical earth and we account for geographic to
geocentric coordinate corrections and for ellipticity when building it. So
when we compute synthetics in the SEM we need to first convert geographic
to geocentric both for source and receiver, and ADD the theoretical
ellipticity effect. Perhaps that needs to be done differently in version 7
and 8, but I'd be surprised if the large amplitude discrepancies come from
an error on the ellipticity correction.
Does this answer your question?
The mesh stretching may be a better candidate, how is it implemented
differently in version 8 compared to version 7? I've compared (but perhaps
not enough) at the way it is done but haven't spotted any differences
between the two versions yet...will try to look at it again, also with
Dorian.
regards
Barbara
…On Fri, Oct 4, 2024 at 1:22 PM Dorian Soergel ***@***.***> wrote:
Just to be clear, the model files in
https://github.com/dosoe/specfem3d_globe/tree/devel are the correct ones
for SemUCB
—
Reply to this email directly, view it on GitHub
<#1718 (reply in thread)>,
or unsubscribe
<https://github.com/notifications/unsubscribe-auth/BJEVZDMUTENMCVHXURTOKNLZZ32INAVCNFSM6AAAAABKSRHCXGVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTAOBUHAZTKNA>
.
You are receiving this because you were mentioned.Message ID:
***@***.***>
|
Beta Was this translation helpful? Give feedback.
-
yes, that is one of the differences between new and old version. I'm currently running some tests with the newly added routines in the devel version and let you know about the progress - hopefully this will be fixed soon. to avoid this lat/lon issue for the Berkeley crustal model I'll use the r/theta/phi positions as arguments which will always be geocentric positions. |
Beta Was this translation helpful? Give feedback.
-
Hi Daniel,
Dorian explained very clearly the status of our benchmarks.
Just to clarify: SEMUCB is radially anisotropic across the whole mantle
(and crust). This is how it's implemented in specfem ver. 7 and it seems
Dorian did the necessary to account for that in the "current" version (i.e.
not sure where you see that the anisotropy stops at 220 km?)
regards
Barbara
…On Wed, Oct 9, 2024 at 5:55 PM Dorian Soergel ***@***.***> wrote:
Hi Daniel,
We have set USE_FULL_TISO_MANTLE to .true. and we have modified the
compute_element_tiso_flag subroutine so that elem_is_tiso is always .true.
in the mantle.
I believe using double precision in model_berkeley.f90 as you suggested
earlier is responsible for most of the improvement, that was a really good
catch!
This is the result of a run with all parameters enabled:
B012796B_section_Z_synth_new_3D_semucb_full_alltrue.png (view on web)
<https://github.com/user-attachments/assets/bdfeab20-0bb6-4eda-8c06-3341b2afee48>
As can be seen, the fit is pretty good but not good enough that we could
substitute one for the other , for example for doing tomography
To rule out effects of oceans, rotation and gravity, I have disabled all
three, which doesn't solve anything as you can see here:
B012796B_section_Z_synth_new_3D_semucb_full_topo_ellip_att.png (view on
web)
<https://github.com/user-attachments/assets/76203b15-b9c5-4783-991a-5afb1f0a567d>
I manually disabled the geocentric to geographic correction when reading
the crust (by hardcoding call
thetaphi_2_geographic_latlon_dble(theta,phi,lat,lon,.false.) in
meshfem3D_models_get3Dcrust_val, it's not the best way for long term but
should work for quick testing), the results are so similar that I would
need a magnifying glass to see any difference, in other worlds, it doesn't
solve our problem:
B012796B_section_Z_synth_new_3D_semucb_full_alltrue_spherical_crust.png
(view on web)
<https://github.com/user-attachments/assets/82f610fb-4d04-4a42-adbd-78000ef23d77>
I also thought it could be due to the topography correction (but kept the
geocentric to geographic correction in the crust), but removing the
topography makes the fit much worse. I don't know if it is a hint, maybe a
similar typing error:
B012796B_section_Z_synth_new_3D_semucb_full_ellip_att.png (view on web)
<https://github.com/user-attachments/assets/51bcb381-1e96-4b8d-b092-ff371e7a6412>
So this is the status of our attempts. I'm looking forward to see if your
tests yield any helpful information.
Regards,
Dorian
—
Reply to this email directly, view it on GitHub
<#1718 (reply in thread)>,
or unsubscribe
<https://github.com/notifications/unsubscribe-auth/BJEVZDL4CYMXNBCWJDDF24LZ2XF6TAVCNFSM6AAAAABKSRHCXGVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTAOBZHAZDOMY>
.
You are receiving this because you were mentioned.Message ID:
***@***.***>
|
Beta Was this translation helpful? Give feedback.
-
Dorian & Barbara, you could now run a comparison with the latest devel version of SPECFEM3D_GLOBE. I added the for comparisons, you would turn on the following flag in file
this should give you identical results as with version 7.0 (which is actually source code files from version 6.0, not 7.0). in the updated devel version, there is a corresponding section in the
apart from the |
Beta Was this translation helpful? Give feedback.
-
Hi Daniel,
Cool!
regards
Barbara
…On Thu, Oct 10, 2024 at 3:18 PM Dorian Soergel ***@***.***> wrote:
Hi Daniel,
Thank you, I will run these comparisons end of this week.
Dorian
—
Reply to this email directly, view it on GitHub
<#1718 (reply in thread)>,
or unsubscribe
<https://github.com/notifications/unsubscribe-auth/BJEVZDMI7A2U74FLUQXXASTZ234KZAVCNFSM6AAAAABKSRHCXGVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTAOJQHEYTMNA>
.
You are receiving this because you were mentioned.Message ID:
***@***.***>
|
Beta Was this translation helpful? Give feedback.
-
hi Dorian, these parameters are optional now. you can still put them into the Par_file in your run and they will override the default values when the UCB heaviside is turned on. On 11 Oct 2024, at 02:11, Dorian Soergel ***@***.***> wrote:
Sorry, I thought I spotted something wrong, but I was wrong myself.
—Reply to this email directly, view it on GitHub, or unsubscribe.You are receiving this because you were mentioned.Message ID: ***@***.***>
|
Beta Was this translation helpful? Give feedback.
-
Beta Was this translation helpful? Give feedback.
-
Beta Was this translation helpful? Give feedback.
-
Beta Was this translation helpful? Give feedback.
-
Beta Was this translation helpful? Give feedback.
-
thanks Dorian! I've added some more steps towards the compatibility to the old version implementation of the SEMUCB model. could you test it again with your example, using the new devel version (with PR SPECFEM/specfem3d_globe#851). to have the backward compatibility, re-run the
this should get you closer to the old SEMUCB seismograms, for all possible parameter settings in please note that I had to re-introduce some old SEMUCB hacks that were probably not correct in how things were implemented in the old version. these hacks are only effective when this flag
was comparing the value of rho1D which was non-dimensionalized (i.e., rho / EARTH_RHOAV) to
so even for GLL points below the actual moho interface, it was always using the "lowest" crustal model velocities instead of the mantle velocities. I'm not sure if this was done on purpose, but the "new" version only assigns now crustal velocities to points at or above the moho interface.
and then the new GLOBE version has a few more improvements independent of SEMUCB:
anyway, please compare with the newest devel version and let me know how your updated comparisons look. |
Beta Was this translation helpful? Give feedback.
-
Hi Daniel,
repeated over and over again. The only thing I modified with respect to my previous run is that:
and here is my
|
Beta Was this translation helpful? Give feedback.
-
this looks more like an issue with the MPI library that is linked to your executable. running the configure script might have taken a wrong MPI library path. check what modules are loaded on your compute node and try to link against the correct ones with MPI_INC and MPI_LIB when running the configure script. |
Beta Was this translation helpful? Give feedback.
-
Beta Was this translation helpful? Give feedback.
-
not too bad - there are two differences in the setup:
to use the same DT in the new version you can just add a line in the Par_file with
having a different DT can affect the stability regime and accuracy of the solution. just to be safe, it's good to compare the solutions with the same DT value set.
this is something I missed in the old source setup, where the start time is only set to zero for force source solutions. the different start time will affect the length of the seismograms, although the source time function values are set to zero for times < 0. I'm not sure how much this would affect the filtering of the traces. anyway, I added this change to the updated devel version now as well (SPECFEM/specfem3d_globe#852). use |
Beta Was this translation helpful? Give feedback.
-
thanks Dorian, it looks like we're doing small steps towards a perfect fit... so based on your figures above, it seems that the main differences are still from attenuation. there must be a small discrepancy probably of the relaxation times between new and old version. given the traces only start differing at larger epicentral distances, these values must be only slightly different which might be due to difference floating point precisions. I'll try to find out more about that. in the meantime, could you attach the |
Beta Was this translation helpful? Give feedback.
-
Hi all,
I'm Dorian Soergel, I'm a postdoc with Barbara Romanowicz (@barbara-brz) at UC Berkeley. I have been trying to integrate SemUCB with its attached parametrisation and crust into SpecFEM3D_globe. For this there are four parts:
I have been working on this on this branch: https://github.com/dosoe/specfem3d_globe/tree/devel and have successfully implemented the source and the 1D model. I'm now working on the crust and the 3D perturbations.
As a first step, I would like to separate the crust and the 3D perturbations to be able to integrate them separately. Currently, I use the following in src/shared/get_model_parameters.f90 to integrate 1D model, 3D perturbations and custom crust:
For a test case that only uses the 1D model (no custom crust or 3D perturbations), I use the following settings in the same file:
If I want only to enable the crust, is it as simple as putting
CRUSTAL = .true.
andONE_CRUST = .true.
?As for the 3D part, I'm a little confused as for what
CASE_3D
andMODEL_3D_MANTLE_PERTUBATIONS
mean, are they redundant? or is it that the former means 'there is a 3D model' when the latter just means 'behave as if you had a 3D model, regardless whether it is true or not' ?Regards,
Dorian
Beta Was this translation helpful? Give feedback.
All reactions