Skip to content

Commit

Permalink
Update fenicsx codes for latest DOLFINx
Browse files Browse the repository at this point in the history
  • Loading branch information
@adeebkor
adeebkor committed Sep 10, 2024
1 parent e049c23 commit 41e8362
Show file tree
Hide file tree
Showing 17 changed files with 231 additions and 105 deletions.
50 changes: 32 additions & 18 deletions cpp/fenicsx/common/Linear.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,7 +51,8 @@ void axpy(la::Vector<T>& r, T alpha, const la::Vector<T>& x, const la::Vector<T>
template <typename T, int P>
class LinearSpectral {
public:
LinearSpectral(std::shared_ptr<mesh::Mesh<T>> Mesh,
LinearSpectral(basix::FiniteElement<T> element,
std::shared_ptr<mesh::Mesh<T>> Mesh,
std::shared_ptr<mesh::MeshTags<std::int32_t>> FacetTags,
std::shared_ptr<fem::Function<T>> speedOfSound,
std::shared_ptr<fem::Function<T>> density, const T& sourceFrequency,
Expand All @@ -76,7 +77,7 @@ class LinearSpectral {

// Define function space
V = std::make_shared<fem::FunctionSpace<T>>(
fem::create_functionspace(functionspace_form_forms_a, "u", mesh));
fem::create_functionspace(mesh, element));

// Define field functions
index_map = V->dofmap()->index_map;
Expand All @@ -95,20 +96,33 @@ class LinearSpectral {
std::fill(u_.begin(), u_.end(), 1.0);

// Compute exterior facets
std::map<fem::IntegralType,
std::vector<std::pair<std::int32_t, std::span<const std::int32_t>>>> fd;
auto facet_domains = fem::compute_integration_domains(
fem::IntegralType::exterior_facet, *V->mesh()->topology_mutable(),
ft->indices(), mesh->topology()->dim() - 1, ft->values());
for (auto& facet : facet_domains) {
std::vector<std::pair<std::int32_t, std::span<const std::int32_t>>> x;
x.emplace_back(facet.first, std::span(facet.second.data(), facet.second.size()));
fd.insert({fem::IntegralType::exterior_facet, std::move(x)});
}
std::vector<std::int32_t> ft_unique(ft->values().size());
std::copy(ft->values().begin(), ft->values().end(), ft_unique.begin());
std::sort(ft_unique.begin(), ft_unique.end());
auto it = std::unique(ft_unique.begin(), ft_unique.end());
ft_unique.erase(it, ft_unique.end());

std::map<fem::IntegralType, std::vector<std::pair<std::int32_t, std::vector<std::int32_t>>>> fd;
std::map<fem::IntegralType, std::vector<std::pair<std::int32_t, std::span<const std::int32_t>>>> fd_view;

std::vector<std::int32_t> facet_domains;
for (auto& tag : ft_unique) {
facet_domains = fem::compute_integration_domains(
fem::IntegralType::exterior_facet, *V->mesh()->topology_mutable(),
ft->find(tag), mesh->topology()->dim()-1);
fd[fem::IntegralType::exterior_facet].push_back(
{tag, facet_domains});
}

for (auto const& [key, val] : fd) {
for (auto const& [tag, vec] : val) {
fd_view[key].push_back({tag, std::span(vec.data(), vec.size())});
}
}

// Define LHS form
a = std::make_shared<fem::Form<T, T>>(
fem::create_form<T, T>(*form_forms_a, {V}, {{"u", u}, {"c0", c0}, {"rho0", rho0}}, {}, {}));
a = std::make_shared<fem::Form<T>>(
fem::create_form<T>(*form_forms_a, {V}, {{"u", u}, {"c0", c0}, {"rho0", rho0}}, {}, {}, {}));

m = std::make_shared<la::Vector<T>>(index_map, bs);
m_ = m->mutable_array();
Expand All @@ -117,9 +131,9 @@ class LinearSpectral {
m->scatter_rev(std::plus<T>());

// Define RHS form
L = std::make_shared<fem::Form<T, T>>(fem::create_form<T, T>(
L = std::make_shared<fem::Form<T>>(fem::create_form<T>(
*form_forms_L, {V}, {{"g", g}, {"u_n", u_n}, {"v_n", v_n}, {"c0", c0}, {"rho0", rho0}}, {},
fd));
fd_view, {}, {}));

b = std::make_shared<la::Vector<T>>(index_map, bs);
b_ = b->mutable_array();
Expand Down Expand Up @@ -148,7 +162,7 @@ class LinearSpectral {
/// @param[out] result Result, i.e. dvn/dtn
void f1(T& t, std::shared_ptr<la::Vector<T>> u, std::shared_ptr<la::Vector<T>> v,
std::shared_ptr<la::Vector<T>> result) {

// Apply windowing
if (t < period * window_length) {
window = 0.5 * (1.0 - cos(freq * M_PI * t / window_length));
Expand Down Expand Up @@ -240,7 +254,7 @@ class LinearSpectral {
// Store solution at start of time step
kernels::copy<T>(*u_, *u0);
kernels::copy<T>(*v_, *v0);

// Runge-Kutta 4th order step
for (int i = 0; i < 4; i++) {
kernels::copy<T>(*u0, *un);
Expand Down
42 changes: 28 additions & 14 deletions cpp/fenicsx/common/Lossy.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,7 +52,8 @@ void axpy(la::Vector<T>& r, T alpha, const la::Vector<T>& x, const la::Vector<T>
template <typename T, int P>
class LossySpectral {
public:
LossySpectral(std::shared_ptr<mesh::Mesh<T>> Mesh,
LossySpectral(basix::FiniteElement<T> element,
std::shared_ptr<mesh::Mesh<T>> Mesh,
std::shared_ptr<mesh::MeshTags<std::int32_t>> FacetTags,
std::shared_ptr<fem::Function<T>> speedOfSound,
std::shared_ptr<fem::Function<T>> density,
Expand All @@ -79,7 +80,7 @@ class LossySpectral {

// Define function space
V = std::make_shared<fem::FunctionSpace<T>>(
fem::create_functionspace(functionspace_form_forms_a, "u", mesh));
fem::create_functionspace(mesh, element));

// Define field functions
index_map = V->dofmap()->index_map;
Expand All @@ -100,22 +101,35 @@ class LossySpectral {
std::fill(u_.begin(), u_.end(), 1.0);

// Compute exterior facets
std::map<fem::IntegralType,
std::vector<std::pair<std::int32_t, std::span<const std::int32_t>>>> fd;
auto facet_domains = fem::compute_integration_domains(
fem::IntegralType::exterior_facet, *V->mesh()->topology_mutable(),
ft->indices(), mesh->topology()->dim() - 1, ft->values());
for (auto& facet : facet_domains) {
std::vector<std::pair<std::int32_t, std::span<const std::int32_t>>> x;
x.emplace_back(facet.first, std::span(facet.second.data(), facet.second.size()));
fd.insert({fem::IntegralType::exterior_facet, std::move(x)});
}
std::vector<std::int32_t> ft_unique(ft->values().size());
std::copy(ft->values().begin(), ft->values().end(), ft_unique.begin());
std::sort(ft_unique.begin(), ft_unique.end());
auto it = std::unique(ft_unique.begin(), ft_unique.end());
ft_unique.erase(it, ft_unique.end());

std::map<fem::IntegralType, std::vector<std::pair<std::int32_t, std::vector<std::int32_t>>>> fd;
std::map<fem::IntegralType, std::vector<std::pair<std::int32_t, std::span<const std::int32_t>>>> fd_view;

std::vector<std::int32_t> facet_domains;
for (auto& tag : ft_unique) {
facet_domains = fem::compute_integration_domains(
fem::IntegralType::exterior_facet, *V->mesh()->topology_mutable(),
ft->find(tag), mesh->topology()->dim()-1);
fd[fem::IntegralType::exterior_facet].push_back(
{tag, facet_domains});
}

for (auto const& [key, val] : fd) {
for (auto const& [tag, vec] : val) {
fd_view[key].push_back({tag, std::span(vec.data(), vec.size())});
}
}

// Define LHS form
a = std::make_shared<fem::Form<T, T>>(fem::create_form<T, T>(
*form_forms_a, {V},
{{"u", u}, {"c0", c0}, {"rho0", rho0}, {"delta0", delta0}},
{}, fd));
{}, fd_view, {}));

m = std::make_shared<la::Vector<T>>(index_map, bs);
m_ = m->mutable_array();
Expand All @@ -133,7 +147,7 @@ class LossySpectral {
{"c0", c0},
{"rho0", rho0},
{"delta0", delta0}},
{}, fd));
{}, fd_view, {}, {}));
b = std::make_shared<la::Vector<T>>(index_map, bs);
b_ = b->mutable_array();
}
Expand Down
42 changes: 28 additions & 14 deletions cpp/fenicsx/common/Westervelt.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -54,7 +54,8 @@ void axpy(la::Vector<T>& r, T alpha, const la::Vector<T>& x, const la::Vector<T>
template <typename T, int P>
class WesterveltSpectral {
public:
WesterveltSpectral(std::shared_ptr<mesh::Mesh<T>> Mesh,
WesterveltSpectral(basix::FiniteElement<T> element,
std::shared_ptr<mesh::Mesh<T>> Mesh,
std::shared_ptr<mesh::MeshTags<std::int32_t>> FacetTags,
std::shared_ptr<fem::Function<T>> speedOfSound,
std::shared_ptr<fem::Function<T>> density,
Expand Down Expand Up @@ -83,7 +84,7 @@ class WesterveltSpectral {

// Define function space
V = std::make_shared<fem::FunctionSpace<T>>(
fem::create_functionspace(functionspace_form_forms_a, "u", mesh));
fem::create_functionspace(mesh, element));

// Define field functions
index_map = V->dofmap()->index_map;
Expand All @@ -104,22 +105,35 @@ class WesterveltSpectral {
std::fill(u_.begin(), u_.end(), 1.0);

// Compute exterior facets
std::map<fem::IntegralType,
std::vector<std::pair<std::int32_t, std::span<const std::int32_t>>>> fd;
auto facet_domains = fem::compute_integration_domains(
fem::IntegralType::exterior_facet, *V->mesh()->topology_mutable(),
ft->indices(), mesh->topology()->dim() - 1, ft->values());
for (auto& facet : facet_domains) {
std::vector<std::pair<std::int32_t, std::span<const std::int32_t>>> x;
x.emplace_back(facet.first, std::span(facet.second.data(), facet.second.size()));
fd.insert({fem::IntegralType::exterior_facet, std::move(x)});
}
std::vector<std::int32_t> ft_unique(ft->values().size());
std::copy(ft->values().begin(), ft->values().end(), ft_unique.begin());
std::sort(ft_unique.begin(), ft_unique.end());
auto it = std::unique(ft_unique.begin(), ft_unique.end());
ft_unique.erase(it, ft_unique.end());

std::map<fem::IntegralType, std::vector<std::pair<std::int32_t, std::vector<std::int32_t>>>> fd;
std::map<fem::IntegralType, std::vector<std::pair<std::int32_t, std::span<const std::int32_t>>>> fd_view;

std::vector<std::int32_t> facet_domains;
for (auto& tag : ft_unique) {
facet_domains = fem::compute_integration_domains(
fem::IntegralType::exterior_facet, *V->mesh()->topology_mutable(),
ft->find(tag), mesh->topology()->dim()-1);
fd[fem::IntegralType::exterior_facet].push_back(
{tag, facet_domains});
}

for (auto const& [key, val] : fd) {
for (auto const& [tag, vec] : val) {
fd_view[key].push_back({tag, std::span(vec.data(), vec.size())});
}
}

// Define LHS form
a = std::make_shared<fem::Form<T, T>>(fem::create_form<T, T>(
*form_forms_a, {V},
{{"u", u}, {"c0", c0}, {"rho0", rho0}, {"delta0", delta0}, {"u_n", u_n}, {"beta0", beta0}},
{}, fd));
{}, fd_view, {}));

m = std::make_shared<la::Vector<T>>(index_map, bs);
m_ = m->mutable_array();
Expand All @@ -138,7 +152,7 @@ class WesterveltSpectral {
{"rho0", rho0},
{"delta0", delta0},
{"beta0", beta0}},
{}, fd));
{}, fd_view, {}, {}));
b = std::make_shared<la::Vector<T>>(index_map, bs);
b_ = b->mutable_array();
}
Expand Down
4 changes: 2 additions & 2 deletions cpp/fenicsx/examples/linear_planewave2d_1/forms.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,8 +35,8 @@
a = inner(u/rho0/c0/c0, v) * dx(metadata=md)

L = - inner(1/rho0*grad(u_n), grad(v)) * dx(metadata=md) \
+ inner(1/rho0*g, v)*ds(1, metadata=md) \
- inner(1/rho0/c0*v_n, v)*ds(2, metadata=md)
+ inner(1/rho0*g, v) * ds(1, metadata=md) \
- inner(1/rho0/c0*v_n, v) * ds(2, metadata=md)

# Define forms to compute the norm
Norm = inner(u_n, u_n) * dx
Expand Down
26 changes: 19 additions & 7 deletions cpp/fenicsx/examples/linear_planewave2d_1/main.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,10 +44,10 @@ int main(int argc, char* argv[]) {
const int degreeOfBasis = 4;

// Read mesh and mesh tags
auto element = fem::CoordinateElement<T>(mesh::CellType::quadrilateral, 1);
auto coord_element = fem::CoordinateElement<T>(mesh::CellType::quadrilateral, 1);
io::XDMFFile fmesh(MPI_COMM_WORLD, "../mesh.xdmf", "r");
auto mesh = std::make_shared<mesh::Mesh<T>>(
fmesh.read_mesh(element, mesh::GhostMode::none, "planewave_2d_1"));
fmesh.read_mesh(coord_element, mesh::GhostMode::none, "planewave_2d_1"));
mesh->topology()->create_connectivity(1, 2);
auto mt_cell = std::make_shared<mesh::MeshTags<std::int32_t>>(
fmesh.read_meshtags(*mesh, "planewave_2d_1_cells"));
Expand All @@ -68,9 +68,21 @@ int main(int argc, char* argv[]) {
MPI_Reduce(&meshSizeMinLocal, &meshSizeMinGlobal, 1, T_MPI, MPI_MIN, 0, MPI_COMM_WORLD);
MPI_Bcast(&meshSizeMinGlobal, 1, T_MPI, 0, MPI_COMM_WORLD);

// Finite element
basix::FiniteElement element = basix::create_element<T>(
basix::element::family::P, basix::cell::type::quadrilateral, degreeOfBasis,
basix::element::lagrange_variant::gll_warped,
basix::element::dpc_variant::unset, false
);

// Define DG function space for the physical parameters of the domain
basix::FiniteElement element_DG = basix::create_element<T>(
basix::element::family::P, basix::cell::type::quadrilateral, 0,
basix::element::lagrange_variant::gll_warped,
basix::element::dpc_variant::unset, true
);
auto V_DG = std::make_shared<fem::FunctionSpace<T>>(
fem::create_functionspace(functionspace_form_forms_a, "c0", mesh));
fem::create_functionspace(mesh, element_DG));
auto c0 = std::make_shared<fem::Function<T>>(V_DG);
auto rho0 = std::make_shared<fem::Function<T>>(V_DG);

Expand Down Expand Up @@ -113,8 +125,8 @@ int main(int argc, char* argv[]) {
}

// Model
auto model = LinearSpectral<T, 4>(mesh, mt_facet, c0, rho0, sourceFrequency,
sourceAmplitude, speedOfSound);
auto model = LinearSpectral<T, 4>(element, mesh, mt_facet, c0, rho0, sourceFrequency,
sourceAmplitude, speedOfSound);

// Solve
common::Timer tsolve("Solve time");
Expand All @@ -134,12 +146,12 @@ int main(int argc, char* argv[]) {
auto u_n = model.u_sol();

// Output to VTX
dolfinx::io::VTXWriter<T> u_out(mesh->comm(), "output_final.bp", {u_n}, "BP4");
dolfinx::io::VTXWriter<T> u_out(mesh->comm(), "output_final.bp", {u_n}, "bp5");
u_out.write(0.0);

// Check norms
auto Norm = std::make_shared<fem::Form<T>>(
fem::create_form<T, T>(*form_forms_Norm, {}, {{"u_n", u_n}}, {}, {}, mesh));
fem::create_form<T, T>(*form_forms_Norm, {}, {{"u_n", u_n}}, {}, {}, {}, mesh));
T norm = fem::assemble_scalar(*Norm);

if (mpi_rank == 0) {
Expand Down
4 changes: 2 additions & 2 deletions cpp/fenicsx/examples/linear_planewave2d_2/forms.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,7 @@
a = inner(u/rho0/c0/c0, v) * dx(metadata=md)

L = - inner(1/rho0*grad(u_n), grad(v)) * dx(metadata=md) \
+ inner(1/rho0*g, v)*ds(1, metadata=md) \
- inner(1/rho0/c0*v_n, v)*ds(2, metadata=md)
+ inner(1/rho0*g, v) * ds(1, metadata=md) \
- inner(1/rho0/c0*v_n, v) * ds(2, metadata=md)

forms = [a, L]
24 changes: 18 additions & 6 deletions cpp/fenicsx/examples/linear_planewave2d_2/main.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,10 +47,10 @@ int main(int argc, char* argv[]) {
const int degreeOfBasis = 4;

// Read mesh and mesh tags
auto element = fem::CoordinateElement<T>(mesh::CellType::quadrilateral, 1);
auto coord_element = fem::CoordinateElement<T>(mesh::CellType::quadrilateral, 1);
io::XDMFFile fmesh(MPI_COMM_WORLD, "../mesh.xdmf", "r");
auto mesh = std::make_shared<mesh::Mesh<T>>(
fmesh.read_mesh(element, mesh::GhostMode::none, "planewave_2d_4"));
fmesh.read_mesh(coord_element, mesh::GhostMode::none, "planewave_2d_4"));
mesh->topology()->create_connectivity(1, 2);
auto mt_cell = std::make_shared<mesh::MeshTags<std::int32_t>>(
fmesh.read_meshtags(*mesh, "planewave_2d_4_cells"));
Expand All @@ -71,9 +71,21 @@ int main(int argc, char* argv[]) {
MPI_Reduce(&meshSizeMinLocal, &meshSizeMinGlobal, 1, T_MPI, MPI_MIN, 0, MPI_COMM_WORLD);
MPI_Bcast(&meshSizeMinGlobal, 1, T_MPI, 0, MPI_COMM_WORLD);

// Finite element
basix::FiniteElement element = basix::create_element<T>(
basix::element::family::P, basix::cell::type::quadrilateral, degreeOfBasis,
basix::element::lagrange_variant::gll_warped,
basix::element::dpc_variant::unset, false
);

// Define DG function space for the physical parameters of the domain
basix::FiniteElement element_DG = basix::create_element<T>(
basix::element::family::P, basix::cell::type::quadrilateral, 0,
basix::element::lagrange_variant::gll_warped,
basix::element::dpc_variant::unset, true
);
auto V_DG = std::make_shared<fem::FunctionSpace<T>>(
fem::create_functionspace(functionspace_form_forms_a, "c0", mesh));
fem::create_functionspace(mesh, element_DG));
auto c0 = std::make_shared<fem::Function<T>>(V_DG);
auto rho0 = std::make_shared<fem::Function<T>>(V_DG);

Expand Down Expand Up @@ -122,8 +134,8 @@ int main(int argc, char* argv[]) {
}

// Model
auto model = LinearSpectral<T, 4>(mesh, mt_facet, c0, rho0, sourceFrequency,
sourceAmplitude, speedOfSoundWater);
auto model = LinearSpectral<T, 4>(element, mesh, mt_facet, c0, rho0, sourceFrequency,
sourceAmplitude, speedOfSoundWater);

// Solve
model.init();
Expand All @@ -133,7 +145,7 @@ int main(int argc, char* argv[]) {
auto u_n = model.u_sol();

// Output to VTX
dolfinx::io::VTXWriter<T> u_out(mesh->comm(), "output_final.bp", {u_n}, "BP4");
dolfinx::io::VTXWriter<T> u_out(mesh->comm(), "output_final.bp", {u_n}, "bp5");
u_out.write(0.0);
}
}
4 changes: 2 additions & 2 deletions cpp/fenicsx/examples/linear_planewave2d_3/forms.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,8 +35,8 @@
a = inner(u/rho0/c0/c0, v) * dx(metadata=md)

L = - inner(1/rho0*grad(u_n), grad(v)) * dx(metadata=md) \
+ inner(1/rho0*g, v)*ds(1, metadata=md) \
- inner(1/rho0/c0*v_n, v)*ds(2, metadata=md)
+ inner(1/rho0*g, v) * ds(1, metadata=md) \
- inner(1/rho0/c0*v_n, v) * ds(2, metadata=md)

# Define forms to compute the norm
Norm = inner(u_n, u_n) * dx
Expand Down
Loading

0 comments on commit 41e8362

Please sign in to comment.