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.
-
thanks for the ellipticity comment, Barbara! If I understand correctly, both crustal and mantle models are thus defining velocities with reference to geocentric latitude/longitude positions. that helps already to set up something correspondingly. Dorian, one of the culprints I found in the new version is this line here in file instead of
you need to define them now as
the try correcting this in the new version, and it should bring your match between versions much closer to each other. most other differences I see are due to how the new version handles the geocentric/geographic conversions for ellipticity. I'll need to rewrite this a bit to allow different (crustal) models with different reference frames. for example, the CRUST1.0/2.0 models use geographic lat/lon positions, while yours is using geocentric lat/lon positions - something to rethink in the new version. another even more subtle point would be to adjust the layer thicknesses when reading in and assigning crustal models that are referencing a geographic (elliptical) framework like the CRUST models. we assign the model velocities while the mesh is still spherical, and then stretch it out afterwards to account for the elliptical shape. this will affect the layer thicknesses slightly in the final mesh and we could in principle correct this such that the thicknesses in the final, elliptical mesh are exactly the ones as defined - this might have a minor effect and has been ignored till now - let me check. or maybe, I'm just wrong about the geographic reference frame of the CRUST1.0/2.0 models and we're all set already. I'll need to go through those model descriptions, again... :) finally, regarding the speed of reading in the model, I'll see what I can do to speed it up. most models read much faster, so I'm not used to wait for so long until the mesher has finished. hopefully, there are some magic tricks like hash tables or precomputed factor arrays to make this fast, i'll check once the "initial" version has been added. |
Beta Was this translation helpful? Give feedback.
-
Hi Daniel and Dorian,
I think the problem is resolved and ***WE"RE DONE *** - no need to scratch
our heads about any implementation of the 3D crust etc..
1) I changed the line in model_berkeley.f90, but more importantly
2) I modified constants.h so that:
logical, parameter :: ATTENUATION_1D_WITH_3D_STORAGE = .true.
(I thought I had modified it earlier, but for some reason it did not
register - my bad! that was the main problem from day 1 and Dorian only
recently spotted it)..
and the fit is now perfect!
I attach a couple of examples at different distances here, but *Dorian
will need to port them to the github and check whether things also work for
a shallow event (mine is deep*):)
Daniel: I hope you find the same with the version that Dorian sent you and
we can officially declare that SEMUCB (and other models in Berkeley 's
format) are now implemented in the latest version of SpecFEM, at least for
CPU's.
Regards
Barbara
…On Sun, Oct 6, 2024 at 1:09 PM daniel peter ***@***.***> wrote:
thanks for the ellipticity comment, Barbara! If I understand correctly,
both crustal and mantle models are thus defining velocities with reference
to geocentric latitude/longitude positions. that helps already to set up
something corresponding.
Dorian, one of the culprints I found in the new version is this line here
in file model_berkeley.f90:
https://github.com/dosoe/specfem3d_globe/blob/41ac1a0e5bce77eca8b2752700d97a6c37381fc8/src/meshfem3D/model_berkeley.f90#L339
instead of
real(kind=4) :: dvsv,dvsh,dvpv,dvph
you need to define them now as
double precision :: dvsv,dvsh,dvpv,dvph
the real(kind=4) parameters have only been used for the S362ANI models in
the past. to separate and distinguish these "special" kind parameters, they
have been renamed to xdvsv,.. etc. in the new version. however, double
precision :: dvsv,.. parameters are now also used for other models (e.g.,
the full spherical harmonic model) in the calling routine in file
meshfem3D_models.F90. so, you can either change the definition in
model_berkeley.f90 to double precision as stated above, or leave as is
and change the calling routine to use xdvsv, xdvsh, xdvpv, xdvph for
those instead. I would prefer the former and keep the single precision
real(kind=4) parameters only for the S362ANI models.
try correcting this in the new version, and it should bring your match
between versions much closer to each other.
most other differences I see are due to how the new version handles the
geocentric/geographic conversions for ellipticity. I'll need to rewrite
this a bit to allow different (crustal) models with different reference
frames. for example, the CRUST1.0/2.0 models use geographic lat/lon
positions, while yours is using geocentric lat/lon positions - something to
rethink in the new version.
another even more subtle point would be to adjust the layer thicknesses
when reading in and assigning crustal models that are referencing a
geographic (elliptical) framework like the CRUST models. we assign the
model velocities while the mesh is still spherical, and then stretch it out
afterwards to account for the elliptical shape. this will affect the layer
thicknesses slightly in the final mesh and we could in principle correct
this such that the thicknesses in the final, elliptical mesh are exactly
the ones as defined - this might have a minor effect and has been ignored
till now - let me check. or maybe, I'm just wrong about the geographic
reference frame of the CRUST1.0/2.0 models and we're all set already. I'll
need to go through those model descriptions, again... :)
finally, regarding the speed of reading in the model, I'll see what I can
do to speed it up. most models read much faster, so I'm not used to wait
for so long until the mesher has finished. hopefully, there are some magic
tricks like hash tables or precomputed factor arrays to make this fast,
i'll check once the "initial" version has been added.
—
Reply to this email directly, view it on GitHub
<#1718 (comment)>,
or unsubscribe
<https://github.com/notifications/unsubscribe-auth/BJEVZDPGA5FM2744OIRFF6TZ2GKGRAVCNFSM6AAAAABKSRHCXGVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTAOBWGA3TMOA>
.
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