Hi,
I’ve been trying to deal with the normalization of a vector \boldsymbol{a} that depends on a discontinuous tensor field \boldsymbol{\gamma} through \boldsymbol{a} = \boldsymbol{X} : \boldsymbol{\gamma}, a_i = \varepsilon_{ijk} \gamma_{jk}, where \boldsymbol{X} is the third-rank permutation tensor. In my specific problem, the initial distribution of \boldsymbol{\gamma} in my domain is \gamma_{13} = 1.0 inside a square in the center, as the figure shows.
The vector \boldsymbol{a} bears a relation with the transport velocity of the field \boldsymbol{\gamma}, whose evolution is given by a transport equation which is solved by the FEM in FEniCS. At each step, the new position of \boldsymbol{\gamma} is calculated, and the resolution introduces some non-zero numerical approximations around the initial square. These approximations results in some spurious values of the field \boldsymbol{a} around the \boldsymbol{\gamma} distribution, as can be seen in the figure below (the values above the rectangle)
As a first attempt to avoid these instabilities due to the sharp changes in the solution of \boldsymbol{\gamma}, I’ve tried to control the normalization of \boldsymbol{a} through some tolerance values, as can be seen in the MWE below.
import dolfin as dlf
import numpy as np
Lbox = 15.0
Hbox = 15.0
Wbox = 100000.0
b = 0.25
r_c = np.sqrt(b)
al_in_X_1 = Lbox/2.
al_in_Y_1 = Hbox/2.
al_in_X_2 = 3./4.*Lbox
al_in_Y_2 = Hbox/2.
gamma_init_1 = dlf.Expression(
"fabs(x[0]-al_in_X) <= rc & fabs(x[1]-al_in_Y) <= rc ? b/(pow(rc, 2)) : 0.",
rc=r_c,
b=b,
al_in_X=al_in_X_1,
al_in_Y=al_in_Y_1,
degree=1,
)
N_el_x, N_el_y, N_el_z = 100, 100, 1
mesh = dlf.BoxMesh(dlf.Point(0.0, 0.0, 0.0), dlf.Point(Lbox, Hbox, Wbox), N_el_x, N_el_y, N_el_z)
V_gamma = dlf.TensorFunctionSpace(mesh, "CG", 1)
V_vec = dlf.VectorFunctionSpace(mesh, "CG", 1)
gamma = dlf.Function(V_gamma, name="gamma")
gamma_init = dlf.Expression(
[
['0', '0', 'gamma_init'],
['0', '0', '0'],
['0', '0', '0']
],
gamma_init=gamma_init_1,
degree=1
)
gamma.assign(dlf.interpolate(gamma_init, V_gamma))
# Define permutation tensor
Id = np.eye(3)
A = lambda *ind: np.linalg.det(np.array([[Id[p, q] for q in range(3)] for p in ind]))
perm_np = np.array(
[
[
[
A(i, j, k) for k in range(3)
] for j in range(3)
]
for i in range(3)
]
)
gamma_vec = gamma.vector()[:]
num_nodes = int(len(gamma_vec) / 9.)
gamma_nodal = gamma_vec.reshape((num_nodes, 3, 3))
a = np.zeros((num_nodes, 3))
a_normalized = np.zeros_like(a)
# Tolerance values to handle the numerical instabilities
tol_a = 2e-2
tol_norm_a = 1e-10
# Get vector values at each node
for ii in range(num_nodes):
a[ii] = np.einsum('ijk,jk', perm_np, gamma_nodal[ii,:,:])
# Check for small values
if np.linalg.norm(a[ii]) < tol_a:
continue
# Normalize the vector
if np.linalg.norm(a[ii]) > tol_norm_a:
a_normalized[ii] = a[ii] / (np.linalg.norm(a[ii]))
a_function = dlf.Function(V_vec, name="a")
a_function.vector()[:] = a.flatten()[:]
a_normalized_function = dlf.Function(V_vec, name="a_normalized")
a_normalized_function.vector()[:] = a_normalized.flatten()[:]
file_results = dlf.XDMFFile("MWE.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True
file_results.write(gamma, 0.)
file_results.write(a_function, 0.)
file_results.write(a_normalized_function, 0.)
Could you give me some ideas on how to better perform the normalization of \boldsymbol{a} in order to handle the fact that \boldsymbol{\gamma} varies too sharply in the domain?
Thanks in advance.
