Mismatch between number of nodes and number of solution data

[Sujit_Kumar_Nath]

Hi there,

I am trying to reproduce the Poisson equation with periodic BC given in the example in

jorgensd/dolfinx_mpc/blob/main/python/demos/demo_periodic_geometrical.py#L51

However, when I am trying to extract the node points in the mesh and the solution value at every node point, there is a mismatch between the number of node points and the number of solution data. In fact, number of solution data is exactly double than the number of node points. Could you please help me to understand why this is happening and how I can extract the solutions at the node points? Below I copy my code and its out put.

####################################

My code is:

####################################

This demo program solves Poisson’s equation

- div grad u(x, y) = f(x, y)

on the unit square with homogeneous Dirichlet boundary conditions

at y = 0, 1 and periodic boundary conditions at x = 0, 1.

Copyright (C) Jørgen S. Dokken 2020-2022.

This file is part of DOLFINX_MPCX.

SPDX-License-Identifier: MIT

from future import annotations

from pathlib import Path
from typing import Union

from mpi4py import MPI
from petsc4py import PETSc

import dolfinx.fem as fem
import numpy as np
import scipy.sparse.linalg
from dolfinx import default_scalar_type
from dolfinx.common import Timer, TimingType, list_timings
from dolfinx.io import XDMFFile
from dolfinx.mesh import create_unit_square, locate_entities_boundary
from ufl import (SpatialCoordinate, TestFunction, TrialFunction, as_vector, dx,
exp, grad, inner, pi, sin)

import dolfinx_mpc.utils
from dolfinx_mpc import LinearProblem, MultiPointConstraint

import pandas as pd

Get PETSc int and scalar types

complex_mode = True if np.dtype(default_scalar_type).kind == ‘c’ else False

Create mesh and finite element

NX = 20
NY = 20
mesh = create_unit_square(MPI.COMM_WORLD, NX, NY)
V = fem.functionspace(mesh, (“Lagrange”, 1, (mesh.geometry.dim, )))
tol = 250 * np.finfo(default_scalar_type).resolution

Listing space points ### *********************************************************

dof_coordinates = V.tabulate_dof_coordinates()
df = pd.DataFrame(dof_coordinates)
df.to_csv(‘outdata/xyzs.csv’, index = None, header=False)

def dirichletboundary(x):
return np.logical_or(np.isclose(x[1], 0, atol=tol), np.isclose(x[1], 1, atol=tol))

Create Dirichlet boundary condition

facets = locate_entities_boundary(mesh, 1, dirichletboundary)
topological_dofs = fem.locate_dofs_topological(V, 1, facets)
zero = np.array([0, 0], dtype=default_scalar_type)
bc = fem.dirichletbc(zero, topological_dofs, V)
bcs = [bc]

def periodic_boundary(x):
return np.isclose(x[0], 1, atol=tol)

def periodic_relation(x):
out_x = np.zeros_like(x)
out_x[0] = 1 - x[0]
out_x[1] = x[1]
out_x[2] = x[2]
return out_x

with Timer(“~PERIODIC: Initialize MPC”):
mpc = MultiPointConstraint(V)
mpc.create_periodic_constraint_geometrical(V, periodic_boundary, periodic_relation, bcs)
mpc.finalize()

Define variational problem

u = TrialFunction(V)
v = TestFunction(V)
a = inner(grad(u), grad(v)) * dx

x = SpatialCoordinate(mesh)
dx_ = x[0] - 0.9
dy_ = x[1] - 0.5
f = as_vector((x[0] * sin(5.0 * pi * x[1]) + 1.0 * exp(-(dx_ * dx_ + dy_ * dy_) / 0.02), 0.3 * x[1]))

rhs = inner(f, v) * dx

Setup MPC system

with Timer(“~PERIODIC: Initialize varitional problem”):
problem = LinearProblem(a, rhs, mpc, bcs=bcs)

solver = problem.solver

Give PETSc solver options a unique prefix

solver_prefix = “dolfinx_mpc_solve_{}”.format(id(solver))
solver.setOptionsPrefix(solver_prefix)

petsc_options: dict[str, Union[str, int, float]]
if complex_mode or default_scalar_type == np.float32:
petsc_options = {“ksp_type”: “preonly”, “pc_type”: “lu”}
else:
petsc_options = {“ksp_type”: “cg”, “ksp_rtol”: 1e-6, “pc_type”: “hypre”, “pc_hypre_type”: “boomeramg”,
“pc_hypre_boomeramg_max_iter”: 1, “pc_hypre_boomeramg_cycle_type”: “v” # ,
# “pc_hypre_boomeramg_print_statistics”: 1
}

Set PETSc options

opts = PETSc.Options() # type: ignore
opts.prefixPush(solver_prefix)
if petsc_options is not None:
for k, v in petsc_options.items():
opts[k] = v
opts.prefixPop()
solver.setFromOptions()

with Timer(“~PERIODIC: Assemble and solve MPC problem”):
uh = problem.solve()
# solver.view()
it = solver.getIterationNumber()
print(“Constrained solver iterations {0:d}”.format(it))

Write solution to file

outdir = Path(“results”)
outdir.mkdir(exist_ok=True, parents=True)

uh.name = “u_mpc”
outfile = XDMFFile(mesh.comm, outdir / “demo_periodic_geometrical.xdmf”, “w”)
outfile.write_mesh(mesh)
outfile.write_function(uh)

Writing the solution to a file

solData = uh.x.array
#finalSolData = np.array([dof_coordinates[i].tolist() + [solData[i]] ] for i in range(len(solData)))
print(“Number of degrees of freedom: “, len(dof_coordinates))
print(dof_coordinates.shape)
print(”=======================================”)
print("Number of solution data points: ",len(solData))
print(solData.shape)
df = pd.DataFrame(solData)
df.to_csv(‘outdata/exampsol.csv’, index = None, header=False)

############################################################

The output is:

############################################################
Constrained solver iterations 4
Number of degrees of freedom: 441
(441, 3)

Number of solution data points: 882
(882,)

Thanks,
Sujit

[dokken]

When you use vector function spaces, you have several degrees of freedom at the same physical point. In your case, you have that the 0-th coordinate from tabulate dof coordinates corresponds to the first two entries in u.x.array, where the first entry is the x component, the second entry is the y component.

[Sujit_Kumar_Nath]

Thanks a lot! Much appreciated.

So, do I need to change

V = fem.functionspace(mesh, (“Lagrange”, 1, (mesh.geometry.dim, )))

and

f = as_vector((x[0] * sin(5.0 * pi * x[1]) + 1.0 * exp(-(dx_ * dx_ + dy_ * dy_) / 0.02), 0.3 * x[1]))

lines in my code?

Also, could you please let me know where I can see what “as_vector()” and “uh.x” return? I tried to find in Dolphinx 0.7.0 documentation but unfortunately couldn’t find the description of these functions in the documentation.

[dokken]

@Sujit_Kumar_Nath

So, do I need to change

V = fem.functionspace(mesh, (“Lagrange”, 1, (mesh.geometry.dim, )))

and

f = as_vector((x[0] * sin(5.0 * pi * x[1]) + 1.0 * exp(-(dx_ * dx_ + dy_ * dy_) / 0.02), 0.3 * x[1]))

lines in my code?

Why would you need to change it? Do you want a periodic Poisson problem of a vector field, or do you want a periodic problem for a scalar quantity?

@Sujit_Kumar_Nath

Also, could you please let me know where I can see what “as_vector()” and “uh.x” return?

as_vector is a ufl command, as indicated by its import Form language — Unified Form Language (UFL) 2021.1.0 documentation

Uh.x: dolfinx.fem — DOLFINx 0.7.2 documentation
and
dolfinx.la — DOLFINx 0.7.2 documentation

[Sujit_Kumar_Nath]

I needed periodic Poisson problem for scalar quantity. I have changed my code accordingly (pasted below). This is working for me. And thank you very much for sending the documentation links.

##########################

This demo program solves Poisson’s equation

- div grad u(x, y) = f(x, y)

on the unit square with homogeneous Dirichlet boundary conditions

at y = 0, 1 and periodic boundary conditions at x = 0, 1.

Copyright (C) Jørgen S. Dokken 2020-2022.

This file is part of DOLFINX_MPCX.

SPDX-License-Identifier: MIT

from future import annotations

from pathlib import Path
from typing import Union

from mpi4py import MPI
from petsc4py import PETSc

import dolfinx.fem as fem
import numpy as np
import scipy.sparse.linalg
from dolfinx import default_scalar_type
from dolfinx.common import Timer, TimingType, list_timings
from dolfinx.io import XDMFFile
from dolfinx.mesh import create_unit_square, locate_entities_boundary
from ufl import (SpatialCoordinate, TestFunction, TrialFunction, as_vector, dx,
exp, grad, inner, pi, sin)

import dolfinx_mpc.utils
from dolfinx_mpc import LinearProblem, MultiPointConstraint

import pandas as pd

Get PETSc int and scalar types

complex_mode = True if np.dtype(default_scalar_type).kind == ‘c’ else False

Create mesh and finite element

NX = 20
NY = 20
mesh = create_unit_square(MPI.COMM_WORLD, NX, NY)
#V = fem.functionspace(mesh, (“Lagrange”, 1, (mesh.geometry.dim, )))
V = fem.functionspace(mesh, (“CG”, 1))
tol = 250 * np.finfo(default_scalar_type).resolution

Listing space points ### *********************************************************

dof_coordinates = V.tabulate_dof_coordinates()
df = pd.DataFrame(dof_coordinates)
df.to_csv(‘outdata/xyzs.csv’, index = None, header=False)

def dirichletboundary(x):
return np.logical_or(np.isclose(x[1], 0, atol=tol), np.isclose(x[1], 1, atol=tol))

# Create Dirichlet boundary condition

facets = locate_entities_boundary(mesh, 1, dirichletboundary)
topological_dofs = fem.locate_dofs_topological(V, 1, facets)
zero = np.array(0, dtype=default_scalar_type)
bc = fem.dirichletbc(zero, topological_dofs, V)
bcs = [bc]

def periodic_boundary(x):
return np.isclose(x[0], 1, atol=tol)

def periodic_relation(x):
out_x = np.zeros_like(x)
out_x[0] = 1 - x[0]
out_x[1] = x[1]
out_x[2] = x[2]
return out_x

with Timer(“~PERIODIC: Initialize MPC”):
mpc = MultiPointConstraint(V)
mpc.create_periodic_constraint_geometrical(V, periodic_boundary, periodic_relation, bcs)
mpc.finalize()

Define variational problem

u = TrialFunction(V)
v = TestFunction(V)
a = inner(grad(u), grad(v)) * dx

x = SpatialCoordinate(mesh)
dx_ = x[0] - 0.9
dy_ = x[1] - 0.5

f = exp(-12 * ((x[0]-0.5)**2 + (x[1] - 0.5)**2))

rhs = f * v * dx

Setup MPC system

with Timer(“~PERIODIC: Initialize varitional problem”):
problem = LinearProblem(a, rhs, mpc, bcs=bcs)

solver = problem.solver

Give PETSc solver options a unique prefix

solver_prefix = “dolfinx_mpc_solve_{}”.format(id(solver))
solver.setOptionsPrefix(solver_prefix)

petsc_options: dict[str, Union[str, int, float]]
if complex_mode or default_scalar_type == np.float32:
petsc_options = {“ksp_type”: “preonly”, “pc_type”: “lu”}
else:
petsc_options = {“ksp_type”: “cg”, “ksp_rtol”: 1e-6, “pc_type”: “hypre”, “pc_hypre_type”: “boomeramg”,
“pc_hypre_boomeramg_max_iter”: 1, “pc_hypre_boomeramg_cycle_type”: “v” # ,
# “pc_hypre_boomeramg_print_statistics”: 1
}

Set PETSc options

opts = PETSc.Options() # type: ignore
opts.prefixPush(solver_prefix)
if petsc_options is not None:
for k, v in petsc_options.items():
opts[k] = v
opts.prefixPop()
solver.setFromOptions()

with Timer(“~PERIODIC: Assemble and solve MPC problem”):
uh = problem.solve()
# solver.view()
it = solver.getIterationNumber()
print(“Constrained solver iterations {0:d}”.format(it))

Write solution to file

outdir = Path(“results”)
outdir.mkdir(exist_ok=True, parents=True)

uh.name = “u_mpc”
outfile = XDMFFile(mesh.comm, outdir / “demo_periodic_geometrical.xdmf”, “w”)
outfile.write_mesh(mesh)
outfile.write_function(uh)

Writing the solution to a file

solData = uh.x.array
finalSolData = np.array([dof_coordinates[i].tolist() + [solData[i]] for i in range(len(solData))])
#print(np.array([dof_coordinates[i].tolist() + [solData[i]] for i in range(0,3)]))
print(“Number of degrees of freedom: “, len(dof_coordinates))
print(dof_coordinates.shape)
print(”=======================================”)
print("Number of solution data points: ",len(solData))
print(solData.shape)
df = pd.DataFrame(finalSolData)
df.to_csv(‘outdata/exampsol.csv’, index = None, header=False)

Thanks a lot!
Sujit

[dokken]

Please format your code with 3x` encapsulation, i.e.

```python

# Add code here

```

[Sujit_Kumar_Nath]

Below I copy my code.

@dokken

Please format your code with 3x` encapsulation, i.e.

```python

# Add code here

```
# This demo program solves Poisson's equation
#
#     - div grad u(x, y) = f(x, y)
#
# on the unit square with homogeneous Dirichlet boundary conditions
# at y = 0, 1 and periodic boundary conditions at x = 0, 1.
#
# Copyright (C) Jørgen S. Dokken 2020-2022.
#
# This file is part of DOLFINX_MPCX.
#
# SPDX-License-Identifier:    MIT
from __future__ import annotations

from pathlib import Path
from typing import Union

from mpi4py import MPI
from petsc4py import PETSc

import dolfinx.fem as fem
import numpy as np
import scipy.sparse.linalg
from dolfinx import default_scalar_type
from dolfinx.common import Timer, TimingType, list_timings
from dolfinx.io import XDMFFile
from dolfinx.mesh import create_unit_square, locate_entities_boundary
from ufl import (SpatialCoordinate, TestFunction, TrialFunction, as_vector, dx,
                 exp, grad, inner, pi, sin)

import dolfinx_mpc.utils
from dolfinx_mpc import LinearProblem, MultiPointConstraint

import pandas as pd 

# Get PETSc int and scalar types
complex_mode = True if np.dtype(default_scalar_type).kind == 'c' else False

# Create mesh and finite element
NX = 20
NY = 20
mesh = create_unit_square(MPI.COMM_WORLD, NX, NY)
#V = fem.functionspace(mesh, ("Lagrange", 1, (mesh.geometry.dim, )))
V = fem.functionspace(mesh, ("CG", 1))
tol = 250 * np.finfo(default_scalar_type).resolution

# Listing space points ### ********************************************************* 
dof_coordinates = V.tabulate_dof_coordinates() 
df = pd.DataFrame(dof_coordinates)
df.to_csv('outdata/xyzs.csv', index = None, header=False) 


def dirichletboundary(x):
    return np.logical_or(np.isclose(x[1], 0, atol=tol), np.isclose(x[1], 1, atol=tol))


# # Create Dirichlet boundary condition
facets = locate_entities_boundary(mesh, 1, dirichletboundary)
topological_dofs = fem.locate_dofs_topological(V, 1, facets)
zero = np.array(0, dtype=default_scalar_type)
bc = fem.dirichletbc(zero, topological_dofs, V)
bcs = [bc]

def periodic_boundary(x):
    return np.isclose(x[0], 1, atol=tol)


def periodic_relation(x):
    out_x = np.zeros_like(x)
    out_x[0] = 1 - x[0]
    out_x[1] = x[1]
    out_x[2] = x[2]
    return out_x


with Timer("~PERIODIC: Initialize MPC"):
    mpc = MultiPointConstraint(V)
    mpc.create_periodic_constraint_geometrical(V, periodic_boundary, periodic_relation, bcs)
    mpc.finalize()

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
a = inner(grad(u), grad(v)) * dx

x = SpatialCoordinate(mesh)
dx_ = x[0] - 0.9
dy_ = x[1] - 0.5

f = exp(-12 * ((x[0]-0.5)**2 + (x[1] - 0.5)**2))

rhs = f * v * dx


# Setup MPC system
with Timer("~PERIODIC: Initialize varitional problem"):
    problem = LinearProblem(a, rhs, mpc, bcs=bcs)

solver = problem.solver

# Give PETSc solver options a unique prefix
solver_prefix = "dolfinx_mpc_solve_{}".format(id(solver))
solver.setOptionsPrefix(solver_prefix)

petsc_options: dict[str, Union[str, int, float]]
if complex_mode or default_scalar_type == np.float32:
    petsc_options = {"ksp_type": "preonly", "pc_type": "lu"}
else:
    petsc_options = {"ksp_type": "cg", "ksp_rtol": 1e-6, "pc_type": "hypre", "pc_hypre_type": "boomeramg",
                     "pc_hypre_boomeramg_max_iter": 1, "pc_hypre_boomeramg_cycle_type": "v"  # ,
                     # "pc_hypre_boomeramg_print_statistics": 1
                     }

# Set PETSc options
opts = PETSc.Options()  # type: ignore
opts.prefixPush(solver_prefix)
if petsc_options is not None:
    for k, v in petsc_options.items():
        opts[k] = v
opts.prefixPop()
solver.setFromOptions()


with Timer("~PERIODIC: Assemble and solve MPC problem"):
    uh = problem.solve()
    # solver.view()
    it = solver.getIterationNumber()
    print("Constrained solver iterations {0:d}".format(it))

# Write solution to file
outdir = Path("results")
outdir.mkdir(exist_ok=True, parents=True)

uh.name = "u_mpc"
outfile = XDMFFile(mesh.comm, outdir / "demo_periodic_geometrical.xdmf", "w")
outfile.write_mesh(mesh)
outfile.write_function(uh)
    
    
### Writing the solution to a file 
solData = uh.x.array
finalSolData = np.array([dof_coordinates[i].tolist() + [solData[i]] for i in range(len(solData))])
#print(np.array([dof_coordinates[i].tolist() + [solData[i]] for i in range(0,3)]))
print("Number of degrees of freedom: ", len(dof_coordinates))
print(dof_coordinates.shape)
print("=======================================")
print("Number of solution data points: ",len(solData))
print(solData.shape)
df = pd.DataFrame(finalSolData)
df.to_csv('outdata/exampsol.csv', index = None, header=False)