Refactor kernels within domain class.

[Rohit-Kakodkar]

No description provided.

[Rohit-Kakodkar]

Okay seems like this is tougher to do then I thought. What I want is be able to define different mediums, material_systems, and element_types at compile time. As we add more medium_tags, property_tags and boundary_tags things can be updated in these arrays.

constexpr auto mediums() {
  constexpr int total_mediums = 2;
  return std::array<medium_tag, total_mediums>{ medium_tag::elastic,
                                                medium_tag::acoustic };
}

constexpr auto material_systems() {
  constexpr int total_material_systems = 3;
  return std::array<std::pair<medium_tag, property_tag>,
                    total_material_systems>{
    std::make_pair(medium_tag::elastic, property_tag::isotropic),
    std::make_pair(medium_tag::elastic, property_tag::anisotropic),
    std::make_pair(medium_tag::acoustic, property_tag::isotropic)
  };
}

constexpr auto element_types() {
  constexpr int total_element_types = 8;
  return std::array<std::tuple<medium_tag, property_tag, total_element_types>,
                    N> {
    std::make_tuple(medium_tag::elastic, property_tag::isotropic,
                    boundary_tag::none),
        std::make_tuple(medium_tag::elastic, property_tag::isotropic,
                        boundary_tag::stacey),
        std::make_tuple(medium_tag::elastic, property_tag::anisotropic,
                        boundary_tag::none),
        std::make_tuple(medium_tag::elastic, property_tag::anisotropic,
                        boundary_tag::stacey),
        std::make_tuple(medium_tag::acoustic, property_tag::isotropic,
                        boundary_tag::none),
        std::make_tuple(medium_tag::acoustic, property_tag::isotropic,
                        boundary_tag::acoustic_free_surface),
        std::make_tuple(medium_tag::acoustic, property_tag::isotropic,
                        boundary_tag::stacey),
        std::make_tuple(medium_tag::acoustic, property_tag::isotropic,
                        boundary_tag::composite_stacey_dirichlet)
  }
}

However, this is only step one. Since our functions take tags at compile time, I'd need a compile time iteration loop - which is hard and can be tough to read. So I started thinking about the interface I wanted first. What I want is something like this:

template <specfem::dimension::type DimensionType,
          specfem::wavefield::simulation_field WavefieldType,
          specfem::medium::medium_tag MediumTag, typename qp_type>
inline void update_wavefields(const int istep) {

  specfem::element::for_each_element_types(
      [&]<specfem::element::medium_tag KernelMediumTag,
          specfem::element::property_tag KernelPropertyTag,
          specfem::element::boundary_tag KernelBoundaryTag>() {
        specfem::domain::compute_source_interaction<
            WavefieldType, DimensionType, KernelMediumTag, KernelBoundaryTag, NGLL>(istep);
        specfem::domain::compute_stiffness_interaction<
            WavefieldType, DimensionType, KernelMediumTag, KernelPropertyTag, KernelBoundaryTag,
            NGLL>();
      })
      .where(MediumTag)
      .execute();

  specfem::element::for_each_medium_types(
      [&]<specfem::element::medium_tag KernelMediumTag>() {
        specfem::domain::divide_mass_matrix<WavefieldType, DimensionType,
                                            KernelMediumTag, NGLL>();
      })
      .where(MediumTag)
      .execute();

  return;
}

Whats happening here is I calling the compute_source_interaction and compute_stiffness_interaction functions on all constexpr element_types where medium matches the MediumTag. execute() finally invokes the functions.

The implementation I have right now is something like this

namespace impl {

struct all_tag {};

} // namespace impl

template <typename Functor> class for_each_element_type {

public:
  for_each_element_type(Functor functor) : functor(functor) {}

  for_each_element_type &where(const specfem::element::medium_tag medium) {
    this->medium = medium;
    return *this;
  }

  for_each_element_type &where(const specfem::element::property_tag property) {
    this->property = property;
    return *this;
  }

  for_each_element_type &where(const specfem::element::boundary_tag boundary) {
    this->boundary = boundary;
    return *this;
  }

  void execute() {
    constexpr auto ntypes = element_types.size();
    invoke_functor(std::make_index_sequence<ntypes>());
  }

private:
  Functor functor;
  std::variant<specfem::element::medium_tag, impl::all_tag> medium =
      impl::all_tag{};
  std::variant<specfem::element::property_tag, impl::all_tag> property =
      impl::all_tag{};
  std::variant<specfem::element::boundary_tag, impl::all_tag> boundary =
      impl::all_tag{};
  constexpr auto element_types = element_types();

  template <std::size_t Index> void invoke_functor() {
    constexpr auto element_type = element_types[Index];
    constexpr auto medium = std::get<0>(element_type);
    constexpr auto property = std::get<1>(element_type);
    constexpr auto boundary = std::get<2>(element_type);
    if ((medium == this->medium || this->medium == impl::all_tag{}) &&
        (property == this->property || this->property == impl::all_tag{}) &&
        (boundary == this->boundary || this->boundary == impl::all_tag{})) {
      functor.template operator()<medium, property, boundary>();
    }
  }

  template <std::integral_sequence<std::size_t...> Sequence>
  void invoke_functor(std::index_sequence<Sequence...>) {
    (invoke_functor<Sequence>(), ...);
  }
}

I think this should do the trick, and will reduce manual work tremendously throughout the code. However it also would require considerable refactor. So it makes sense we do this after anisotopy.