search

FEniCS / dolfinx

Next generation FEniCS problem solving environment
https://fenicsproject.org
GNU Lesser General Public License v3.0
776 stars 181 forks source link

[BUG]: Erroneous data in Expression interpolation #3040

Closed nate-sime closed 9 months ago

nate-sime commented 9 months ago

Summarize the issue

Attempting to interpolate the following UFL form yields unexpected data.

How to reproduce the bug

Consider the following example. I'd expect the output to yield In the first four tests: result: (1.0, 1.0) In the second four tests: result: (100.0, 100.0)

Further removing the variable a yields correct data.

Minimal Example (Python)

import dolfinx
import ufl
from mpi4py import MPI

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 1, 1)

V = dolfinx.fem.FunctionSpace(mesh, ("CG", 1))
p = dolfinx.fem.Function(V)
p.interpolate(lambda x: x[0] + x[1])

left = dolfinx.fem.Constant(mesh, (0.0, 0.0))
right = dolfinx.fem.Constant(mesh, (0.0, 0.0))
a = dolfinx.fem.Constant(mesh, 1.0)

q = dolfinx.fem.Function(dolfinx.fem.FunctionSpace(
    mesh, ("DG", 0, (2,))))
expression = dolfinx.fem.Expression(
    a * ufl.grad(p - ufl.dot(left, right)),
    q.function_space.element.interpolation_points())

def try_values(vals):
    left.value[:] = vals
    right.value[:] = vals
    q.interpolate(expression)
    print(f"trying: {vals}, result: {q.x.array.min(), q.x.array.max()}")

print(f"Now setting a.value = {a.value}")
try_values((0.0, 0.0))
try_values((0.0, 1.0))
try_values((1.0, 0.0))
try_values((1.0, 1.0))

a.value = 100.0
print(f"Now setting a.value = {a.value}")
try_values((0.0, 0.0))
try_values((0.0, 1.0))
try_values((1.0, 0.0))
try_values((1.0, 1.0))

Output (Python)

Now setting a.value = 1.0
trying: (0.0, 0.0), result: (0.0, 0.0)
trying: (0.0, 1.0), result: (0.0, 0.0)
trying: (1.0, 0.0), result: (1.0, 1.0)
trying: (1.0, 1.0), result: (1.0, 1.0)
Now setting a.value = 100.0
trying: (0.0, 0.0), result: (0.0, 0.0)
trying: (0.0, 1.0), result: (0.0, 0.0)
trying: (1.0, 0.0), result: (1.0, 1.0)
trying: (1.0, 1.0), result: (1.0, 1.0)

Version

0.7.3

DOLFINx git commit

25db9a743ef95a78f686f05e38607aa18281b0d6

Installation

From source.

Additional information

No response

jorgensd commented 9 months ago

The generated code messes up the constant numbering, as the constant inside the gradient should be zeroed. Thus, the numbering is wrong:

import dolfinx
import ufl
from mpi4py import MPI
from pathlib import Path

cache_dir = f"{str(Path.cwd())}/.cache"
jit_options = {"cffi_extra_compile_args": ["-O3", "-march=native"],
               "cache_dir": cache_dir, "cffi_libraries": ["m"]}

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 1, 1)

V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))
p = dolfinx.fem.Function(V)
p.interpolate(lambda x: x[0] + x[1])

left = dolfinx.fem.Constant(mesh, (0.0, 0.0))
a = dolfinx.fem.Constant(mesh, 1.0)

q = dolfinx.fem.Function(dolfinx.fem.functionspace(
    mesh, ("DG", 0, (2,))))
expression = dolfinx.fem.Expression(
    a * ufl.grad(p - ufl.dot(left, left)),
    q.function_space.element.interpolation_points(), jit_options=jit_options)

def try_values(vals):
    left.value[:] = vals
    q.interpolate(expression)
    print(f"trying: {vals}, result: {q.x.array.min(), q.x.array.max()}")

print(f"Now setting a.value = {a.value}")
try_values((0.0, 0.0))
try_values((0.0, 1.0))
try_values((1.0, 0.0))
try_values((1.0, 1.0))

a.value = 100.0
print(f"Now setting a.value = {a.value}")
try_values((0.0, 0.0))
try_values((0.0, 1.0))
try_values((1.0, 0.0))
try_values((1.0, 1.0))

gives generated code

void tabulate_tensor_expression_a6098be59ff4a4e76a22744d94f352be9ee2c0c3(double* restrict A,
                                    const double* restrict w,
                                    const double* restrict c,
                                    const double* restrict coordinate_dofs,
                                    const int* restrict entity_local_index,
                                    const uint8_t* restrict quadrature_permutation)
{
// Precomputed values of basis functions
// FE* dimensions: [entities][points][dofs]
static const double FE0_C0_D10_Q083[1][1][1][3] = {{{{-1.0, 1.0, 0.0}}}};
static const double FE1_C1_D01_Q083[1][1][1][3] = {{{{-1.0, 0.0, 1.0}}}};
// Unstructured piecewise computations
double w0_d10 = 0.0;
for (int ic = 0; ic < 3; ++ic)
{
  w0_d10 += w[ic] * FE0_C0_D10_Q083[0][0][0][ic];
}
double J_c3 = 0.0;
for (int ic = 0; ic < 3; ++ic)
{
  J_c3 += coordinate_dofs[(ic) * 3 + 1] * FE1_C1_D01_Q083[0][0][0][ic];
}
double J_c0 = 0.0;
for (int ic = 0; ic < 3; ++ic)
{
  J_c0 += coordinate_dofs[(ic) * 3] * FE0_C0_D10_Q083[0][0][0][ic];
}
double J_c1 = 0.0;
for (int ic = 0; ic < 3; ++ic)
{
  J_c1 += coordinate_dofs[(ic) * 3] * FE1_C1_D01_Q083[0][0][0][ic];
}
double J_c2 = 0.0;
for (int ic = 0; ic < 3; ++ic)
{
  J_c2 += coordinate_dofs[(ic) * 3 + 1] * FE0_C0_D10_Q083[0][0][0][ic];
}
double w0_d01 = 0.0;
for (int ic = 0; ic < 3; ++ic)
{
  w0_d01 += w[ic] * FE1_C1_D01_Q083[0][0][0][ic];
}
double sp_0 = J_c0 * J_c3;
double sp_1 = J_c1 * J_c2;
double sp_2 = -sp_1;
double sp_3 = sp_0 + sp_2;
double sp_4 = J_c3 / sp_3;
double sp_5 = w0_d10 * sp_4;
double sp_6 = -J_c2;
double sp_7 = sp_6 / sp_3;
double sp_8 = w0_d01 * sp_7;
double sp_9 = sp_5 + sp_8;
double sp_10 = c[0] * sp_9;
double sp_11 = J_c0 / sp_3;
double sp_12 = w0_d01 * sp_11;
double sp_13 = -J_c1;
double sp_14 = sp_13 / sp_3;
double sp_15 = w0_d10 * sp_14;
double sp_16 = sp_12 + sp_15;
double sp_17 = c[0] * sp_16;
for (int iq = 0; iq < 1; ++iq)
{
  A[2 * iq + 0] += sp_10;
  A[2 * iq + 1] += sp_17;
}

}

where we only use c[0]. but c[0] in our code is the left constant's first value.

jorgensd commented 9 months ago

Slightly more compact example:

import dolfinx
import ufl
from mpi4py import MPI
from pathlib import Path

cache_dir = f"{str(Path.cwd())}/.cache"
jit_options = {"cffi_extra_compile_args": ["-O3", "-march=native"],
               "cache_dir": cache_dir, "cffi_libraries": ["m"]}

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 1, 1)

V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))
p = dolfinx.fem.Function(V)
p.interpolate(lambda x: 2*x[0])

b = dolfinx.fem.Constant(mesh, 10.0)
a = dolfinx.fem.Constant(mesh, 1.0)
q = dolfinx.fem.Function(dolfinx.fem.functionspace(
    mesh, ("DG", 0)))
expression = dolfinx.fem.Expression(
    a * ufl.Dx(p - b, 0),
    q.function_space.element.interpolation_points(), jit_options=jit_options)

def try_values(val):
    b.value = val
    q.interpolate(expression)
    print(f"trying: {val}, result: {q.x.array.min(), q.x.array.max()}")

print(f"Now setting a.value = {a.value}")
try_values(3.0)
try_values(2.0)

a.value = 100.0
print(f"Now setting a.value = {a.value}")
try_values(3.0)
try_values(2.0)
Now setting a.value = 1.0
trying: 3.0, result: (6.0, 6.0)
trying: 2.0, result: (4.0, 4.0)
Now setting a.value = 100.0
trying: 3.0, result: (6.0, 6.0)
trying: 2.0, result: (4.0, 4.0)

Generated code

void tabulate_tensor_expression_0d130e2ee62787808de92693faef788adb9b77d2(double* restrict A,
                                    const double* restrict w,
                                    const double* restrict c,
                                    const double* restrict coordinate_dofs,
                                    const int* restrict entity_local_index,
                                    const uint8_t* restrict quadrature_permutation)
{
// Precomputed values of basis functions
// FE* dimensions: [entities][points][dofs]
static const double FE0_C0_D10_Q083[1][1][1][3] = {{{{-1.0, 1.0, 0.0}}}};
static const double FE1_C1_D01_Q083[1][1][1][3] = {{{{-1.0, 0.0, 1.0}}}};
// Unstructured piecewise computations
double w0_d10 = 0.0;
for (int ic = 0; ic < 3; ++ic)
{
  w0_d10 += w[ic] * FE0_C0_D10_Q083[0][0][0][ic];
}
double J_c3 = 0.0;
for (int ic = 0; ic < 3; ++ic)
{
  J_c3 += coordinate_dofs[(ic) * 3 + 1] * FE1_C1_D01_Q083[0][0][0][ic];
}
double J_c0 = 0.0;
for (int ic = 0; ic < 3; ++ic)
{
  J_c0 += coordinate_dofs[(ic) * 3] * FE0_C0_D10_Q083[0][0][0][ic];
}
double J_c1 = 0.0;
for (int ic = 0; ic < 3; ++ic)
{
  J_c1 += coordinate_dofs[(ic) * 3] * FE1_C1_D01_Q083[0][0][0][ic];
}
double J_c2 = 0.0;
for (int ic = 0; ic < 3; ++ic)
{
  J_c2 += coordinate_dofs[(ic) * 3 + 1] * FE0_C0_D10_Q083[0][0][0][ic];
}
double w0_d01 = 0.0;
for (int ic = 0; ic < 3; ++ic)
{
  w0_d01 += w[ic] * FE1_C1_D01_Q083[0][0][0][ic];
}
double sp_0 = J_c0 * J_c3;
double sp_1 = J_c1 * J_c2;
double sp_2 = -sp_1;
double sp_3 = sp_0 + sp_2;
double sp_4 = J_c3 / sp_3;
double sp_5 = w0_d10 * sp_4;
double sp_6 = -J_c2;
double sp_7 = sp_6 / sp_3;
double sp_8 = w0_d01 * sp_7;
double sp_9 = sp_5 + sp_8;
double sp_10 = c[0] * sp_9;
for (int iq = 0; iq < 1; ++iq)
{
  A[iq + 0] += sp_10;
}

}