Dear all,
I have a time dependent problem which I need to solve multiple times within a larger iteration loop, so I need the solution to be as fast as possible. I therefore turned to the c++ interface. As a starter, i implemented a simple heat equation with a gaussian initial condition, but it seems that the new solution function is not being updated into the previous solution at each time step. Here is the c++ code:
//#include <dolfinx.h>
#include "heat_equation.h" // From FFCX
#include <basix/finite-element.h>
#include <cmath>
#include <dolfinx.h>
#include <dolfinx/fem/Constant.h>
#include <dolfinx/fem/petsc.h>
#include <dolfinx/la/petsc.h>
#include <petscmat.h>
#include <petscsys.h>
#include <petscsystypes.h>
#include <utility>
#include <dolfinx/io/XDMFFile.h>
#include <vector>
using namespace dolfinx;
using T = PetscScalar;
using U = typename dolfinx::scalar_value_type<T>;
int main(int argc, char* argv[])
{ // Step 1: Create
dolfinx::init_logging(argc, argv);
PetscInitialize(&argc, &argv, nullptr, nullptr);
// Create mesh and function space
auto part = mesh::create_cell_partitioner(mesh::GhostMode::shared_facet);
auto mesh = std::make_shared<mesh::Mesh<double>>(
mesh::create_rectangle(MPI_COMM_WORLD, {{{0.0, 0.0}, {2.0, 2.0}}},
{32, 16}, mesh::CellType::triangle, part));
// Create finite element
auto element = basix::create_element<double>(
basix::element::family::P, basix::cell::type::triangle, 1,
basix::element::lagrange_variant::unset,
basix::element::dpc_variant::unset, false); // Continuity (continuous Galerkin)
// Create function space
auto V = std::make_shared<fem::FunctionSpace<double>>(
fem::create_functionspace(mesh, element));
// Create initial condition function (e.g., Gaussian bump)
//auto u0 = fem::Function<double>(V);
auto u = std::make_shared<fem::Function<T>>(V);
auto u_n = std::make_shared<fem::Function<T>>(V); // Previous solution
u_n->interpolate(
[](auto x) -> std::pair<std::vector<T>, std::vector<std::size_t>>
{
std::vector<T> f;
for (std::size_t p = 0; p < x.extent(1); ++p)
{
auto dx = (x(0, p) - 1.0) * (x(0, p) - 1.0);
auto dy = (x(1, p) - 1.0) * (x(1, p) - 1.0);
f.push_back(10 * std::exp(-(dx + dy) / 0.02));
}
return {f, {f.size()}};
});
// Interpolate the initial condition function into the solution function u_n
//u_n->interpolate(initial_condition);
// auto u0_const = std::make_shared<fem::Constant<T>>(1.0);
// Interpolate initial condition into u_n
//u_n->interpolate([](auto x) -> T { return /* initial condition expression */ 1.0; });
// Define alpha as a fem::Constant
auto alpha = std::make_shared<fem::Constant<T>>(20.0);
double deltaT = 0.01;
auto dt = std::make_shared<fem::Constant<T>>(deltaT);
std::vector<std::shared_ptr<const fem::Function<T>>> coeffs = {u_n};
auto a = fem::create_form<T>(*form_heat_equation_a, {V, V}, {},{{"dt",dt}}, {}, {});
auto L_form = fem::create_form<T>(*form_heat_equation_L, {V},{{"u_n",u_n}}, {{"dt",dt}}, {}, {});
//u_n->x()->set(1.0);
// Assemble matrix once if it’s time-invariant
// Correcting matrix assembly
la::petsc::Matrix A(fem::petsc::create_matrix(a), false); // Create the PETSc matrix
MatZeroEntries(A.mat());
fem::assemble_matrix(la::petsc::Matrix::set_block_fn(A.mat(), ADD_VALUES),
a, {});
MatAssemblyBegin(A.mat(), MAT_FLUSH_ASSEMBLY);
MatAssemblyEnd(A.mat(), MAT_FLUSH_ASSEMBLY);
MatAssemblyBegin(A.mat(), MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(A.mat(), MAT_FINAL_ASSEMBLY);
//fem::assemble_matrix(A, a);
//A.assemble();
// Pass the form directly without dereferencing
//fem::assemble_matrix(A, a); // a_form is already of type Form<T>, no need for *a
int num_steps = int(1.0/deltaT);
// Create linear solver
//auto solver = std::make_unique<la::PETScKrylovSolver>(MPI_COMM_WORLD);
//solver->set_operator(A);
// Create PETSc linear solver
//KSP ksp;
//KSPCreate(V->mesh()->comm(), &ksp);
//KSPSetOperators(ksp, A.mat(), A.mat());
la::petsc::KrylovSolver lu(MPI_COMM_WORLD);
la::petsc::options::set("ksp_type", "cg");
la::petsc::options::set("pc_type", "lu");
lu.set_from_options();
lu.set_operator(A.mat());
// Create vector for RHS
la::Vector<T> b(L_form.function_spaces()[0]->dofmap()->index_map,
L_form.function_spaces()[0]->dofmap()->index_map_bs());
for (int n = 0; n < num_steps; ++n)
{
double t = (n + 1) * deltaT;
std::cout << "Time step " << n + 1 << ", Time = " << t << std::endl;
// Assemble right-hand side vector using u_n
b.set(0.0);
fem::assemble_vector(b.mutable_array(), L_form);
b.scatter_rev(std::plus<T>());
//bc.set(b.mutable_array(), std::nullopt)
//fem::apply_lifting(b, {*a_form}, {{u}}, {}, 1.0);
//b.ghostUpdate(la::GhostMode::add);
//fem::set_bc(b, {}); // Pass Dirichlet BCs if any
// Solve linear system A * u = b
// Solve
//KSPSolve(ksp, b.vec(), u->x()->vec());
la::petsc::Vector _u(la::petsc::create_vector_wrap(*u->x()), false);
la::petsc::Vector _b(la::petsc::create_vector_wrap(b), false);
lu.solve(_u.vec(), _b.vec());
u->x()->scatter_fwd();
std::cout << "u[176] = " << (u_n->x()->array())[176] << std::endl;
// Update previous solution
// u_n->set(*u);
//u_n->x()->copy_from(*u->x());
//VecCopy(la::petsc::create_vector_wrap(*u_n->x()), _u.vec());
std::ranges::copy(u->x()->array(), u_n->x()->mutable_array().begin());
u_n->x()->scatter_fwd();
// Update previous solution
//std::swap(u, u_n);
//for (std::size_t i = 0; i < (u_n->x()->array()).size(); ++i)
//{
//std::cout << "u[" << i << "] = " << (u_n->x()->array())[i] << std::endl;
//}
}
io::VTKFile file_u(mesh->comm(), "u.pvd", "w");
file_u.write<T>({*u}, 0.0);
//io::XDMFWriter file(MPI_COMM_WORLD, "u.xdmf", "w");
//file.write<T>({*u}, 0.0)
// io::XDMFFile file_sigma(mesh->comm(), "u.xdmf", "w");
// file_sigma.write_mesh(*mesh);
// file_sigma.write_function(u, 0.0);
// std::shared_ptr<io::XDMFFile> xdmf_file = std::make_shared<io::XDMFFile>(MPI_COMM_WORLD, "solution.xdmf", io::FileMode::WRITE);
// xdmf_file->write_mesh(*mesh);
// xdmf_file->write_function(*u, "solution");
return 0;
}
As well as the ufl code that defines the relevant forms:
from basix.ufl import element
from ufl import (
Coefficient,
Constant,
FunctionSpace,
Mesh,
TestFunction,
TrialFunction,
ds,
dx,
grad,
inner,
dot,
)
from dolfinx import fem, cpp
e = element("Lagrange", "triangle", 1)
coord_element = element("Lagrange", "triangle", 1, shape=(2,))
mesh = Mesh(coord_element)
V = FunctionSpace(mesh, e)
u = TrialFunction(V)
v = TestFunction(V)
u_n = Coefficient(V)
alpha = Constant(mesh)
dt = Constant(mesh)
a = u*v*dx + dt*alpha*dot(grad(u), grad(v))*dx/2.0
L = u_n*v*dx - dt*alpha*dot(grad(u_n), grad(v))*dx/2.0
form_a = fem.form(a)
form_L = fem.form(L)
and the CmakeLists:
cmake_minimum_required(VERSION 3.16)
project(DolfinxTimeDependent)
set(CMAKE_CXX_STANDARD 17)
if(NOT TARGET dolfinx)
find_package(DOLFINX REQUIRED)
endif()
#find_package(ufcx REQUIRED HINTS $ENV{CONDA_PREFIX}/lib/cmake/ufcx)
include_directories(SYSTEM /dolfinx-env/lib/python3.12/site-packages/ffcx/codegeneration)
#include_directories(/usr/local/include)
#find_package(DOLFINX REQUIRED)
#find_package(DOLFINX REQUIRED HINTS
# $ENV{CONDA_PREFIX}/lib/cmake/dolfinx
# /usr/lib/cmake/dolfinx
# /usr/local/lib/cmake/dolfinx
#)
# Add C source from ffcx
add_library(forms STATIC
${CMAKE_CURRENT_SOURCE_DIR}/heat_equation.c
)
target_include_directories(forms PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})
add_executable(time_sim main.cpp)
target_link_libraries(time_sim PRIVATE forms dolfinx)
I am quite new to the c++ interface of dolfinx so any help would be deeply appreciated.