Time dependent problem in c++ not updating

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.

So your question is, why doesn’t this work?
Have you tried printing the entries of u->x()->array() pre and post copy (same with u_n->x()->array())?

(As Im currently not at my computer, I can’t run the code)

I just realized there was a typo in the code which improperly defined the weak form, after fixing now the code propagates the solution just fine. I just wasnt sure if I was doing something wrong because its my first time using the c++ interface. Thanks for taking a look though!

1 Like