I am following from the methodology laid out in this paper to solve a special problem, and I am stuck on creating a vector field, \textbf{V}_h, such that it satisfies certain conditions. In order for \textbf{V}_h to be used in subsequent steps, it has to satisfy 2 conditions
\tag{5.13}\oint \textbf{V}_h \cdot \hat{\textbf{n}} ds = 0 i.e. the net flux on the boundary is 0, and
\tag{5.14}\oint \textbf{V}_h \cdot \hat{\textbf{s}}ds = 0 i.e. the net circulation on the boundary is 0.
If the obtained vector field \textbf{V}_h does not satisfy these conditions, then according to the paper, âsome adjustment must be made to the value of \textbf{V}_h so that the integral constraints (5.13) and (5.14) are satisfied.â And that gets us to what I want to implement. According to this user on stackexchange, I can try to convolute the divergence of my vector field with a Coulomb field in 2D, and the logic behind it is that it will create a vector field with the same divergence as the original vector field, and I can then subtract that vector field off of my \textbf{V}_h so that it becomes divergence-free, and then doing something similar to subsequently get rid of the curl. So, what I want to do, is take divV_h = ufl.div(V_h) and convolve it with a Coulomb field C = ufl.as_vector( (x[0]/(x[0]**2+x[1]**2), x[1]/(x[0]**2+x[1]**2)) ). Or I could create Dolfinx interpolated quantities of these ufl expressions and work with those, the point is that I want to see how to perform the convolution integrals between the scalar field with the vector field; the convolution will create a new vector field
\textbf{V}_{new} = \nabla \cdot \textbf{V}_h * \Big(\frac{x}{x^2 + y^2}\hat{\textbf{i}} + \frac{y}{x^2 + y^2}\hat{\textbf{j}}\Big)
. Additionally, it is stated that we can make this process much faster and simpler by taking advantage of convolution theorem by using a Fourier transform
F\Big\{\nabla \cdot \textbf{V}_h * \Big(\frac{x}{x^2+y^2}\hat{\textbf{i}} + \frac{y}{x^2+y^2}\hat{\textbf{j}}\Big)\Big\} = F\{\nabla \cdot \textbf{V}_h\} F\Big\{\Big(\frac{x}{x^2+y^2}\hat{\textbf{i}} + \frac{y}{x^2+y^2}\hat{\textbf{j}}\Big)\Big\}
as this lets us simply multiply the two transforms together rather than perform integration, and then we can subsequently inverse transform the product. I would appreciated it if someone could show me how to implement either the convolution, or taking the Fourier Transform of these variables.
To preface my code a bit, it is following the methodology in the paper (where I have excluded the iterative procedure part as it is not relevant to this question), except instead of using the Harmonic-Sine method, I am using Dolfinx for everything involving solving PDEs. The formula solved in the code is \textbf{V}_h = \textbf{V} - \textbf{V}_I, and the code is used to compute \textbf{V}_I by first solving 2 Poisson Equations with homogeneous boundary conditions to obtain \psi_I and \chi_I, and \textbf{V}_I = \textbf{V}_{\psi_I} + \textbf{V}_{\chi_I}= \hat{\textbf{k}}\times \nabla \psi_I + \nabla \chi_I. Here is my code for the creation of \textbf{V}_h given a vector field \textbf{V} = 3(x-y)\hat{\textbf{i}} + (3x-y)\hat{\textbf{j}} on the domain R:[-4\le x \le 4, -4\le y \le 4]
from mpi4py import MPI
import pyvista as pv
import dolfinx as df
import ufl
from petsc4py import PETSc
import numpy as np
import matplotlib.pyplot as plt
Nx,Ny = 60,60 # Number of points in x and y
family = "Lagrange"
order = 1
# overlapping_facets = True # Uncertainty of whether or not degrees of freedom should overlap when setting boundary conditions for each boundary
Lx, Ly = 4,4
mesh = df.mesh.create_rectangle(MPI.COMM_WORLD, [[-Lx,-Ly],[Lx,Ly]] ,[Nx, Ny])
tdim = mesh.topology.dim
fdim = tdim-1
mesh.topology.create_connectivity(tdim-1, tdim)
# Define a true functionspace for later usage
S = df.fem.functionspace(mesh, (family, order)) # scalar space for grids and projection of scalar quantities
# Define total wind field
V = df.fem.functionspace(mesh, (family, order, (2,))) # 2D Vector function space
wind = df.fem.Function(V)
wind.interpolate(lambda x: (3*(x[0]-x[1]),3*x[0]-x[1]))
# Extract vertical relative vorticity from wind field (onto scalar space IMPORTANT)
Z_expr = df.fem.Expression(wind[1].dx(0) - wind[0].dx(1),S.element.interpolation_points)
Z = df.fem.Function(S)
Z.interpolate(Z_expr) # Evaluate vort epression
# Define boundary condition
boundary_facets = df.mesh.exterior_facet_indices(mesh.topology)
boundary_dofs = df.fem.locate_dofs_topological(S, tdim-1, boundary_facets)
bc = df.fem.dirichletbc(df.default_scalar_type(0),boundary_dofs,S)
u = ufl.TrialFunction(S)
v = ufl.TestFunction(S)
a = ufl.inner(ufl.grad(u),ufl.grad(v)) * ufl.dx
# f = Vorticity
L = -Z*v*ufl.dx
from dolfinx.fem.petsc import LinearProblem
problem = LinearProblem(a=a,L=L,bcs = [bc], petsc_options_prefix="Streamfunction_I",petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
psi_I = problem.solve()
# Extract horizontal divergence from wind field (onto scalar space)
d_expr = df.fem.Expression(ufl.div(wind),S.element.interpolation_points)
d = df.fem.Function(S)
d.interpolate(d_expr) # Evaluate div epression
# f = Divergence
L = -d*v*ufl.dx
problem = LinearProblem(a=a,L=L,bcs = [bc], petsc_options_prefix="Velocity_potential_I",petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
chi_I = problem.solve()
# Move to 2D vector space to store V_I into
V_I = df.fem.Function(V)
V_I_expr = df.fem.Expression(ufl.perp(ufl.grad(psi_I)) + ufl.grad(chi_I), V.element.interpolation_points)
V_I.interpolate(V_I_expr)
V_h_true = df.fem.Function(V) # Unique harmonic wind vector
V_h_expr = df.fem.Expression(wind - V_I, V.element.interpolation_points)
V_h_true.interpolate(V_h_expr)
# V_h = V_h_true.copy() # This will be modified in the iterative procedure.
V_h_values = np.zeros([V.dofmap.index_map.size_local, 3])
V_h_values[:,:mesh.geometry.dim] = V_h_true.x.array.reshape(V.dofmap.index_map.size_local, V.dofmap.index_map_bs)
grid = pv.UnstructuredGrid(*df.plot.vtk_mesh(V))
grid["V_h"] = V_h_values
glyphs = grid.glyph(orient = "V_h", scale = "V_h", factor = .075, tolerance=0.025)
grid.set_active_vectors("V_h")
plotter = pv.Plotter()
plotter.add_mesh(grid, show_edges= False, color = 'white')
plotter.add_mesh(glyphs)
plotter.view_xy()
plotter.show()
where we want to modify the vector field V_h_true seen plotted in the code so that it then becomes divergence-free and curl-free. Additionally, here is the code for computing the net flux and circulation as well as the individual net fluxes and circulations along each of the 4 boundaries
boundaries = [
(1, lambda x: np.isclose(x[1], -Ly)),
(2, lambda x: np.isclose(x[0], Lx)),
(3, lambda x: np.isclose(x[1], Ly)),
(4, lambda x: np.isclose(x[0], -Lx)),
]
facet_indices, facet_markers = [], []
fdim = tdim - 1
for marker, locator in boundaries:
facets = df.mesh.locate_entities(mesh, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = df.mesh.meshtags(
mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets]
)
ds = ufl.Measure("ds", domain=mesh, subdomain_data=facet_tag)
boundary_Vh_flux = []
boundary_Vh_circ = []
boundary_psih_flux = []
for boundary in boundaries:
boundary_marker = boundary[0]
compiled_form_n = df.fem.form(ufl.dot(ufl.FacetNormal(mesh),V_h)* ds(boundary_marker))
bf = mesh.comm.allreduce(df.fem.assemble_scalar(compiled_form_n), op = MPI.SUM)
boundary_Vh_flux.append(bf)
compiled_form_t = df.fem.form(ufl.dot(ufl.perp(ufl.FacetNormal(mesh)),V_h)* ds(boundary_marker))
bc = mesh.comm.allreduce(df.fem.assemble_scalar(compiled_form_t), op = MPI.SUM)
boundary_Vh_circ.append(bc)
#compiled_form = df.fem.form(ufl.dot(ufl.FacetNormal(mesh),ufl.grad(psi_h_0))* ds(boundary_marker))
#psih_f = mesh.comm.allreduce(df.fem.assemble_scalar(compiled_form), op = MPI.SUM)
#boundary_psih_flux.append(psih_f)
print("fluxes of V_h along 4 boundaries = ", boundary_Vh_flux)
print("net V_h flux = ", sum(boundary_Vh_flux))
print("circulation of $V_h$ along 4 boundaries = ", boundary_Vh_circ)
print("net V_h circulation = ", sum(boundary_Vh_circ))
#print("fluxes of $psi_h$ along 4 boundaries = ", boundary_psih_flux)
#print("net psi_h flux = ", sum(boundary_psih_flux))
where we see that it fails to satisfy both conditions for the net flux and circulation conditions. Once we are finished modifying V_h by subtracting off the convolution (or convolution theorem using Fourier Transform) of the divergence and Coulomb field, we want the resulting vector field to be divergence-free resulting in net flux to approximately vanish. Then, performing the same procedure to make the vector field curl-free should cause both the net flux and net circulation to vanish.
Just as an aside, something else I was considering was modifying the variational form to incorporate these net flux and circulation conditions, but it was a bit too advanced for me to figure out how to create the correct form. If I could incorporate these as conditions, then V_h would satisfy these conditions from the get-go.
