Add example for 2D penetrable circular scatterer

cpp/fenicsx-sf-naive/examples/linear_penetrable2d_1/CMakeLists.txt

@@ -0,0 +1,32 @@
+cmake_minimum_required(VERSION 3.16)
+
+set(PROJECT_NAME linear_penetrable2d_1)
+project(${PROJECT_NAME} LANGUAGES C CXX)
+
+# Set standard
+set(CMAKE_CXX_STANDARD 20)
+set(CMAKE_CXX_STANDARD_REQUIRED ON)
+set(CMAKE_CXX_EXTENSIONS OFF)
+
+# Set flags
+set(CMAKE_CXX_FLAGS "-Ofast -march=native ${CMAKE_CXX_FLAGS} -Wall")
+set(CMAKE_C_FLAGS "-Ofast -march=native ${CMAKE_C_FLAGS} -Wall")
+
+if(NOT TARGET dolfinx)
+  find_package(DOLFINX REQUIRED)
+endif()
+
+add_custom_command(
+  OUTPUT ${CMAKE_CURRENT_SOURCE_DIR}/forms.c
+  COMMAND ffcx ${CMAKE_CURRENT_SOURCE_DIR}/forms.py -o
+          ${CMAKE_CURRENT_SOURCE_DIR}/
+  VERBATIM
+  DEPENDS forms.py
+  COMMENT "Compile forms.py using FFCx")
+
+set(CMAKE_INCLUDE_CURRENT_DIR ON)
+
+include_directories("../../common" ".")
+
+add_executable(${PROJECT_NAME} forms.c main.cpp)
+target_link_libraries(${PROJECT_NAME} dolfinx)

cpp/fenicsx-sf-naive/examples/linear_penetrable2d_1/forms.py

@@ -0,0 +1,39 @@
+import basix
+from basix.ufl import element
+from ufl import (Coefficient, FunctionSpace, Mesh, TestFunction,
+                 ds, dx, inner)
+
+P = 4  # Degree of polynomial basis
+Q = 5  # Number of quadrature points
+
+# Define mesh and finite element
+geom_order = 2
+coord_element = element("Lagrange", "quadrilateral", geom_order, shape=(2, ))
+mesh = Mesh(coord_element)
+e = element(basix.ElementFamily.P, basix.CellType.quadrilateral, P,
+    basix.LagrangeVariant.gll_warped)
+e_DG = element(basix.ElementFamily.P, basix.CellType.quadrilateral, 0,
+    basix.LagrangeVariant.gll_warped, basix.DPCVariant.unset, True)
+
+# Define function spaces
+V = FunctionSpace(mesh, e)
+V_DG = FunctionSpace(mesh, e_DG)
+
+c0 = Coefficient(V_DG)
+rho0 = Coefficient(V_DG)
+
+u = Coefficient(V)
+u_n = Coefficient(V)
+v_n = Coefficient(V)
+g = Coefficient(V)
+v = TestFunction(V)
+
+# Map from quadrature points to basix quadrature degree
+qdegree = {3: 4, 4: 5, 5: 6, 6: 8, 7: 10, 8: 12, 9: 14, 10: 16}
+md = {"quadrature_rule": "GLL", "quadrature_degree": qdegree[Q]}
+
+# Define forms
+a = inner(u/rho0/c0/c0, v) * dx(metadata=md)
+
+L = inner(1/rho0*g, v) * ds(1, metadata=md) \
+    - inner(1/rho0/c0*v_n, v) * ds(2, metadata=md)

cpp/fenicsx-sf-naive/examples/linear_penetrable2d_1/main.cpp

@@ -0,0 +1,168 @@
+//
+// Linear solver for the 2D circular penetrable scatterer
+// - circular scatterer
+// - wavelength < scatterer radius
+// ======================================================
+// Copyright (C) 2024 Adeeb Arif Kor
+
+#include "Linear.hpp"
+#include "forms.h"
+
+#include <cmath>
+#include <dolfinx.h>
+#include <dolfinx/fem/Constant.h>
+#include <dolfinx/io/XDMFFile.h>
+#include <iomanip>
+#include <iostream>
+
+#define T_MPI MPI_DOUBLE
+using T = double;
+
+int main(int argc, char* argv[]) {
+  dolfinx::init_logging(argc, argv);
+  PetscInitialize(&argc, &argv, nullptr, nullptr);
+  {
+    // MPI
+    int mpi_rank, mpi_size;
+    MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank);
+    MPI_Comm_size(MPI_COMM_WORLD, &mpi_size);
+
+    // Material parameters
+    const T speedOfSound1 = 1500; // (m/s)
+    const T density1 = 1000;      // (kg/m^3)
+    const T speedOfSound2 = 3500; // (m/s)
+    const T density2 = 1900;      // (kg/m^3)
+
+    // Source parameters
+    const T sourceFrequency = 2000;      // (Hz)
+    const T sourceSpeed = 1.0;  // (m/s)
+    const T sourceAmplitude = density1 * speedOfSound1 * sourceSpeed; // (Pa)
+    const T period = 1 / sourceFrequency; // (s)
+
+    // Domain parameters
+    const T wavelength = speedOfSound1 / sourceFrequency;
+    const T scattererRadius = 1.0 * wavelength;
+    const T domainScale = 10.0;
+    const T domainLength =  1.5 * domainScale * scattererRadius + 2 * scattererRadius; // (m)
+
+    // FE parameters
+    const int degreeOfBasis = 4;
+
+    // Read mesh and mesh tags
+    int geom_order = 2;
+    auto coord_element = fem::CoordinateElement<T>(mesh::CellType::quadrilateral, geom_order);
+    io::XDMFFile fmesh(MPI_COMM_WORLD, "../mesh.xdmf", "r");
+    auto mesh = std::make_shared<mesh::Mesh<T>>(
+        fmesh.read_mesh(coord_element, mesh::GhostMode::none, "circular_scatterer_penetrable_2d_1"));
+    mesh->topology()->create_connectivity(1, 2);
+    auto mt_cell = std::make_shared<mesh::MeshTags<std::int32_t>>(
+        fmesh.read_meshtags(*mesh, "circular_scatterer_penetrable_2d_1_cells"));
+    auto mt_facet = std::make_shared<mesh::MeshTags<std::int32_t>>(
+        fmesh.read_meshtags(*mesh, "circular_scatterer_penetrable_2d_1_facets"));
+
+    // Mesh parameters
+    const int tdim = mesh->topology()->dim();
+    const int num_cell = mesh->topology()->index_map(tdim)->size_local();
+    std::vector<int> num_cell_range(num_cell);
+    std::iota(num_cell_range.begin(), num_cell_range.end(), 0.0);
+    std::vector<T> mesh_size_local = mesh::h(*mesh, num_cell_range, tdim);
+    std::vector<T>::iterator min_mesh_size_local
+        = std::min_element(mesh_size_local.begin(), mesh_size_local.end());
+    int mesh_size_local_idx = std::distance(mesh_size_local.begin(), min_mesh_size_local);
+    T meshSizeMinLocal = mesh_size_local.at(mesh_size_local_idx);
+    T meshSizeMinGlobal;
+    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(mesh, element_DG));
+    auto c0 = std::make_shared<fem::Function<T>>(V_DG);
+    auto rho0 = std::make_shared<fem::Function<T>>(V_DG);
+
+    auto cells_1 = mt_cell->find(1);
+    auto cells_2 = mt_cell->find(2);
+
+    std::span<T> c0_ = c0->x()->mutable_array();
+    std::for_each(cells_1.begin(), cells_1.end(),
+                  [&](std::int32_t& i) { c0_[i] = speedOfSound1; });
+    std::for_each(cells_2.begin(), cells_2.end(),
+                  [&](std::int32_t& i) { c0_[i] = speedOfSound2; });
+    c0->x()->scatter_fwd();
+
+    std::span<T> rho0_ = rho0->x()->mutable_array();
+    std::for_each(cells_1.begin(), cells_1.end(),
+                  [&](std::int32_t& i) { rho0_[i] = density1; });
+    std::for_each(cells_2.begin(), cells_2.end(),
+                  [&](std::int32_t& i) { rho0_[i] = density2; });
+    rho0->x()->scatter_fwd();
+
+    // Temporal parameters
+    const T CFL = 0.9;
+    T timeStepSize = CFL * meshSizeMinGlobal / (speedOfSound2 * degreeOfBasis * degreeOfBasis);
+    const int stepPerPeriod = period / timeStepSize + 1;
+    timeStepSize = period / stepPerPeriod;
+    const T startTime = 0.0;
+    const T finalTime = domainLength / speedOfSound1 + 8.0 / sourceFrequency;
+    const int numberOfStep = (finalTime - startTime) / timeStepSize + 1;
+
+    if (mpi_rank == 0) {
+      std::cout << "Problem type: Planewave 2D"
+                << "\n";
+      std::cout << "Scatterer radius: " << scattererRadius << "\n";
+      std::cout << "Wavelength: " << wavelength << "\n";
+      std::cout << "Speed of sound (1): " << speedOfSound1 << "\n";
+      std::cout << "Speed of sound (2): " << speedOfSound2 << "\n";
+      std::cout << "Density (1): " << density1 << "\n";
+      std::cout << "Density (2): " << density2 << "\n";
+      std::cout << "Source frequency: " << sourceFrequency << "\n";
+      std::cout << "Source amplitude: " << sourceAmplitude << "\n";
+      std::cout << "Domain length: " << domainLength << "\n";
+      std::cout << "Polynomial basis degree: " << degreeOfBasis << "\n";
+      std::cout << "Minimum mesh size: ";
+      std::cout << std::setprecision(2) << meshSizeMinGlobal << "\n";
+      std::cout << "CFL number: " << CFL << "\n";
+      std::cout << "Time step size: " << timeStepSize << "\n";
+      std::cout << "Final time: " << finalTime << "\n";
+      std::cout << "Number of steps per period: " << stepPerPeriod << "\n";
+      std::cout << "Total number of steps: " << numberOfStep << "\n";
+    }
+
+    // Model
+    auto model = LinearSpectral2D<T, degreeOfBasis>(element, mesh, mt_facet, c0, rho0, sourceFrequency,
+                                                    sourceAmplitude, speedOfSound1);
+
+    // Solve
+    common::Timer tsolve("Solve time");
+
+    model.init();
+
+    tsolve.start();
+    model.rk4(startTime, finalTime, timeStepSize);
+    tsolve.stop();
+
+    if (mpi_rank == 0) {
+      std::cout << "Solve time: " << tsolve.elapsed()[0] << std::endl;
+      std::cout << "Time per step: " << tsolve.elapsed()[0] / numberOfStep << std::endl;
+    }
+
+    // Final solution
+    auto u_n = model.u_sol();
+
+    // Output to VTX
+    dolfinx::io::VTXWriter<T> u_out(mesh->comm(), "output_final.bp", {u_n}, "bp5");
+    u_out.write(0.0);
+  }
+}

cpp/fenicsx-sf-naive/examples/linear_penetrable2d_1/mesh.xdmf

@@ -0,0 +1,32 @@
+<?xml version="1.0"?>
+<!DOCTYPE Xdmf SYSTEM "Xdmf.dtd" []>
+<Xdmf Version="3.0" xmlns:xi="https://www.w3.org/2001/XInclude">
+  <Domain>
+    <Grid Name="circular_scatterer_penetrable_2d_1" GridType="Uniform">
+      <Topology TopologyType="Quadrilateral_9" NumberOfElements="18556" NodesPerElement="9">
+        <DataItem Dimensions="18556 9" NumberType="Int" Format="HDF">mesh.h5:/Mesh/circular_scatterer_penetrable_2d_1/topology</DataItem>
+      </Topology>
+      <Geometry GeometryType="XY">
+        <DataItem Dimensions="74769 2" Format="HDF">mesh.h5:/Mesh/circular_scatterer_penetrable_2d_1/geometry</DataItem>
+      </Geometry>
+    </Grid>
+    <Grid Name="circular_scatterer_penetrable_2d_1_cells" GridType="Uniform">
+      <xi:include xpointer="xpointer(/Xdmf/Domain/Grid[@Name='circular_scatterer_penetrable_2d_1']/Geometry)" />
+      <Topology TopologyType="Quadrilateral_9" NumberOfElements="18556" NodesPerElement="9">
+        <DataItem Dimensions="18556 9" NumberType="Int" Format="HDF">mesh.h5:/MeshTags/circular_scatterer_penetrable_2d_1_cells/topology</DataItem>
+      </Topology>
+      <Attribute Name="circular_scatterer_penetrable_2d_1_cells" AttributeType="Scalar" Center="Cell">
+        <DataItem Dimensions="18556 1" Format="HDF">mesh.h5:/MeshTags/circular_scatterer_penetrable_2d_1_cells/Values</DataItem>
+      </Attribute>
+    </Grid>
+    <Grid Name="circular_scatterer_penetrable_2d_1_facets" GridType="Uniform">
+      <xi:include xpointer="xpointer(/Xdmf/Domain/Grid[@Name='circular_scatterer_penetrable_2d_1']/Geometry)" />
+      <Topology TopologyType="Edge_3" NumberOfElements="544" NodesPerElement="3">
+        <DataItem Dimensions="544 3" NumberType="Int" Format="HDF">mesh.h5:/MeshTags/circular_scatterer_penetrable_2d_1_facets/topology</DataItem>
+      </Topology>
+      <Attribute Name="circular_scatterer_penetrable_2d_1_facets" AttributeType="Scalar" Center="Cell">
+        <DataItem Dimensions="544 1" Format="HDF">mesh.h5:/MeshTags/circular_scatterer_penetrable_2d_1_facets/Values</DataItem>
+      </Attribute>
+    </Grid>
+  </Domain>
+</Xdmf>