Merge branch 'main' of https://github.com/tduretz/SeismicQ

Project.toml

@@ -5,8 +5,12 @@ version = "0.1.0"
 
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 SymbolicNumericIntegration = "78aadeae-fbc0-11eb-17b6-c7ec0477ba9e"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

README.md

@@ -1,6 +1,6 @@
-# SeismicQ
-
 <h1> <img src="docs/src/assets/logo.png" alt="SeismicQ.jl" width="50"> SeismicQ </h1>
 
 [![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://tduretz.github.io/SeismicQ/dev/)
 [![Build Status](https://github.com/tduretz/SeismicQ/workflows/CI/badge.svg)](https://github.com/tduretz/SeismicQ/actions)
+
+<h1> <img src="docs/src/assets/visuel_ecocup2.png" alt="SeismicQ.jl" width="500">  </h1>

docs/src/index.md

@@ -8,4 +8,7 @@ Ricker
 ```
 ```@docs
 f_bulk
+Spec
+ComputeQgraph
+Getfreq
 ```

examples/FausseTrace.jl

@@ -0,0 +1,89 @@
+# Function that generates a seismic trace (Ricker function), attenuating along time, 
+# and plot a receiver gather (received waves at different geophone positions)
+
+using SeismicQ, Plots, SpecialFunctions, LinearAlgebra, Printf
+
+function main()
+
+    # Functtion that create a vector of time and 
+    function FausseTrace(x,Δt,Nt,t,Vp,Vs,αp,αs,𝑓₀,t₀)
+
+        tₓp = x/Vp
+        tₓs = x/Vs
+
+        # Storage
+        time = zeros(Nt)
+        acc = zeros(Nt)
+
+        # Time loop
+        for it=1:Nt
+
+            # Compute Ricker function
+            t += Δt
+            a = 0.5 * exp(-αp*x)*Ricker(t, tₓp+t₀, 𝑓₀)+ 0.5 * exp(-αs*x)*Ricker(t, tₓs+t₀, 𝑓₀)
+
+            # For visualisation purpose
+            time[it] = t
+            acc[it]  = a
+        end
+
+        return(time',acc')
+
+    end
+
+    # Geophone position [m]
+    listₓ = 0:100:5000;
+
+    # Time domain
+    Δt  = 1e-3
+    Nt  = 2000
+    # Central frequency of the source [Hz]
+    𝑓₀  = 10.
+    t₀ = 1.0/𝑓₀
+    t   = -t₀
+
+    # Velocities
+    Vp = 7000 # m/s
+    Vs = 4000 
+    αp = 2e-4
+    αs = 4e-4 # (π * f)/ (Q * V)
+
+    time_axis = t:(Δt*Nt-t)/(Nt-1):Δt*Nt ;
+    time_vec = zeros(size(listₓ,1),Nt);
+    acc_vec  = zeros(size(listₓ,1),Nt);
+    time_vec[1,:],acc_vec[1,:]= FausseTrace(listₓ[1,1],Δt,Nt,t,Vp,Vs,αp,αs,𝑓₀,t₀);
+    dist_axis = 0:1: size(listₓ,1)-1; # distances des geophones pour l instant 1,2,3,...
+
+    for i=1:size(listₓ,1)-1
+        time_vec[i+1,:],acc_vec[i+1,:] = FausseTrace(listₓ[i+1,1],Δt,Nt,t,Vp,Vs,αp,αs,𝑓₀,t₀)
+    end
+
+    seismic_matrix = hcat(listₓ,acc_vec);
+    @show size(seismic_matrix)
+    #=
+    # Visualisation si peu de positions
+    fig1   = plot(layout = (size(listₓ,1),1)) 
+    p1     = plot!(fig1[1],time_vec[1,:], acc_vec[1,:], xlabel="t", ylabel="a")
+
+    for i=1:size(listₓ,1)-1
+        p1 = plot!(fig1[i+1],time_vec[i+1,:], acc_vec[i+1,:], xlabel="t", ylabel="a")
+    end
+    display(fig1)
+    =#
+
+    # Visualisation of receiver gather
+    #p2 = heatmap(dist_axis,time_axis,acc_vec', 
+    p2 = heatmap(listₓ,time_axis,acc_vec', 
+        color=palette(:RdBu, 100, rev=true),
+        clim=(-0.5, 0.5),
+        cbar=true,
+        label=" ", 
+        yflip=true,
+        title = "Receiver gather", 
+        xlabel = "Geophone position [m]",
+        ylabel = "Time [s]")
+    display(plot(p2)) # equivalent du drawnow
+
+
+end
+main()

examples/GetQTest.jl

@@ -0,0 +1,86 @@
+using SeismicQ,Plots
+
+function main()
+
+#P and S velocities
+Vp=7000
+Vs=4000
+
+SeismicMatrix=GenMatrix(Vp,Vs)
+
+#distances from source
+d1=3000
+d2=5000
+
+#Number of samples and time increment of traces
+Nt=2000
+dt=1e-3
+
+#min and max frequencies to compute Q
+fqmin=5.0
+fqmax=20.0
+
+Tmax=dt*Nt
+time=[dt*i for i=1:Nt]
+
+println("Trace length is $Tmax s using $Nt samples")
+
+ind1=Int(trunc(d1/100))
+ind2=Int(trunc(d2/100))
+
+rec1=SeismicMatrix[ind1,2:Nt+1]
+rec2=SeismicMatrix[ind2,2:Nt+1]
+
+println("Trace 1 at offset $d1 m is found at index $ind1")
+println("Trace 2 at offset $d2 m is found at index $ind2")
+
+ampmin=-0.15
+ampmax=0.3
+
+p1=plot(time,rec1,xlabel="time (s)",ylabel="acc. (m.s⁻²)",ylim=(ampmin,ampmax),c=:"blue")
+p2=plot(time,rec2,xlabel="time (s)",ylabel="acc. (m.s⁻²)",ylim=(ampmin,ampmax),c=:"red")
+
+# second part rec1:
+println("Trace 1:")
+
+ampP,ampS,frequP,frequS=Spec(rec1,0.1,dt,Nt,0.4,0.75,0.35)
+
+ymax=maximum(ampP)
+
+p3=plot(frequP,ampP,title="P",xlabel="f(Hz)",ylabel="Amp",ylim=(0,ymax),c=:"blue")
+p4=plot(frequS,ampS,title="S",xlabel="f(Hz)",ylabel="Amp",ylim=(0,ymax),c=:"blue")
+
+#plot(p1,p3,p4,layout=(2,2))
+
+# second part rec2:
+println("Trace 2:")
+
+ampP2,ampS2,frequP2,frequS2=Spec(rec2,0.1,dt,Nt,0.7,1.25,0.35)
+
+p5=plot(frequP2,ampP2,title="P2",xlabel="f(Hz)",ylabel="Amp",ylim=(0,ymax),c=:"red")
+p6=plot(frequS2,ampS2,title="S2",xlabel="f(Hz)",ylabel="Amp",ylim=(0,ymax),c=:"red")
+
+
+@show frequP
+indf=findall(x->x>fqmin && x<fqmax,frequP)
+@show indf
+
+#Qp=zeros(length(indf))
+
+xp,yp=ComputeQgraph(ampP,ampP2,frequP,indf,d1,d2,Vp)
+
+p7=plot(xp,yp,xlabel="f(Hz)",ylim=(-1.0,0.),title="Qgraph for P phase",c=:"green")
+
+xs,ys=ComputeQgraph(ampS,ampS2,frequS,indf,d1,d2,Vs)
+
+@show xs,ys
+
+p8=plot(xs,ys,xlabel="f(Hz)",ylim=(-1.0,0.),title="Qgraph for S phase",c=:"green")
+
+
+plot(p1,p2,p3,p5,p4,p6,p7,p8,layout=(4,2),legend=false)
+#plot(p7,legend=false)
+
+end
+
+main()

examples/Wave1D_Sismo_Day4.jl

@@ -0,0 +1,128 @@
+using SeismicQ, Plots
+
+function MainSource()
+
+    # Spatial extent
+    Lx   = 50.0
+
+    # Mechanical parameters 
+    ρ₀   = 1500.0
+    K₀   = 1.e9
+    G₀   = 1.e8
+    c₀   = sqrt((K₀+4/3*G₀)/ρ₀) 
+
+    # Discretization
+    Ncx = 200
+    Δx  = Lx/Ncx
+    xv  = LinRange(0,Lx,Ncx+1)
+    xc  = LinRange(0-Δx/2,Lx+Δx/2,Ncx+2)
+
+    # Source parameters
+    𝑓₀   = 200     # Central frequency of the source [Hz]
+    t₀   = 1.2/𝑓₀
+    isrc = Int((Ncx/2)+1)
+
+    # Time domain
+    Δt   = min(1e10, Δx/c₀) # Courant criteria from wavespeed
+    Nt   = 100
+    Nout = 10
+    t    = 0.0#-t₀
+
+    # Parameters for Sismo.
+    Xs   = 25:0.5:50   # x_coordinates [m]
+    Ns   = size(Xs,1)
+    ds   = zeros(size(Xs))
+    @. ds   = abs(Xs-xv[isrc])
+    velocity_matrix = zeros(Ns, Nt)
+    time = zeros(Nt)
+
+    # Storage on centers # +2 for ghost nodes for BCs
+    szv   = (Ncx+1,)
+    szc   = (Ncx+2,)
+    # Storage on centroids 
+    K     = ones(szc)*K₀ 
+    G     = ones(szc)*G₀
+    ε̇     = ( xx=zeros(szc), yy=zeros(szc), zz=zeros(szc), xy=zeros(szc), yz=zeros(szc), xz=zeros(szc) )  
+    ∇V    = zeros(Ncx+2)
+    P     = zeros(Ncx+2)
+    τ     = ( xx=zeros(szc), yy=zeros(szc), zz=zeros(szc), xy=zeros(szc), yz=zeros(szc), xz=zeros(szc) )  
+    ∂Vx∂x = zeros(szc)
+    # Storage on vertices
+    V     = ( x=zeros(szv), y=zeros(szv), z=zeros(szv))
+    ρ     = ones(szv)*ρ₀ 
+    f_ext = zeros(szv)
+
+    # BC
+    Lbc        = 2.
+    bc_filtW_v = 1.0 .- exp.(-(xv.-0Lx).^2/Lbc.^2)
+    bc_filtW_c = 1.0 .- exp.(-(xc.-0Lx).^2/Lbc.^2)
+    bc_filtE_v = 1.0 .- exp.(-(xv.- Lx).^2/Lbc.^2)
+    bc_filtE_c = 1.0 .- exp.(-(xc.- Lx).^2/Lbc.^2)
+
+    # Time loop
+    @time for it=1:Nt
+
+        # Compute Ricker function
+        t          += Δt
+        a           = Ricker(t, t₀, 𝑓₀)
+        f_ext[isrc] = ρ[isrc]*a
+
+        # Velocity gradient components
+        @. ∂Vx∂x[2:end-1] = (V.x[2:end] - V.x[1:end-1])/Δx
+
+        # Divergence
+        @. ∇V   = ∂Vx∂x
+
+        # Deviatoric strain rate 
+        @. ε̇.xx = ∂Vx∂x - 1/3*∇V
+
+        # Stress update
+        @. τ.xx = f_shear(G)*Δt*(ε̇.xx) + f_relax(G)*τ.xx
+
+        # Pressure update 
+        @. P    = P - Δt*f_bulk(K)*∇V
+
+        # Linear momentum balance
+        @. V.x[2:end-1] = V.x[2:end-1] + Δt/ρ[2:end-1]*((τ.xx[3:end-1]-τ.xx[2:end-2])/Δx - (P[3:end-1]-P[2:end-2])/Δx - f_ext[2:end-1])
+
+        # Absorbing boundary Cerjean et al. (1985)
+        @.  V.x  = V.x  * bc_filtW_v 
+        @.  P    = P    * bc_filtW_c 
+        @.  τ.xx = τ.xx * bc_filtW_c 
+        @.  V.x  = V.x  * bc_filtE_v 
+        @.  P    = P    * bc_filtE_c 
+        @.  τ.xx = τ.xx * bc_filtE_c 
+
+        # Visualisation
+        if mod(it, Nout)==0
+            display(plot(xv, V.x, ylim=(-2e-4, 2e-4)))
+            sleep(0.1)
+        end
+
+        # Extract sismo data:
+        time[it] = t
+        @. velocity_matrix[:,it] = V.x[Int(Xs[:]/Δx)+1]
+
+    end
+
+    # Visualization Receiver Gather:
+    valim = max(abs(maximum(velocity_matrix)),abs(minimum(velocity_matrix)))
+    p = heatmap(ds,time,velocity_matrix',color=palette(:RdBu,100,rev=true),
+    title="Receiver gather", xlabel="distance to the source [m]",
+    ylabel="time [s]",yflip=true,clim=(-valim,+valim))
+    display(p)
+
+end
+
+function f_bulk(K) 
+   return K
+end
+
+function f_shear(G)
+    return 2*G
+end
+function f_relax(G)
+    return 1.
+end
+
+MainSource()

examples/Wave2D_Day3.jl

@@ -1,4 +1,4 @@
-using SeismicQ, Plots
+using SeismicQ, Plots,Classical
 
 function MainSource()
 
@@ -22,7 +22,7 @@ function MainSource()
 
 
     # Source parameters
-    𝑓₀   = 200     # Central frequency of the source [Hz]
+    𝑓₀   = 50     # Central frequency of the source [Hz]
     t₀   = 1.2/𝑓₀
     isrc = Int((Nc.x/2)+1)
     jsrc = Int((Nc.y/2)+1)
@@ -87,10 +87,10 @@ function MainSource()
         @. L.i.yx[:,2:end-1] = (V.c.y[2:end,2:end-1] - V.c.y[1:end-1,2:end-1])/Δ.x
         @. L.j.yx[2:end-1,:] = (V.v.y[2:end,:] - V.v.y[1:end-1,:])/Δ.x
 
-        @. L.i.yy[:,2:end-1] = (V.v.y[:,2:end] - V.v.y[:,2:end])/Δ.y
+        @. L.i.yy[:,2:end-1] = (V.v.y[:,2:end] - V.v.y[:,1:end-1])/Δ.y
         @. L.j.yy[2:end-1,:] = (V.c.y[2:end-1,2:end] - V.c.y[2:end-1,1:end-1])/Δ.y
 
-        @. L.i.xy[:,2:end-1] = (V.v.x[:,2:end] - V.v.x[:,2:end])/Δ.y
+        @. L.i.xy[:,2:end-1] = (V.v.x[:,2:end] - V.v.x[:,1:end-1])/Δ.y
         @. L.j.xy[2:end-1,:] = (V.c.x[2:end-1,2:end] - V.c.x[2:end-1,1:end-1])/Δ.y
 
 
@@ -149,14 +149,14 @@ function MainSource()
                                     *((τ.j.xy[3:end-1,2:end-1]-τ.j.xy[2:end-2,2:end-1])/Δ.x
                                     + (τ.i.yy[2:end-1,3:end-1]-τ.i.yy[2:end-1,2:end-2])/Δ.y 
                                     - (P.i[2:end-1,3:end-1]-P.i[2:end-1,2:end-2])/Δ.y 
-                                    - 0.0*f_ext.v[2:end-1,2:end-1]))
+                                    - f_ext.v[2:end-1,2:end-1]))
 
         @. V.c.y[2:end-1,2:end-1] = (V.c.y[2:end-1,2:end-1] 
                                     + Δt/ρ.c[2:end-1,2:end-1]
                                     *((τ.i.xy[2:end,2:end-1]-τ.i.xy[1:end-1,2:end-1])/Δ.x
                                     + (τ.j.yy[2:end-1,2:end]-τ.j.yy[2:end-1,1:end-1])/Δ.y 
                                     - (P.j[2:end-1,2:end]-P.j[2:end-1,1:end-1])/Δ.y 
-                                    - 0.0*Xf_ext.c[2:end-1,2:end-1]))   
+                                    - f_ext.c[2:end-1,2:end-1]))   
 
 
     #     # Absorbing boundary Cerjean et al. (1985)