From 79cae09a66a76fa2f42ff504a0837af6928940ee Mon Sep 17 00:00:00 2001 From: Daniel O'Malley Date: Sun, 29 Sep 2024 23:55:51 -0600 Subject: [PATCH] work on 3d fracture example --- .../monte_carlo_fracture_3D.jl | 15 +- .../weeds_fracture_3D.jl | 138 ++++++------------ 2 files changed, 52 insertions(+), 101 deletions(-) diff --git a/examples/pressure_management_inj/monte_carlo_fracture_3D.jl b/examples/pressure_management_inj/monte_carlo_fracture_3D.jl index 75f2fc0..96935ac 100644 --- a/examples/pressure_management_inj/monte_carlo_fracture_3D.jl +++ b/examples/pressure_management_inj/monte_carlo_fracture_3D.jl @@ -2,26 +2,27 @@ import Random Random.seed!(1) include("weeds_fracture_3D.jl") x_test = randn(3 * num_eigenvectors + 2) #the input to the function should follow an N(0, I) distribution -- the eigenvectors and eigenvalues are embedded inside f -@show f(x_test) #compute the pressure at the critical location -@show extrema(grad_f(x_test)) #compute the gradient of the pressure at the critical location with respect to the eigenvector coefficients (then only show the largest/smallest values so it doesn't print a huge array) +@time @show f(x_test) #compute the pressure at the critical location +@time @show extrema(grad_f(x_test)) #compute the gradient of the pressure at the critical location with respect to the eigenvector coefficients (then only show the largest/smallest values so it doesn't print a huge array) #Quick demo of how you could do Monte Carlo to get the distribution of pressures at the critical location -- hopefully the method Georg described can help us understand the tail distribution better -numruns = 10 ^ 3 +numruns = 10 ^ 2 fs = Float64[] @time for i = 1:numruns x_test = randn(3 * num_eigenvectors + 2) - push!(fs, f(x_test)) + @time push!(fs, f(x_test)) @show i end @show sum(fs) / length(fs) @show extrema(fs) -fig, axs = PyPlot.subplots(1, 2) +fig, axs = PyPlot.subplots(1, 2, figsize=(16, 9)) axs[1].hist(fs, bins=30) axs[1].set(xlabel="Pressure [MPa] at critical location", ylabel="Count (out of $(numruns))") axs[2].hist(log.(fs), bins=30) axs[2].set(xlabel="log(Pressure [MPa]) at critical location", ylabel="Count (out of $(numruns))") -#display(fig); println() -#fig.show() +fig.tight_layout() +display(fig); println() +fig.show() fig.savefig("distribution_of_pressure.png") PyPlot.close(fig) diff --git a/examples/pressure_management_inj/weeds_fracture_3D.jl b/examples/pressure_management_inj/weeds_fracture_3D.jl index 96c3e47..b60d43a 100644 --- a/examples/pressure_management_inj/weeds_fracture_3D.jl +++ b/examples/pressure_management_inj/weeds_fracture_3D.jl @@ -3,7 +3,7 @@ import GaussianRandomFields import PyPlot import Zygote -plotstuff = false +plotstuff = true module FractureMaker mutable struct FractureTip @@ -37,18 +37,24 @@ end function timestep!(m::Material, tips::Array{FractureTip}) for tip in tips if inside(m, tip) && tip.alive - if tip.di ==0 + if tip.di == 0 + m.isfractured[tip.i, tip.j, :] .= true + #= m.isfractured[tip.i, tip.j, tip.k] = true m.isfractured[tip.i, tip.j, tip.k+1] = true m.isfractured[tip.i, tip.j, tip.k-1] = true m.isfractured[tip.i, tip.j, tip.k+2] = true m.isfractured[tip.i, tip.j, tip.k-2] = true + =# else + m.isfractured[tip.i, tip.j, :] .= true + #= m.isfractured[tip.i, tip.j, tip.k] = true m.isfractured[tip.i, tip.j, tip.k+1] = true m.isfractured[tip.i, tip.j, tip.k-1] = true m.isfractured[tip.i, tip.j, tip.k+2] = true m.isfractured[tip.i, tip.j, tip.k-2] = true + =# end end if tip.alive @@ -71,7 +77,7 @@ end end # module FractureMaker # Set up the mesh -ns = [101, 101, 50]#mesh size -- make it at least 100x200 to have the resolution for the fractures +ns = [101, 101, 201]#mesh size steadyhead = 0e0 #meters width = 1000 #meters heights = [200, 100, 200] @@ -91,30 +97,9 @@ for i = 1:size(coords, 2) end masks = map(x->reshape(x, reverse(ns)...), masks) -masksPlot = zeros(prod(ns),1) - -for i = 1:size(coords, 2) - if coords[3, i] > heights[1] + heights[2] - masksPlot[i] = 1 - elseif coords[3, i] > heights[1] - masksPlot[i] =2 - else - masksPlot[i] = 3 - end -end -# Create a new figure for the 3D plot -if plotstuff - fig = PyPlot.figure() - ax = fig[:add_subplot](111, projection="3d") - scatter_plot = ax[:scatter](coords[1,:],coords[2,:],coords[3,:], c=masksPlot, cmap="viridis") - display(fig) - println() - PyPlot.close(fig) -end #setup the fracture network material = FractureMaker.Material(reverse(ns)...) #Logical array size equal to total number of cell - -true_ind = findall( masks[2]) +true_ind = findall(masks[2]) true_ks=[] for i=1:length(true_ind) global true_ks=push!(true_ks,true_ind[i][1]) @@ -126,8 +111,8 @@ for frac = 1:num_horizontal_fractures i = rand(true_ks) j = rand(1:ns[2]) k = rand(1:ns[1]) - tips[2 * frac - 1] = FractureMaker.FractureTip(i,j,k, 0, 1, 0, true) - tips[2 * frac] = FractureMaker.FractureTip(i,j,k, 0, -1, 0, true) + tips[2 * frac - 1] = FractureMaker.FractureTip(i, j, k, 0, 1, 0, true) + tips[2 * frac] = FractureMaker.FractureTip(i, j, k, 0, -1, 0, true) end FractureMaker.simulate(material, tips, 9 * ns[1] ÷ 20)#the time steps let the fracture grow to be at most a little less than the width of the domain horizontal_fracture_mask = copy(material.isfractured) .* masks[2] @@ -138,47 +123,20 @@ for frac = 1:num_vertical_fractures i = rand(true_ks) j = rand(1:ns[2]) k = rand(1:ns[1]) - tips[2*frac-1] = FractureMaker.FractureTip(i,j,k,1,0,0, true) - tips[2*frac] = FractureMaker.FractureTip(i, j,k,-1,0,0, true) + tips[2 * frac - 1] = FractureMaker.FractureTip(i, j, k, 1, 0, 0, true) + tips[2 * frac] = FractureMaker.FractureTip(i, j, k, -1, 0, 0, true) end FractureMaker.simulate(material, tips, ns[3] * heights[2] ÷ (3 * sum(heights))) vertical_fracture_mask = copy(material.isfractured) .* masks[2] - horizontal_fracture_mask -horizontal_fracture_mask_plot=reshape(horizontal_fracture_mask,prod(ns)) -vertical_fracture_mask_plot=reshape(vertical_fracture_mask,prod(ns)) +horizontal_fracture_mask_plot=reshape(horizontal_fracture_mask, prod(ns)) +vertical_fracture_mask_plot=reshape(vertical_fracture_mask, prod(ns)) -if plotstuff - fig = PyPlot.figure() - ax = fig[:add_subplot](111, projection="3d") - using Colors - category_colors = [ - [1, 0, 0, 0.0], - [0, 1, 0, 1.0], - ] - colors = [category_colors[c+1] for c in Int.(horizontal_fracture_mask_plot)] - scatter_plot = ax[:scatter](coords[1,:],coords[2,:],coords[3,:], c=colors) - display(fig) - println() - PyPlot.close(fig) - - - category_colors = [ - [0, 0, 0, 0.0], - [1, 0, 0, 1.0], - ] - fig = PyPlot.figure() - ax = fig[:add_subplot](111, projection="3d") - colors = [category_colors[c+1] for c in Int.(vertical_fracture_mask_plot)] - scatter_plot = ax[:scatter](coords[1,:],coords[2,:],coords[3,:], c=colors) - display(fig) - println() - PyPlot.close(fig) -end # Set up the eigenvector parameterization of the geostatistical log-permeability field num_eigenvectors = 200 sigma = 1.0 lambda = 50 -mean_log_conductivities = [log(1e-4), log(1e-8), log(1e-3)] #log(1e-4 [m/s]) -- similar to a pretty porous oil reservoir +mean_log_conductivities = [log(1e-3), log(1e-8), log(1e-3)] #log(1e-4 [m/s]) -- similar to a pretty porous oil reservoir cov = GaussianRandomFields.CovarianceFunction(2, GaussianRandomFields.Matern(lambda, 1; σ=sigma)) x_pts = range(mins[1], maxs[1]; length=ns[1]) y_pts = range(mins[2], maxs[2]; length=ns[2]) @@ -186,37 +144,26 @@ grf = GaussianRandomFields.GaussianRandomField(cov, GaussianRandomFields.Karhune parameterization = copy(grf.data.eigenfunc) sigmas = copy(grf.data.eigenval) - logKs = zeros(reverse(ns)...) parameterization3D=zeros(ns[3],ns[1]*ns[2],num_eigenvectors) sigmas3D=zeros(ns[3],num_eigenvectors) logKs2d = GaussianRandomFields.sample(grf)'#generate a random realization of the log-conductivity field for i = 1:ns[3]#copy the 2d field to each of the 3d layers - if i< true_ks[1] - j=1 - elseif i>true_ks[end] - j=3 + if i < true_ks[1] + j = 1 + elseif i > true_ks[end] + j = 3 else - j=2 + j = 2 end - t=view(parameterization3D, i, :,:) - t.=parameterization - m=view(sigmas3D, i, :) - m.=sigmas + t = view(parameterization3D, i, :,:) + t .= parameterization + m = view(sigmas3D, i, :) + m .= sigmas v = view(logKs, i, :, :) v .= mean_log_conductivities[j].+logKs2d end -if plotstuff - #plot the log-conductivity - fig, ax = PyPlot.subplots() - img = ax.imshow(logKs[:, 50, :], origin="lower") - ax.title.set_text("Conductivity Field") - fig.colorbar(img) - display(fig) - println() - PyPlot.close(fig) -end #play with these four numbers to make the fractures more or less likely to enable lots of flow from the lower reservoir to the upper reservoir horizontal_fracture_mean_conductivity = log(1e-10) vertical_fracture_mean_conductivity = log(1e-10) @@ -231,13 +178,14 @@ function x2logKs(x) end logKs = x2logKs(randn(3 * num_eigenvectors + 2)) if plotstuff - for i=1:101 - fig, axs = PyPlot.subplots(1, 1) - im2=axs.imshow(logKs[:, i, :], origin="lower") - display(fig); println() - PyPlot.close(fig) - @show i - end + #plot the log-conductivity + fig, ax = PyPlot.subplots() + img = ax.imshow(logKs[:, 50, :], origin="lower") + ax.title.set_text("Conductivity Field with fractures") + fig.colorbar(img) + display(fig) + println() + PyPlot.close(fig) end # Make the boundary conditions be dirichlet with a fixed head on the left and right boundaries with zero flux on the top and bottom dirichletnodes = Int[] @@ -257,7 +205,7 @@ function get_nodes_near(coords, obslocs) return obsnodes end -injection_node, critical_point_node = get_nodes_near(coords, [[width / 2,width / 2, heights[1] / 2], [width / 2, width / 2, heights[1] + heights[2] + 0.5 * heights[3]]]) #put the critical point near (-80, -80) and the injection node near (80, 80) +injection_node, critical_point_node = get_nodes_near(coords, [[width / 2,width / 2, 4 * heights[1] / 5], [width / 2, width / 2, heights[1] + heights[2] + 0.1 * heights[3]]]) injection_nodes = [injection_node] # Set the Qs for the whole field to zero, then set it to the injection rate divided by the number of injection nodes at the injection nodes @@ -270,7 +218,7 @@ function solveforh(logKs) #@show length(logKs) @assert length(logKs) == length(Qs) Ks_neighbors = logKs2Ks_neighbors(logKs) - result = DPFEHM.groundwater_steadystate(Ks_neighbors, neighbors, areasoverlengths, dirichletnodes, dirichleths, Qs; linear_solver=DPFEHM.cholesky_solver) + result = DPFEHM.groundwater_steadystate(Ks_neighbors, neighbors, areasoverlengths, dirichletnodes, dirichleths, Qs; linear_solver=DPFEHM.amg_solver) return result end @@ -282,12 +230,15 @@ grad_f(x) = Zygote.gradient(f, x)[1] if plotstuff import Random Random.seed!(1) + x_test = randn(3 * num_eigenvectors + 2) #the input to the function should follow an N(0, I) distribution -- the eigenvectors and eigenvalues are embedded inside f - @show x_test[1:3] - @show f(x_test) #compute the pressure at the critical location + @show x_test[[1, end - 1, end]] + #@show f(x_test) #compute the pressure at the critical location #compute the gradient of the pressure at the critical location with respect to the eigenvector coefficients (then only show the largest/smallest values so it doesn't print a huge array) - @show extrema(grad_f(x_test)) - P_mat=reshape(solveforh(x2logKs(x_test)), reverse(ns)...) + #@show extrema(grad_f(x_test)) + + @time P_mat=reshape(solveforh(x2logKs(x_test)), reverse(ns)...) + @show P_mat[critical_point_node] fig, ax = PyPlot.subplots() img = ax.imshow(P_mat[:, :,50], origin="lower") ax.title.set_text("pressure Field") @@ -295,5 +246,4 @@ if plotstuff display(fig) println() PyPlot.close(fig) - end - +end