Hackathon

[brendanjmeade]

@mallickrishg

Conceptual

• Comparisons between CS and quadratic elements
  • Wavy fault
  • Wavy fault with harmonic slip
  • Triple junctions
• Decide on a primary triple junction approach
  • Are there ways to further reduce the number of needed equations?
• Write out algebra for effective body force (Bannerjee/some textbook???)
• Where are stresses/tractions discontinuous?
• Is energy calculation in a volume bounded and meaningful? Does it converge with mesh refinement and an increasing number of observation coordinates?
• Is there and FEM comparison we can do for a wavy fault model?

Direct efforts

• RSF for quasi-dynamic
  • T. Ben Thompson direct DAE solve?
  • ODE formulation?
  • sinh regularization?
• Arbitrary traction free-surface solution
• Bookkeeping problem (see bem2d.jl)
• Dealing with spatially variable elastic moduli. Two layers first.
• HRF
  • I’ve got a semi-digitized version of the Hubbard geometry…needs more work
  • We can tackle the RSF problem on a non-planar fault. We could use the Hubbard geometry - I have started collaborating with Sharadha Sathiakumar recently and she has worked extensively with Judith on ramp-flat structures relevant to Nepal in 2-d and 3-d. I’m guessing she has text files of discretized MHT.
  • Is the strain energy calculation meaningful?
  • Has Ghorka coseismic to afterslip been done (https://www.nature.com/articles/s43247-021-00153-3)
• Integration of Kelvin point source over a right triangle

Papers possibilites

• Write abstract/outline/list figures for theory paper
• Write abstract/outline/list figures for HRF energy paper
• Write abstract/outline/list figures for Ghorka coseismic/afterslip energy paper
• Other targets:
  • CFS on non-planar faults
  • LA cross section (PHT/EPT) as imbricated example
  • Stress state/RSF evolution in rough/quasi-fractal fault systems

[mallickrishg]

I just emailed Sharadha requesting her 2-d geometry of the MHT from this paper

Looks like we have much to discuss and work on for the hackathon. I think we need to have some high level discussions about how we structure the theory paper and MHT paper, and most importantly what figures should go in. I think we will also need to sort out our ideas about triple junctions and non-planar faults.

Apart from that, I think the RSF problem on planar/non-planar faults and the heterogeneous elastic materials problems could be ideal targets for our limited time window.

Regarding Gorkha coafterslip: I need to show you some ideas about how I have been pursuing about coseismic + postseismic deformation problems in a rheology-agnostic manner, which is philosophically similar but mathematically different from your recent preprints on geometric moment conservation over the eq cycle. I think both of our approaches might reveal (at least i hope it will) something interesting about the spatial and temporal behaviour of Gorkha coseismic + afterslip.

[brendanjmeade]

I just emailed Sharadha requesting her 2-d geometry of the MHT from this paper

Great!

I think we will also need to sort out our ideas about triple junctions and non-planar faults.

This may be the most critical thing to work out together.

Apart from that, I think the RSF problem on planar/non-planar faults and the heterogeneous elastic materials problems could be ideal targets for our limited time window.

Decisions to be made and bookkeeping to get organized.

Regarding Gorkha coafterslip: I need to show you some ideas about how I have been pursuing about coseismic + postseismic deformation problems in a rheology-agnostic manner, which is philosophically similar but mathematically different from your recent preprints on geometric moment conservation over the eq cycle. I think both of our approaches might reveal (at least i hope it will) something interesting about the spatial and temporal behavior of Gorkha coseismic + afterslip.

This is one of the most interesting/intriguing clues I've read in a long time. I think there are real alternatives to what's been traditionally done. Looking forward to learning more.

[mallickrishg]

Wavy fault low & high-res done - https://github.com/brendanjmeade/bemcs/blob/main/notebooks/example_wavy_fault.ipynb

When you have time, could you push your new get_displacement_stress_kernel_constant() function? I will try to make the figures tonight, if i have time.

[brendanjmeade]

Nice with the wavy fault! Lets to the constant slip version of that first thing tomorrow and get some nice visualization out of it. The new constant slip function is now pushed.

[brendanjmeade]

Figures for theory/concept paper:

1. Rishav's slip basic flat fault figure
2. Constructing slip with three quadratic functions
3. Construction of smooth slip with quadratic elements
4. Quadratic vs. constant slip discretization
5. Quadratic vs. constant slip stresses at centroids and residual
6. Flat fault with smooth (Guassian?) slip distribution
   • CFS/Strain energy
   • Constant slip volume figure
   • Smooth slip volume figure
7. Convergence comparison
8. Non-planar fault smooth slip concept (?)
9. Wavy fault smooth vs. constant
   • Strain energy
   • CFS
10. Triple junctions (constant vs. smooth)

[mallickrishg]

Adding a note here to remind myself tomorrow regarding the final figure of this notebook - https://github.com/brendanjmeade/bemcs/blob/main/notebooks/example_constant_approximation_to_quadratic.ipynb

The values of $\sigma_{yy}$ are small but non-zero for the quadratic and constant slip elements because it is evaluated at a non-zero distance away from the fault. If we evaluated then exactly on the fault, both methods should give 0.

IMO, we should plot $\sigma_{xx}$ just as you did in your bem2d.jl paper to show that 1 component of the stress tensor is completely incorrect in the constant slip element calculation.

[mallickrishg]

@brendanjmeade I am not sure if this is what you foresaw coming out of the strain energy calculations - but voila!

[brendanjmeade]

This is a great figure!

I see how you've done the on-fault strain energy calculation. Looking at this, I can see that it can be positive or negative. However, the simple strain energy in the volume calculation, $U\sim \epsilon^2$, is always positive (at least by convention, we could imagine some slip event taking away strain energy). In the case you've calculated the sign changes are because sometimes the shear stress on the fault elements is positive and sometimes negative?

[mallickrishg]

In the case you've calculated the sign changes are because sometimes the shear stress on the fault elements is positive and sometimes negative?

Yes, exactly! I think the shear stress is mostly negative on the fault. That is expected because slip $sx$ results in negative $\sigma{xy}$ on-fault. I will check this using the strain energy in the volume calculation either tonight or tomorrow and get back to you.

You can see the on-fault $\sigma_{xy}$ here: https://github.com/brendanjmeade/bemcs/blob/main/notebooks/figure_convergence_strainenergy_comparison.ipynb

[brendanjmeade]

I get that, but what I don't get is how this is consistent with the general definition of strain energy with $U \sim \epsilon^2$ which is always positive. Something is inconsistent somewhere...or I'm misunderstanding something.

[mallickrishg]

However, the simple strain energy in the volume calculation, $U\sim\epsilon^2$, is always positive (at least by convention, we could imagine some slip event taking away strain energy)

I think this is true for purely axial strains, for 1-d beams and such. But does the strain energy density function also have to be positive? I am not sure.

Writing down stress and strain in plane strain, $\sigma{xy} = \frac{E}{(1+\nu)} \epsilon{xy}$, here we see that $\sigma{xy} \epsilon{xy} \geq 0$. So that is not interesting.

However, when we examine $\sigma{xx},\sigma{yy}$, there is no need for $\sigma$ and $\epsilon$ to have the same sign. $\sigma{xx} = \frac{E}{(1+\nu)(1-2\nu)}\left((1-\nu)\epsilon{xx} + \nu\epsilon{yy}\right)$ If $\epsilon{xx} <0$ and $\epsilon{yy}>0$, it is possible that $\sigma{xx}>0$ if $|\nu\epsilon{yy}| > |(1-\nu)\epsilon{xx}|$, thus the strain energy contribution would then be negative!

[mallickrishg]

Actually doing this more rigorously has me confused too! The strain energy density contributions from $\epsilon{xx},\epsilon{yy}$ terms are:

$W \propto (1-\nu)\left(\epsilon{xx}^2 + \epsilon{yy}^2\right) + 2\nu\epsilon{xx}\epsilon{yy}$ The first term is positive while the second term can be negative. If we assume $\nu = 0.25$, then we get $W \propto 0.75\left(\epsilon{xx}^2 + \epsilon{yy}^2\right) + 0.5\epsilon{xx}\epsilon{yy}$, and there is the identity that for any real numbers $a,b$, $a^2+b^2 \geq 2ab$ which implies that $W\geq 0$!

Such a roundabout way of saying you were right, and now I don't know what the $\int\tau_x s_x dx$ calculation's sign means!

[brendanjmeade]

This is certainly interesting. I think the $\tau_x sx$ is something, but I'm not sure it's elastic strain energy in the Malvern or Reddy continuum mechanics sense. Can you do the $\sigma{ij} \epsilon_{ij}$ calculation on the fault?

[brendanjmeade]

Back to basics:

$\sigma{ij} = \lambda \delta{ij} \epsilon{kk} + 2\mu\epsilon{ij}$

So strain energy is: $U = (\lambda \delta{ij} \epsilon{kk} + 2 \mu \epsilon{ij}) \epsilon{ij}$

So we're not squaring strains in the general case. Maybe it can go negative?

[mallickrishg]

Using that formula for plane strain, $U = (\lambda +2\mu)(\epsilon{xx}^2 + \epsilon{yy}^2) + 2\lambda\epsilon{xx}\epsilon{yy} + 2\mu\epsilon{xy}^2)$, the only non-squared term is $2\lambda\epsilon{xx}\epsilon_{yy}$.

And as i noted in a post above, $(\epsilon{xx}^2 + \epsilon{yy}^2) \geq 2\epsilon{xx}\epsilon{yy}$, so that proves $U\geq 0$.

[mallickrishg]

I was just looking through some older articles about work partitioning and strain energy calculations for earthquakes, and realized that the sign convention for strain energy is $U = (\tau_i - \tau_f) \frac{s}{2}$ where $\tau_i$ is the initial stress and $\tau_f$ is stress at the end of the slip event. So $\tau_i - \tauf$ is the opposite sign of stress change that we calculate! So maybe we just need to calculate $-\int{\partial\Omega}\tau_x s_x + \tau_y s_y dx$?

I just played around with this calculation, and used random numbers for slip and the result is always positive!

[mallickrishg]

This is certainly interesting. I think the τxsx is something, but I'm not sure it's elastic strain energy in the Malvern or Reddy continuum mechanics sense. Can you do the σijϵij calculation on the fault?

Just checked. The results appear to be a rescaled version of the $\int \tau s dx$ value.

[brendanjmeade]

So does that mean the sign still flip positive to negative?

[mallickrishg]

$\int\tau.s \partial\Omega$ is always negative so $-\int\tau.s \partial\Omega$ is always positive. This is applicable to continuous slip elements. With constant slip, it really depends on where we sample - the sign can flip based on that. Basically the strain energy calculation for constant slip elements at points apart from the patch centers is not meaningful.

[mallickrishg]

Also this strain energy calculation exercise has convinced me that we should use $\Delta\tau$ and $\Delta\sigma$ instead of $\tau,\sigma$, because we are only dealing with changes in stress. Specifically it is $\sigma{final}-\sigma{initial}$.

[mallickrishg]

To clear any doubts, I evaluated stress carefully at points that were away from the singularities for constant slip elements and then calculated strain energy. It turns out with this sort of sampling, the strain energy calculation from constant slip tends to what we get from quadratic elements. There is still large error in the constant slip element calculations, but at least it doesn't oscillate.

Below is the stress change evaluated at those points.

Strain energy comparison

[brendanjmeade]

points that were away from the singularities

Where are they? Still on the fault at centroids or somewhere in the volume around the fault?

Also, here's a version of figure you made with fancy line styles. Not sure I like. Just playing around.

[mallickrishg]

The points I evaluate are exactly ON the fault, just not on the end points of each fault element. I chose the x values of observations shifted away from els.x1, els.x2 end points as follows:

for j in range(n_els): x_obs[j * npts : (j + 1) * npts] = np.linspace(0.1*els.half_lengths[j] + els.x1[j], els.x2[j] - 0.1*els.half_lengths[j],npts)

About the figure, I like that you switched the x axis to log. My request is that we use the same colors for constant slip and continuous slip elements for all our figures. In our other figures we use nice_blue and nice_orange, and I vote that we stick with those. Of course if you found better colors, I will change my mind.

[brendanjmeade]

Fixed the color and pushed.

[mallickrishg]

Going over stuff for the theory paper

1. Rishav's slip basic flat fault figure .. done
2. Constructing slip with three quadratic functions .. done
3. Construction of smooth slip with quadratic elements .... done (need to add second subplot to show stress is smooth too)
4. Quadratic vs. constant slip discretization .... done
5. Quadratic vs. constant slip stresses at centroids and residual .... done
6. Flat fault with smooth (Guassian?) slip distribution .... done
7. CFS ~/Strain energy~ .... done
8. Constant slip volume figure .... done
9. Smooth slip volume figure .... done
10. Convergence comparison..... done for strain energy
11. Non-planar fault smooth slip concept (?) ....... need to discuss and not done yet
12. Wavy fault smooth vs. constant ..... done
13. ~Strain energy~ done $I_2$ instead
14. CFS part of a figure already
15. Triple junctions (constant vs. smooth) ....................... not done yet, and perhaps we don't need to do constant slip? We have bashed constant slip a lot.

We got a lot of things done last week! Kudos!

I have started adding some text to the paper, mainly the method. I will continue adding what I can as I find time. The last main figure is the triple junction stuff. I think there is an additional cartoon needed for the triple junction to show how we are dealing with the node shared by 3 elements. Thoughts?

[brendanjmeade]

11. Non-planar fault smooth slip concept (?) ....... need to discuss and not done yet

I can't even remember what this meant. Do you?

15. Triple junctions (constant vs. smooth) ....................... not done yet, and perhaps we don't need to do constant slip?

I think we do have the figures for both already in: https://github.com/brendanjmeade/bemcs/blob/main/notebooks/example_triple_junctions_120.ipynb I can go either way with showing results for the constant slip case. Happy to leave it out.

I think there is an additional cartoon needed for the triple junction to show how we are dealing with the node shared by 3 elements.

Maybe quiver plot with slip evaluated on the elements. Something like 10 evaluation points on each of the three segments. It might not look great but it's pretty clear.

Also, I've added a new notebook that creates sort of honeycomb geometry and has a set of slips (in $x$, $y$ coordinates). The purpose of this geometry is to show off the use of triple junctions for an applied problem. We mentioned this problem in section 3 of the Overleaf doc. Do you think you can run a forward model with this and do the volume evaluation?

[mallickrishg]

21. Non-planar fault smooth slip concept (?) ....... need to discuss and not done yet

I can't even remember what this meant. Do you?

I didn't yesterday, but as I started writing things in overleaf, this point made sense. I think we need a figure like this, for non-planar faults. The problem is I don't know how to visualize this. I am imagining a figure where the non-planar fault geometry is on the x-y plane, and somehow on top of that we show slip?

Also, I've added a new notebook that creates sort of honeycomb geometry and has a set of slips (in coordinates). The purpose of this geometry is to show off the use of triple junctions for an applied problem. We mentioned this problem in section 3 of the Overleaf doc. Do you think you can run a forward model with this and do the volume evaluation?

Awesome and amazing! Thanks for doing this - it's exactly as you drew it in your office! I am looking forward to this task :) - this will be my target for tomorrow morning.

[brendanjmeade]

I didn't yesterday, but as I started writing things in overleaf, this point made sense. I think we need a figure like this, for non-planar faults. The problem is I don't know how to visualize this. I am imagining a figure where the non-planar fault geometry is on the x-y plane, and somehow on top of that we show slip?

I see. What about doing this as quiver plots (along a line or curved fault) as I suggested for the triple junction? Could be the best way to do this efficiently. What do you think?

Awesome and amazing! Thanks for doing this - it's exactly as you drew it in your office! I am looking forward to this task :) - this will be my target for tomorrow morning.

Sounds awesome. We're going to have all the figures done soon!

[mallickrishg]

What about doing this as quiver plots (along a line or curved fault) as I suggested for the triple junction? Could be the best way to do this efficiently. What do you think?

I will try it out. It sounds pretty good to me.

[mallickrishg]

Good news is that the automatic node labelling system works correctly, and I think we have a semi-automatic method of constructing and solving BEM problems, at least for static linear problems: https://github.com/brendanjmeade/bemcs/blob/main/notebooks/example_honeycomb.ipynb

I am still working out what the best way of computing and plotting fault slip values for the honeycomb is. Our functions give us shear and tensile slip for a horizontal fault, so I guess I need to rotate each of these from (s,t) to (x,y) coordinates? So I discretize each fault segment with N points, and calculate the (x,y) slip for each point, and then repeat for next segment........? That makes sense right?

[brendanjmeade]

Good news is that the automatic node labelling system works correctly, and I think we have a semi-automatic method of constructing and solving BEM problems, at least for static linear problems:

That's fantastic news. Whether automated bookkeeping works is always unknown!

I am still working out what the best way of computing and plotting fault slip values for the honeycomb is. Our functions give us shear and tensile slip for a horizontal fault, so I guess I need to rotate each of these from (s,t) to (x,y) coordinates?

I think the thing to do is project the (x, y) slip coordinates to (s, t) and then we should be able to use those values as arguments for to the quadratic element functions which work for segments of any orientation.

So I discretize each fault segment with N points, and calculate the (x,y) slip for each point, and then repeat for next segment........? That makes sense right?

Not sure we need to discretize. First I'd run this at the resolution we have.

[mallickrishg]

I think the thing to do is project the (x, y) slip coordinates to (s, t) and then we should be able to use those values as arguments for to the quadratic element functions which work for segments of any orientation.

So I discretize each fault segment with N points, and calculate the (x,y) slip for each point, and then repeat for next segment........? That makes sense right?

Not sure we need to discretize. First I'd run this at the resolution we have.

I am now plotting slip at the center and end points of each element. Was that what you suggested?

Also, I can now see how slip is divided among the 3 fault elements at a triple junction. That's pretty neat!

[brendanjmeade]

The TJ slip partitioning is so nice to see! Help me to understand what slip is here. In the first honeycomb notebook slip was given in x, y coordinates with no y component. What are these? Is this some representation of shear and normal slip?

[mallickrishg]

I changed the slip to have non zero x and y components while being symmetric. The stress components look really cool. I've pushed these updates to example_honeycomb.ipynb

[brendanjmeade]

Not sure about this. See notes in the attached image. Lot's to learn!

[mallickrishg]

I changed the x_slips,y_slips from the values you had initially set to some new values that had variations in both x and y directions just so that I could visualize things more clearly. That's why everything looks so different from your expectation that there should be a gradient in x_slips while y_slips=0. I hope I understood your concern correctly. Anyway, later tonight I will revert back to the values you had imposed and update the plot and post here.

[mallickrishg]

Here are the new figures with the original slip.

[brendanjmeade]

@mallickrishg This is looking great. What file is the implemented in? I've pulled the latest version example_honeycomb.ipynb but it seems to still have the older slip rates in it.

[mallickrishg]

@brendanjmeade - oops, I had forgotten to push changes. Done now.

[mallickrishg]

@brendanjmeade I also added some functionality for plotting slip on a non-planar fault in map/cross-sectional view for an arbitrary number of points per element. So far we could only do it for the end points and center. Do you think we should add this to bemcs monolith as a plotting function?

[brendanjmeade]

Jaw. Dropped.

Simply awesome. I think that's the best figure we have for the paper. Definitely add it to the bemcs.py!

[mallickrishg]

Perhaps it would be more concise if we combined plot with both slip and tractions? I still need your help with choosing colors and some magic like you did with the continuous slip for 2 elements figure.

The nonplanar fault notebook is here - https://github.com/brendanjmeade/bemcs/blob/main/examples/example_non_planar_fault_plotting.ipynb

[brendanjmeade]

I think the quiver one is enough. It's simply great and requires no interpretation. I think the colors are great too! I've adjusted the ticks a bit but other than that I think it's good to go.

[mallickrishg]

@brendanjmeade - question about figures in the paper.

I had initially thought to add the quiver plot above as a figure, but then I realised you already made a really nice example of constant slip on a wavy fault. That wavy fault example has really nice figures of displacements, $I_2$ and $\Delta$CFS. I would like to make a quiver plot for that geometry and use that in the paper instead of the 3-segment one from above. Could you provide the name of the notebook that makes this figure: https://github.com/brendanjmeade/bemcs/blob/main/figures/six_panel_wavy_constant_slip.pdf

[brendanjmeade]

Sounds good. The file you're looking for is: example_wavy_fault_volume_figures.ipynb. Currently, it's set to not be wavy (Y_AMPLITUDE = 0.0), but if you just increase the magnitude of the global variable Y_AMPLITUDE.