Julia tutorial - GeoEnergyLab-EPFL/BigWham GitHub Wiki

Griffith Crack (2D)

# Julia code blocks use ##
## %%
using Pkg
Pkg.add("Plots")
Pkg.add("PyCall")
Pkg.add("LinearAlgebra")
Pkg.add("SparseArrays")
Pkg.add("LaTeXStrings")

## %%
using PyCall
using LinearAlgebra
using SparseArrays
using Plots
using LaTeXStrings

## %% # Add Python path for the bigwham4py module
home = ENV["HOME"]
bigwham_module_path = home * "/geolab/BigWham/build/interfaces/julia"
include(bigwham_module_path * "/bigwham4jl.jl")
import .BigWham:BEMatrix, build_hierarchical_matrix


## %% Material properties
G = 1.0;
nu = 0.25;
E = (2 * G) * (1 + nu);


## %% Mesh
a = 10;
nelts = 10000;
coor1D = LinRange(-a, a, nelts + 1);
coor = hcat(coor1D, zeros(nelts + 1));
conn = [i + j for i in 0:nelts-1, j in 0:1];

## Hmatrix
# H-matrix parameter
max_leaf_size = 100
eta = 3.0
eps_aca = 1.e-4

kernel = "2DP0-H"
elas_prop = [E, nu]
num_threads = 8
hmat = BEMatrix(coor, conn, kernel, elas_prop);
build_hierarchical_matrix(hmat, max_leaf_size, eta, eps_aca)
col_pts = hmat.col_pts

## Test hmat matvec timing
dd = rand(h.shape[2]);
num_iter = 10000
# Timing the loop
total_time = @elapsed begin
    for i in 1:num_iter
        hmat * dd
    end
end

println("Total time for $num_iter iterations: $total_time seconds")

Scaling of Griffith crack with OMP threads

## %%
using Pkg
Pkg.add("Plots")
Pkg.add("PyCall")
Pkg.add("LinearAlgebra")
Pkg.add("SparseArrays")
Pkg.add("LaTeXStrings")

## %%
using PyCall
using LinearAlgebra
using SparseArrays
using Plots
using LaTeXStrings


home = ENV["HOME"]
bigwham_module_path = home * "/geolab/BigWham/build/interfaces/julia"
include(bigwham_module_path * "/bigwham4jl.jl")
import .BigWham: BEMatrix, build_hierarchical_matrix

## %% Material properties
G = 1.0;
nu = 0.25;
E = (2 * G) * (1 + nu);


## %% Mesh
a = 10;
nelts = 10000;
coor1D = LinRange(-a, a, nelts + 1);
coor = hcat(coor1D, zeros(nelts + 1));
conn = [i + j for i in 0:nelts-1, j in 0:1];

## Hmatrix
# H-matrix parameter
max_leaf_size = 100
eta = 3.0
eps_aca = 1.e-4

kernel = "2DP0-H"
elas_prop = [E, nu]
# for num_threads in [4, 8, 16, 32, 64, 96, 128]
for num_threads in [4, 8, 16, 32, 64, 96, 128]
# for num_threads in [32]
    # println("Number of threads: $num_threads")
    h = BEMatrix(coor, conn, kernel, elas_prop, num_threads)
    build_hierarchical_matrix(h, max_leaf_size, eta, eps_aca)

    ## Test hmat matvec timing
    dd = rand(size(h)[2]);
    num_iter = 10000
    # Timing the loop
    t = similar(dd)
    mul!(t, h, dd)
    total_time = @elapsed begin
        for i in 1:num_iter
            t = similar(dd)
            mul!(t, h, dd)
        end
    end

    println("Total time for $num_iter iterations: $total_time [s] with $num_threads threads")
end

Griffith Crack (2D): Inverting for slip_rate given shear stressing rate

using PyCall
using LinearAlgebra
using SparseArrays
using Plots
using LaTeXStrings
using Krylov
using Printf

home = ENV["HOME"]
bigwham_module_path = home * "/geolab/BigWham/build/interfaces/julia"
include(bigwham_module_path * "/bigwham4jl.jl")
import .BigWham: BEMatrix, build_hierarchical_matrix


## %% Material properties
G = 30e9;
nu = 0.25;
E = (2 * G) * (1 + nu);


## %% Mesh
a = 10;
nelts = 1000;
coor1D = LinRange(-a, a, nelts + 1);
coor = hcat(coor1D, zeros(nelts + 1));
conn = [i + j for i in 0:nelts-1, j in 0:1];
W = 2 * a

## Hmatrix
# H-matrix parameter
max_leaf_size = 100
eta = 3.0
eps_aca = 1.e-4

kernel = "2DP0-H"
elas_prop = [E, nu]
num_threads = 16
h = BEMatrix(coor, conn, kernel, elas_prop, num_threads)
build_hierarchical_matrix(h, max_leaf_size, eta, eps_aca)
dd = rand(size(h)[2]);
t = similar(dd)
mul!(t, h, dd)

# shear stressing rate
t = ones(size(h)[2]) * 1e-3
t[2:2:end] .= 0

# Without preconditioning
x, stats = bicgstab(h, t, history=true)
r = t - h * x
@printf("[Without preconditioning] Residual norm: %8.1e\n", norm(r))
@printf("[Without preconditioning] Number of iterations: %3d\n", length(stats.residuals) - 1)

col_pts = h.col_pts
vplate_scale = t[1] * (W / G)
plot(col_pts[1:2:end], x[1:2:end] ./ vplate_scale, label=false, xlabel="x [m]", 
            ylabel="v / v_plate", 
            title="slip rate along fault")