Skip to content

Commit ea3f52e

Browse files
schnellerhaseschnellerhase
committed
Add geometric twogrid demo
1 parent b21cc50 commit ea3f52e

File tree

4 files changed

+354
-0
lines changed

4 files changed

+354
-0
lines changed

‎cpp/demo/CMakeLists.txt

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -17,6 +17,7 @@ endmacro(add_demo_subdirectory)
1717

1818
# Add demos
1919
add_demo_subdirectory(custom_kernel)
20+
add_demo_subdirectory(geometric-multigrid)
2021
add_demo_subdirectory(poisson)
2122
add_demo_subdirectory(poisson_matrix_free)
2223
add_demo_subdirectory(hyperelasticity)

‎cpp/demo/geometric-twogrid/CMakeLists.txt

Lines changed: 67 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,67 @@
1+
# This file was generated by running
2+
#
3+
# python cmake/scripts/generate-cmakefiles from dolfinx/cpp
4+
#
5+
cmake_minimum_required(VERSION 3.19)
6+
7+
set(PROJECT_NAME demo_geometric-twogrid)
8+
project(${PROJECT_NAME} LANGUAGES C CXX)
9+
10+
# Set C++20 standard
11+
set(CMAKE_CXX_STANDARD 20)
12+
set(CMAKE_CXX_STANDARD_REQUIRED ON)
13+
set(CMAKE_CXX_EXTENSIONS OFF)
14+
15+
if(NOT TARGET dolfinx)
16+
find_package(DOLFINX REQUIRED)
17+
endif()
18+
19+
include(CheckSymbolExists)
20+
set(CMAKE_REQUIRED_INCLUDES ${PETSC_INCLUDE_DIRS})
21+
check_symbol_exists(PETSC_USE_COMPLEX petscsystypes.h PETSC_SCALAR_COMPLEX)
22+
check_symbol_exists(PETSC_USE_REAL_DOUBLE petscsystypes.h PETSC_REAL_DOUBLE)
23+
24+
# Add target to compile UFL files
25+
if(PETSC_SCALAR_COMPLEX EQUAL 1)
26+
if(PETSC_REAL_DOUBLE EQUAL 1)
27+
set(SCALAR_TYPE "--scalar_type=complex128")
28+
else()
29+
set(SCALAR_TYPE "--scalar_type=complex64")
30+
endif()
31+
else()
32+
if(PETSC_REAL_DOUBLE EQUAL 1)
33+
set(SCALAR_TYPE "--scalar_type=float64")
34+
else()
35+
set(SCALAR_TYPE "--scalar_type=float32")
36+
endif()
37+
endif()
38+
add_custom_command(
39+
OUTPUT poisson.c
40+
COMMAND ffcx ${CMAKE_CURRENT_SOURCE_DIR}/poisson.py ${SCALAR_TYPE}
41+
VERBATIM
42+
DEPENDS poisson.py
43+
COMMENT "Compile poisson.py using FFCx"
44+
)
45+
46+
set(CMAKE_INCLUDE_CURRENT_DIR ON)
47+
48+
add_executable(${PROJECT_NAME} main.cpp ${CMAKE_CURRENT_BINARY_DIR}/poisson.c)
49+
target_link_libraries(${PROJECT_NAME} dolfinx)
50+
51+
# Do not throw error for 'multi-line comments' (these are typical in rst which
52+
# includes LaTeX)
53+
include(CheckCXXCompilerFlag)
54+
check_cxx_compiler_flag("-Wno-comment" HAVE_NO_MULTLINE)
55+
set_source_files_properties(
56+
main.cpp
57+
PROPERTIES
58+
COMPILE_FLAGS
59+
"$<$<BOOL:${HAVE_NO_MULTLINE}>:-Wno-comment -Wall -Wextra -pedantic -Werror>"
60+
)
61+
62+
# Test targets (used by DOLFINx testing system)
63+
set(TEST_PARAMETERS2 -np 2 ${MPIEXEC_PARAMS} "./${PROJECT_NAME}")
64+
set(TEST_PARAMETERS3 -np 3 ${MPIEXEC_PARAMS} "./${PROJECT_NAME}")
65+
add_test(NAME ${PROJECT_NAME}_mpi_2 COMMAND "mpirun" ${TEST_PARAMETERS2})
66+
add_test(NAME ${PROJECT_NAME}_mpi_3 COMMAND "mpirun" ${TEST_PARAMETERS3})
67+
add_test(NAME ${PROJECT_NAME}_serial COMMAND ${PROJECT_NAME})

‎cpp/demo/geometric-twogrid/main.cpp

Lines changed: 242 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,242 @@
1+
// Copyright (C) 2024 Paul Kühner
2+
//
3+
// This file is part of DOLFINX (https://www.fenicsproject.org)
4+
//
5+
// SPDX-License-Identifier: LGPL-3.0-or-later
6+
7+
/// @cond TODO
8+
9+
#include <cassert>
10+
#include <cmath>
11+
#include <cstddef>
12+
#include <cstdint>
13+
#include <memory>
14+
#include <numbers>
15+
16+
#include <petscksp.h>
17+
#include <petsclog.h>
18+
#include <petscmat.h>
19+
#include <petscpc.h>
20+
#include <petscpctypes.h>
21+
#include <petscsys.h>
22+
#include <petscsystypes.h>
23+
#include <petscvec.h>
24+
#include <petscviewer.h>
25+
26+
#include <mpi.h>
27+
28+
#include <basix/finite-element.h>
29+
30+
#include <dolfinx/common/MPI.h>
31+
#include <dolfinx/fem/DirichletBC.h>
32+
#include <dolfinx/fem/FunctionSpace.h>
33+
#include <dolfinx/fem/petsc.h>
34+
#include <dolfinx/fem/utils.h>
35+
#include <dolfinx/io/ADIOS2Writers.h>
36+
#include <dolfinx/io/VTKFile.h>
37+
#include <dolfinx/io/XDMFFile.h>
38+
#include <dolfinx/la/SparsityPattern.h>
39+
#include <dolfinx/la/petsc.h>
40+
#include <dolfinx/mesh/cell_types.h>
41+
#include <dolfinx/mesh/generation.h>
42+
#include <dolfinx/mesh/utils.h>
43+
#include <dolfinx/refinement/refine.h>
44+
#include <dolfinx/transfer/transfer_matrix.h>
45+
46+
#include "poisson.h"
47+
48+
using namespace dolfinx;
49+
using T = PetscScalar;
50+
using U = typename dolfinx::scalar_value_type_t<T>;
51+
52+
struct PetscEnv
53+
{
54+
PetscEnv(int argc, char** argv) { PetscInitialize(&argc, &argv, NULL, NULL); }
55+
56+
~PetscEnv() { PetscFinalize(); }
57+
};
58+
59+
// recommended to run with
60+
// ./demo_geometric-multigrid -pc_mg_log -all_ksp_monitor -ksp_converged_reason
61+
// -ksp_view
62+
int main(int argc, char** argv)
63+
{
64+
PetscEnv petscEnv(argc, argv);
65+
// PetscLogDefaultBegin();
66+
67+
auto element = basix::create_element<U>(
68+
basix::element::family::P, basix::cell::type::tetrahedron, 1,
69+
basix::element::lagrange_variant::unset,
70+
basix::element::dpc_variant::unset, false);
71+
auto part = mesh::create_cell_partitioner(mesh::GhostMode::shared_vertex);
72+
73+
auto mesh_coarse = std::make_shared<mesh::Mesh<U>>(dolfinx::mesh::create_box(
74+
MPI_COMM_WORLD, {{{0, 0, 0}, {1, 1, 1}}}, {10, 10, 10},
75+
mesh::CellType::tetrahedron, part));
76+
77+
auto V_coarse = std::make_shared<fem::FunctionSpace<U>>(
78+
fem::create_functionspace<U>(mesh_coarse, element, {}));
79+
80+
// refinement routine requires edges to be initialized
81+
mesh_coarse->topology()->create_entities(1);
82+
auto [mesh, parent_cells, parent_facets]
83+
= dolfinx::refinement::refine(*mesh_coarse, std::nullopt);
84+
85+
auto V = std::make_shared<fem::FunctionSpace<U>>(fem::create_functionspace<U>(
86+
std::make_shared<mesh::Mesh<U>>(mesh), element, {}));
87+
88+
auto f_ana = [](auto x) -> std::pair<std::vector<T>, std::vector<std::size_t>>
89+
{
90+
std::vector<T> f;
91+
for (std::size_t p = 0; p < x.extent(1); ++p)
92+
{
93+
auto x0 = x(0, p);
94+
auto x1 = x(1, p);
95+
auto x2 = x(2, p);
96+
auto pi = std::numbers::pi;
97+
f.push_back(2 * pi * pi * std::sin(pi * x0) * std::sin(pi * x1)
98+
* std::sin(pi * x2));
99+
}
100+
return {f, {f.size()}};
101+
};
102+
auto f = std::make_shared<fem::Function<T>>(V);
103+
f->interpolate(f_ana);
104+
105+
{
106+
io::VTKFile file(MPI_COMM_WORLD, "f.pvd", "w");
107+
file.write<T>({*f}, 0.0);
108+
}
109+
110+
auto a = std::make_shared<fem::Form<T>>(
111+
fem::create_form<T>(*form_poisson_a, {V, V}, {}, {}, {}));
112+
auto a_coarse = std::make_shared<fem::Form<T>>(
113+
fem::create_form<T>(*form_poisson_a, {V_coarse, V_coarse}, {}, {}, {}));
114+
auto L = std::make_shared<fem::Form<T>>(
115+
fem::create_form<T>(*form_poisson_L, {V}, {{"f", f}}, {}, {}));
116+
la::petsc::Matrix A(fem::petsc::create_matrix(*a), true);
117+
la::petsc::Matrix A_coarse(fem::petsc::create_matrix(*a_coarse), true);
118+
119+
la::Vector<T> b(L->function_spaces()[0]->dofmap()->index_map,
120+
L->function_spaces()[0]->dofmap()->index_map_bs());
121+
122+
// TOOD: this somehow breaks?!?
123+
// V->mesh()->topology_mutable()->create_connectivity(2, 3);
124+
// mesh::exterior_facet_indices(*V->mesh()->topology())
125+
auto bc = std::make_shared<const fem::DirichletBC<T>>(
126+
0.0,
127+
mesh::locate_entities_boundary(
128+
*V->mesh(), 0,
129+
[](auto x) { return std::vector<std::int8_t>(x.extent(1), true); }),
130+
V);
131+
132+
// V_coarse->mesh()->topology_mutable()->create_connectivity(2, 3);
133+
auto bc_coarse = std::make_shared<const fem::DirichletBC<T>>(
134+
0.0,
135+
mesh::locate_entities_boundary(
136+
*V_coarse->mesh(), 0,
137+
[](auto x) { return std::vector<std::int8_t>(x.extent(1), true); }),
138+
V_coarse);
139+
140+
// assemble A
141+
MatZeroEntries(A.mat());
142+
fem::assemble_matrix(la::petsc::Matrix::set_block_fn(A.mat(), ADD_VALUES), *a,
143+
{bc});
144+
MatAssemblyBegin(A.mat(), MAT_FLUSH_ASSEMBLY);
145+
MatAssemblyEnd(A.mat(), MAT_FLUSH_ASSEMBLY);
146+
fem::set_diagonal<T>(la::petsc::Matrix::set_fn(A.mat(), INSERT_VALUES), *V,
147+
{bc});
148+
MatAssemblyBegin(A.mat(), MAT_FINAL_ASSEMBLY);
149+
MatAssemblyEnd(A.mat(), MAT_FINAL_ASSEMBLY);
150+
151+
// assemble b
152+
b.set(0.0);
153+
fem::assemble_vector(b.mutable_array(), *L);
154+
fem::apply_lifting<T, U>(b.mutable_array(), {a}, {{bc}}, {}, T(1));
155+
b.scatter_rev(std::plus<T>());
156+
fem::set_bc<T, U>(b.mutable_array(), {bc});
157+
158+
// assemble A_coarse
159+
MatZeroEntries(A_coarse.mat());
160+
fem::assemble_matrix(
161+
la::petsc::Matrix::set_block_fn(A_coarse.mat(), ADD_VALUES), *a_coarse,
162+
{bc_coarse});
163+
MatAssemblyBegin(A_coarse.mat(), MAT_FLUSH_ASSEMBLY);
164+
MatAssemblyEnd(A_coarse.mat(), MAT_FLUSH_ASSEMBLY);
165+
fem::set_diagonal<T>(la::petsc::Matrix::set_fn(A_coarse.mat(), INSERT_VALUES),
166+
*V_coarse, {bc_coarse});
167+
MatAssemblyBegin(A_coarse.mat(), MAT_FINAL_ASSEMBLY);
168+
MatAssemblyEnd(A_coarse.mat(), MAT_FINAL_ASSEMBLY);
169+
170+
KSP ksp;
171+
KSPCreate(MPI_COMM_WORLD, &ksp);
172+
KSPSetType(ksp, "cg");
173+
174+
PC pc;
175+
KSPGetPC(ksp, &pc);
176+
KSPSetFromOptions(ksp);
177+
PCSetType(pc, "mg");
178+
179+
PCMGSetLevels(pc, 2, NULL);
180+
PCMGSetType(pc, PC_MG_MULTIPLICATIVE);
181+
PCMGSetCycleType(pc, PC_MG_CYCLE_V);
182+
183+
// do not rely on coarse grid operators to be generated by
184+
// restriction/prolongation
185+
PCMGSetGalerkin(pc, PC_MG_GALERKIN_NONE);
186+
PCMGSetOperators(pc, 0, A_coarse.mat(), A_coarse.mat());
187+
188+
// PCMGSetNumberSmooth(pc, 1);
189+
PCSetFromOptions(pc);
190+
191+
mesh.topology_mutable()->create_connectivity(0, 1);
192+
mesh.topology_mutable()->create_connectivity(1, 0);
193+
auto inclusion_map = transfer::inclusion_mapping(*mesh_coarse, mesh);
194+
la::SparsityPattern sp
195+
= transfer::create_sparsity_pattern(*V_coarse, *V, inclusion_map);
196+
la::petsc::Matrix restriction(MPI_COMM_WORLD, sp);
197+
transfer::assemble_transfer_matrix<double>(
198+
la::petsc::Matrix::set_block_fn(restriction.mat(), INSERT_VALUES),
199+
*V_coarse, *V, inclusion_map,
200+
[](std::int32_t d) -> double { return d == 0 ? 1. : .5; });
201+
202+
MatAssemblyBegin(restriction.mat(), MAT_FINAL_ASSEMBLY);
203+
MatAssemblyEnd(restriction.mat(), MAT_FINAL_ASSEMBLY);
204+
205+
// PETSc figures out to use transpose by dimensions!
206+
PCMGSetInterpolation(pc, 1, restriction.mat());
207+
PCMGSetRestriction(pc, 1, restriction.mat());
208+
209+
// MatView(A.mat(), PETSC_VIEWER_STDOUT_SELF);
210+
KSPSetOperators(ksp, A.mat(), A.mat());
211+
KSPSetUp(ksp);
212+
213+
auto u = std::make_shared<fem::Function<T>>(V);
214+
215+
la::petsc::Vector _u(la::petsc::create_vector_wrap(*u->x()), false);
216+
la::petsc::Vector _b(la::petsc::create_vector_wrap(b), false);
217+
218+
KSPSolve(ksp, _b.vec(), _u.vec());
219+
u->x()->scatter_fwd();
220+
221+
{
222+
io::VTKFile file(MPI_COMM_WORLD, "u.pvd", "w");
223+
file.write<T>({*u}, 0.0);
224+
225+
io::XDMFFile file_xdmf(mesh.comm(), "u_xdmf.xdmf", "w");
226+
file_xdmf.write_mesh(mesh);
227+
file_xdmf.write_function(*u, 0.0);
228+
file_xdmf.close();
229+
230+
#ifdef HAS_ADIOS2
231+
io::VTXWriter<U> vtx_writer(mesh.comm(), std::filesystem::path("u_vtx.bp"),
232+
{u}, io::VTXMeshPolicy::reuse);
233+
vtx_writer.write(0);
234+
#endif
235+
}
236+
237+
KSPDestroy(&ksp);
238+
239+
// PetscLogView(PETSC_VIEWER_STDOUT_SELF);
240+
}
241+
242+
/// @endcond

‎cpp/demo/geometric-twogrid/poisson.py

Lines changed: 44 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,44 @@
1+
# The first step is to define the variational problem at hand. We define
2+
# the variational problem in UFL terms in a separate form file
3+
# {download}`demo_poisson/poisson.py`. We begin by defining the finite
4+
# element:
5+
6+
from basix.ufl import element
7+
from ufl import (
8+
Coefficient,
9+
FunctionSpace,
10+
Mesh,
11+
TestFunction,
12+
TrialFunction,
13+
dx,
14+
grad,
15+
inner,
16+
)
17+
18+
e = element("Lagrange", "tetrahedron", 1)
19+
20+
# The first argument to :py:class:`FiniteElement` is the finite element
21+
# family, the second argument specifies the domain, while the third
22+
# argument specifies the polynomial degree. Thus, in this case, our
23+
# element `element` consists of first-order, continuous Lagrange basis
24+
# functions on triangles (or in order words, continuous piecewise linear
25+
# polynomials on triangles).
26+
#
27+
# Next, we use this element to initialize the trial and test functions
28+
# ($u$ and $v$) and the coefficient functions ($f$ and $g$):
29+
30+
coord_element = element("Lagrange", "tetrahedron", 1, shape=(3,))
31+
mesh = Mesh(coord_element)
32+
33+
V = FunctionSpace(mesh, e)
34+
35+
u = TrialFunction(V)
36+
v = TestFunction(V)
37+
f = Coefficient(V)
38+
g = Coefficient(V)
39+
40+
# Finally, we define the bilinear and linear forms according to the
41+
# variational formulation of the equations:
42+
43+
a = inner(grad(u), grad(v)) * dx
44+
L = inner(f, v) * dx

0 commit comments

Comments
 (0)