import numpy as np
import dolfinx
from mpi4py import MPI
from petsc4py import PETSc
from ufl import TrialFunction, TestFunction, dot, grad, dx, inner
from dolfinx.fem import Function, functionspace, dirichletbc, locate_dofs_geometrical
from dolfinx.mesh import create_unit_square
from scipy.spatial.distance import cdist
from scipy.linalg import eigh
from dolfinx.fem.petsc import LinearProblem
from dolfinx import fem
from dolfinx import default_scalar_type
from dolfinx.fem import Function
from petsc4py import PETSc
from scipy.interpolate import RegularGridInterpolator
import numpy as np
from scipy.optimize import fsolve
import random
def solve_darcy(k_sample, f,mesh, V):
u_D = Function(V)
u_D.interpolate(lambda x:-1/(4*(np.pi**2))*np.sin(2*np.pi*x[0])-1/(4 *(np.pi**2))*np.cos(2*np.pi*x[1]))
left_boundary = locate_dofs_geometrical(V, lambda x: np.isclose(x[0], 0))
right_boundary = locate_dofs_geometrical(V, lambda x: np.isclose(x[0], 1))
bottom_boundary = locate_dofs_geometrical(V, lambda x: np.isclose(x[1], 0))
top_boundary = locate_dofs_geometrical(V, lambda x: np.isclose(x[1], 1))
bcs = [dirichletbc(u_D, left_boundary),
dirichletbc(u_D, right_boundary),
dirichletbc(u_D, bottom_boundary),
dirichletbc(u_D, top_boundary)]
dx_mesh = dx(domain=mesh)
u = TrialFunction(V)
v = TestFunction(V)
a = -dot(k_sample * grad(u), grad(v)) * dx_mesh
L = f * v * dx_mesh
u_sol= Function(V)
problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "cg", "pc_type": "hypre"})
u_sol = problem.solve()
return u_sol
from dolfinx.fem import assemble_scalar
import dolfinx
from mpi4py import MPI
from petsc4py import PETSc
from mpi4py import MPI
def compute_flux(u_sol, k_sample, mesh, V):
q = -k_sample * ufl.grad(u_sol)
facets = dolfinx.mesh.locate_entities_boundary(
mesh, mesh.topology.dim - 1, lambda x: np.isclose(x[0], 0.0)
)
if len(facets) == 0:
print("Warning: No facets found on the left boundary!")
facet_tags = dolfinx.mesh.meshtags(
mesh, mesh.topology.dim - 1, facets, np.full_like(facets, 1, dtype=np.int32)
)
n = ufl.FacetNormal(mesh)
ds_left = ufl.Measure("ds", domain=mesh, subdomain_data=facet_tags, subdomain_id=1)
flux_expr = ufl.dot(q, n) * ds_left
flux = mesh.comm.allreduce(fem.assemble_scalar(fem.form(flux_expr)), op=MPI.SUM)
return flux
def solve_transcendental_equation(lambda_value, num_modes):
def equation(w):
return np.tan(w) - (2 * lambda_value * w) / (lambda_value**2 * w**2 - 1)
w_n = [1.72066]
for n in range(1, num_modes ):
initial_guess = n * np.pi
w_n.append(fsolve(equation, initial_guess)[0])
return w_n
import numpy as np
import scipy
from dolfinx import fem
lam =0.5
dom = np.linspace(0, 1, 1000)
wn = solve_transcendental_equation(0.5, 100)
wn = np.array(wn)
thetan_1D = 2 * lam / (lam ** 2 * wn ** 2 + 1)
bn = np.zeros((len(dom), len(wn)))
An = np.zeros(len(wn))
for j in range(len(wn)):
for i in range(len(dom)):
bn[i, j] = np.sin(wn[j] * dom[i]) + lam * wn[j] * np.cos(wn[j] * dom[i])
An[j] = 1/np.sqrt(np.trapz(bn[:, j] ** 2, dom))
def generate_permeability_sample(norm_sample, wn, N_modes,mesh,V, sigma=1.0, lam=0.5 ):
wn = np.array(wn)
points = mesh.geometry.x
num_points = points.shape[0]
thetan_2D = np.zeros(N_modes)
i = 0
rowcolsum = 0
rownum = 0
colnum = rowcolsum
k_array = np.zeros(num_points)
while i < N_modes:
if rownum > rowcolsum:
rowcolsum += 1
rownum = 0
colnum = rowcolsum
thetan_2D[i] = thetan_1D[rownum] * thetan_1D[colnum]
scale_factor = norm_sample[i] * np.sqrt(thetan_2D[i])
b1 = An[rownum] * (np.sin(wn[rownum] * points[:, 0]) + lam * wn[rownum] * np.cos(wn[rownum] * points[:, 0]))
b2 = An[colnum] * (np.sin(wn[colnum] * points[:, 1]) + lam * wn[colnum] * np.cos(wn[colnum] * points[:, 1]))
k_array += scale_factor * b1 * b2
rownum += 1
colnum -= 1
i += 1
k_array = np.exp(k_array)
k_sample = fem.Function(V)
k_sample.x.array[:] = k_array
return k_sample
from dolfinx import geometry
import numpy as np
import ufl
lambda_value =0.5
sigma = 10**(-2)
N_modes = 20
np.random.seed(45)
mesh = create_unit_square(MPI.COMM_WORLD, 128, 128)
f =fem.Constant(mesh, default_scalar_type(-1))
V = functionspace(mesh, ("Lagrange", 1))
sol = np.random.normal(0, 1, N_modes)
lambda_value = 0.5
w_n = solve_transcendental_equation(lambda_value, N_modes)
k = generate_permeability_sample(sol, w_n, N_modes, mesh,V)
p_cor = solve_darcy(k, f,mesh,V)
sol_QOI = compute_flux(p_cor, k, mesh,V)
Sorry, I tested it and this is a reproducible example. I don’t think the issue is in the permeability but I think I fixed it by using a different way of defining my boundaries through which I compute the flux (see the darcy flux function). However, I am now wondering if there is some way to optimize the code since computing the solution to the darcy problem and computing the flux is a very expensive operation. Thank you so much!
