Skip to content

Public Interface

Index

Interfaces

# JuEQ.DieterichStateLaw — Type.

$\frac{\mathrm{d}θ}{\mathrm{d}t} = 1 - \frac{V θ}{L}$

# JuEQ.PrzStateLaw — Type.

$\frac{\mathrm{d}θ}{\mathrm{d}t} = 1 - (\frac{V θ}{2L})^2$

# JuEQ.RuinaStateLaw — Type.

$\frac{\mathrm{d}θ}{\mathrm{d}t} = -\frac{V θ}{L} * \log{\frac{V θ}{L}}$

# JuEQ.DECallbackSaveToFile — Method.

1
DECallbackSaveToFile(iot::IOStream, iou::IOStream)

Construct a functional callback to write ODESolution (t & u) into file. The reason to separate t and u is for more easily reshape u w.r.t grids specification. It right now falls on users' memory on what the type of solution is for accurately retrieving results.

Arguments

  • iot::IOStream: stream pointing to solution of time
  • iou::IOStream: stream pointing to solution of domain

Note It is strongly not recommended to use "skipping" scheme (by defining thrd and dts(a) for each case) when solution is too oscillated.

# JuEQ.EarthquakeCycleProblem — Method.

1
EarthquakeCycleProblem(p::PlaneMaterialProperties, u0, tspan; se=DieterichStateLaw(), fform=CForm())

Return an ODEProblem that encapsulate all the parameters and functions required for simulation. For the entailing usage, please refer DifferentialEquations.jl

Arguments

  • gd::BoundaryElementGrid: grids for fault domain.
  • p::PlaneMaterialProperties: material profile.
  • u0::AbstractArray: initial condition, should be organized such that the first of last dim is velocity while the 2nd of last dim is state.
  • tspan::NTuple: time interval to be simulated.
  • se::StateEvolutionLaw: state evolution law to be applied.
  • fform::FrictionLawForm: forms of frictional law to be applied.

# JuEQ.dc3d_okada — Method.

Calculate displacements and gradient of displacements due to a dislocation in an elastic isotropic halfspace. See dc3d for details.

test/test_okada.dat is obtained using DC3dfortran

An example wrapper for DC3D in julia as below:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
function dc3d_fortran(x::T, y::T, z::T, α::T, dep::T, dip::T, al1::T, al2::T, aw1::T, aw2::T,
    disl1::T, disl2::T, disl3::T) where {T <: AbstractFloat}

    # initial return values
    # `RefValue{T}` may be also viable other than `Array{T, 1}`
    ux = Array{Float64}(1)
    uy = Array{Float64}(1)
    uz = Array{Float64}(1)
    uxx = Array{Float64}(1)
    uyx = Array{Float64}(1)
    uzx = Array{Float64}(1)
    uxy = Array{Float64}(1)
    uyy = Array{Float64}(1)
    uzy = Array{Float64}(1)
    uxz = Array{Float64}(1)
    uyz = Array{Float64}(1)
    uzz = Array{Float64}(1)
    iret = Array{Int64}(1)

    # call okada's code which is renamed as "__dc3d__" (see binding rename shown below)
    # input args tuple must be syntactically written instead of a variable assigned
    # macros could be used to simplify this in the future
    ccall((:__dc3d__, "dc3d.so"), Void,
        (
            Ref{Float64},
            Ref{Float64}, Ref{Float64}, Ref{Float64}, Ref{Float64},
            Ref{Float64}, Ref{Float64}, Ref{Float64}, Ref{Float64},
            Ref{Float64}, Ref{Float64}, Ref{Float64}, Ref{Float64},
            Ptr{Array{Float64,1}}, Ptr{Array{Float64,1}}, Ptr{Array{Float64,1}},
            Ptr{Array{Float64,1}}, Ptr{Array{Float64,1}}, Ptr{Array{Float64,1}},
            Ptr{Array{Float64,1}}, Ptr{Array{Float64,1}}, Ptr{Array{Float64,1}},
            Ptr{Array{Float64,1}}, Ptr{Array{Float64,1}}, Ptr{Array{Float64,1}},
            Ref{Int64},
        ),
        α, x, y, z, dep, dip, al1, al2, aw1, aw2, disl1, disl2, disl3,
        ux, uy, uz, uxx, uyx, uzx, uxy, uyy, uzy, uxz, uyz, uzz,
        iret,
    )

    # results valid iff iret[1] == 0
    return (
        iret[1],
        ux[1], uy[1], uz[1],
        uxx[1], uyx[1], uzx[1],
        uxy[1], uyy[1], uzy[1],
        uxz[1], uyz[1], uzz[1]
    )
end

The corresponding fortran module is:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
MODULE okada
  USE, INTRINSIC :: iso_c_binding
  IMPLICIT NONE
CONTAINS

  SUBROUTINE dc3d_wrapper(&
       & alpha, &
       & x, y, z, &
       & depth, dip, &
       & al1, al2, &
       & aw1, aw2, &
       & disl1, disl2, disl3, &
       & ux, uy, uz, &
       & uxx, uyx, uzx, &
       & uxy, uyy, uzy, &
       & uxz, uyz, uzz, &
       & iret) BIND(C, NAME='__dc3d__')

    REAL*8 :: &
         & alpha, &
         & x, y, z, &
         & depth, dip, &
         & al1, al2, &
         & aw1, aw2, &
         & disl1, disl2, disl3, &
         & ux, uy, uz, &
         & uxx, uyx, uzx, &
         & uxy, uyy, uzy, &
         & uxz, uyz, uzz

    INTEGER*8 :: iret

    CALL dc3d(&
         & alpha, &
         & x, y, z, &
         & depth, dip, &
         & al1, al2, &
         & aw1, aw2, &
         & disl1, disl2, disl3, &
         & ux, uy, uz, &
         & uxx, uyx, uzx, &
         & uxy, uyy, uzy, &
         & uxz, uyz, uzz, &
         & iret)

  END SUBROUTINE dc3d_wrapper

END MODULE okada

A sample of makefile is as below:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
# Build Okada's code for calculating deformation due to a fault model
#
CC = gfortran
CFLAGS = -fPIC -w -O3
LDFLAGS = -shared

SRCS = dc3d.f okada.f90
OBJS = $(SRCS:.c=.o)

TARGET = dc3d.so

$(TARGET): $(OBJS)
    $(CC) $(CFLAGS) $(LDFLAGS) -o $(TARGET) $(OBJS)

# JuEQ.discretize — Method.

1
discretize(fa::PlaneFaultDomain{ftype, 1}, Δξ::T; ax_ratio=12.5)

Generate the grid for given 1D fault domain. The grids will be forced to start at (x=0, y=0, z=0).

Arguments

  • Δξ: grid space along-downdip
  • ax_ratio::Number: ration of along-strike length agsinst along-downdip length for mimicing an extended 2d (x & ξ) fault represented by 1d (ξ) domain. Default ax_ratio=12.5 is more than enough for producing consistent results.

# JuEQ.discretize — Method.

1
discretize(fa::PlaneFaultDomain{ftype, 2}, Δx, Δξ; buffer=:auto) where {ftype <: PlaneFault}

Generate the grid for given 2D fault domain. The grids will be forced to start at (z=0) and spread symmetrically along x-axis w.r.t y-z plane. By such setting, we would be able to utilize the symmetry properties of stiffness tensor for performance speed up.

Arguments

  • Δx, Δξ: grid space along-strike and along-downdip respectively
  • `buffer_ratio::Number: ration of buffer size against along-strike length for introducing zero-dislocation area at along-strike edges of defined fault domain.

# JuEQ.fault — Method.

1
fault(ftype::Type{<:PlaneFault}, dip, span)

Generate a fault given the fault type, dip angle and its spatial span.

Arguments

  • ftype::Type{<:PlaneFault}: type of plane fault
  • dip: dip angle in degree
  • span: spatial span of fault size

# JuEQ.friction — Method.

1
friction(::FrictionLawForm, v::T, θ::T, L::T, a::T, b::T, f0::T, v0::T) where {T<:Number}

Calculate friction given by the form of fomula as well as other necessary parameters.

  • Conventional Form:

  • Regularized Form:

# JuEQ.max_velocity — Method.

1
max_velocity(t::AbstractVector, u::AbstractArray, getu::Function)

Return max velocity across the fault at each time step. A number of convenient interfaces for common output are implemented.

Arguments

  • t::AbstractVector: vector of time steps
  • u::AbstractArray: array of solution
  • getu::Function: method for retrieving velocity section at each time step

# JuEQ.moment_magnitude — Method.

Calculate moment magnitude.

# JuEQ.properties — Method.

1
properties(fa::PlaneFaultDomain, gd::BoundaryElementGrid{dim}; _kwargs...) where {dim}

Establishing a material-properties-profile given by the fault domain and grids. User must provide the necessary parameters in according to the grid size specified or just a scalar for broadcasting.

Arguments that are required:

  • a: contrib from velocity.
  • b: contrib from state.
  • L: critical distance.
  • σ: effective normal stress.
  • η: radiation damping. It is recommended to set as $μ / 2\mathrm{Vs}$ where $μ$ is shear modulus and $\mathrm{Vs}$ shear wave velocity.
  • vpl: plate rate.
  • f0: ref. frictional coeff.
  • v0: ref. velocity.

Arguments that need options

  • k: stiffness tensor.

    (1) Providing shear modulus denoted as μ and Lamé's first parameter denoted as λ (same as μ if missing), then calculate it based on grid and fault domain, choosing parallel scheme if nprocs() != 1. (2) an AbstractArray represent the pre-calculated stiffness tensor. No verification will be performed here.

# JuEQ.stiffness_tensor — Method.

1
stiffness_tensor(fa::PlaneFaultDomain, gd::BoundaryElementGrid, ep::HomogeneousElasticProperties)

Calculate the reduced stiffness tensor. For 2D fault, the final result will be dimensionally reduced to a 3D array due to the translational & reflective & perodic symmetry, such that the tensor contraction will be equivalent to convolution, hence we could use FFT for better performace.

Note

  • Faults are originated from surface and extends downwards, thus dep = 0