Skip to content

Issue 1593 part 3 - Reads attenuation database #1628

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
Rohit-Kakodkar wants to merge 8 commits into
base: devel
Choose a base branch
from
Open

Issue 1593 part 3 - Reads attenuation database #1628

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
d23889b
Issue 1593 - 3 - Read the database
Rohit-Kakodkar Feb 3, 2026
f0b8cff
Updated database
Rohit-Kakodkar Feb 3, 2026
bde990c
Fixed tests | added poroelastic test
Rohit-Kakodkar Feb 3, 2026
9e6b3dc
Refactor implementation macros
Rohit-Kakodkar Feb 3, 2026
decc4b7
Merge branch 'devel' into issue-1593-part-3
Rohit-Kakodkar Feb 3, 2026
d73e820
Merge pull request #1629 from PrincetonUniversity/refactor-impl-macro
Rohit-Kakodkar Feb 5, 2026
d579eeb
Merge branch 'devel' into issue-1593-part-3
Rohit-Kakodkar Feb 5, 2026
8bffee0
Fixed compilation
Rohit-Kakodkar Feb 5, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion core/specfem/io/mesh/impl/fortran/dim2/mesh.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -194,7 +194,7 @@ specfem::io::read_2d_mesh(
MEDIUM_TAG(ELASTIC_PSV, ELASTIC_SH, ACOUSTIC, POROELASTIC, ELASTIC_PSV_T,
ELECTROMAGNETIC_TE),
PROPERTY_TAG(ISOTROPIC, ANISOTROPIC, ISOTROPIC_COSSERAT),
ATTENUATION_TAG(NONE)),
ATTENUATION_TAG(NONE, CONSTANT_ISOTROPIC)),
{
total_materials_read += mesh.materials
.get_container<_medium_tag_, _property_tag_,
Expand Down
516 changes: 184 additions & 332 deletions core/specfem/io/mesh/impl/fortran/dim2/read_material_properties.cpp
View file
Open in desktop

Large diffs are not rendered by default.

39 changes: 34 additions & 5 deletions core/specfem/io/mesh/impl/fortran/dim3/read_materials.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -52,7 +52,16 @@ specfem::io::mesh::impl::fortran::dim3::read_materials(std::ifstream &stream,
"materials.");
}

if ((std::abs(Qmu - 9999.0) < 1e-6) || (std::abs(Qmu) < 1e-6)) {
if (!((std::abs(Qmu - 9999.0) < 1e-6) || (std::abs(Qmu) < 1e-6))) {
std::ostringstream error_message;
error_message
<< "Qmu should be set to 9999 or 0 for acoustic materials. "
<< "Found Qmu = " << Qmu << " for material index " << imat
<< "." << "[" << __FILE__ << ":" << __LINE__ << "]\n";
throw std::runtime_error(error_message.str());
}

if ((std::abs(Qkappa - 9999.0) < 1e-6) || (std::abs(Qkappa) < 1e-6)) {

specfem::medium_container::material<
specfem::element::dimension_tag::dim3,
Expand All @@ -66,8 +75,18 @@ specfem::io::mesh::impl::fortran::dim3::read_materials(std::ifstream &stream,
specfem::element::attenuation_tag::none, index,
imat });
} else {
throw std::runtime_error(
"Attenuation not yet supported for acoustic materials in 3D");
specfem::medium_container::material<
specfem::element::dimension_tag::dim3,
specfem::element::medium_tag::acoustic,
specfem::element::property_tag::isotropic,
specfem::element::attenuation_tag::constant_isotropic>
material(rho, vp, Qkappa, static_cast<type_real>(0.0));
const int index = materials.add_material(material);
mapping.push_back(
{ specfem::element::medium_tag::acoustic,
specfem::element::property_tag::isotropic,
specfem::element::attenuation_tag::constant_isotropic, index,
imat });
}
} else if (vs > 0.0) {
// Isotropic elastic material
Expand All @@ -92,8 +111,18 @@ specfem::io::mesh::impl::fortran::dim3::read_materials(std::ifstream &stream,
specfem::element::attenuation_tag::none, index,
imat });
} else {
throw std::runtime_error(
"Attenuation not yet supported for elastic materials in 3D");
specfem::medium_container::material<
specfem::element::dimension_tag::dim3,
specfem::element::medium_tag::elastic,
specfem::element::property_tag::isotropic,
specfem::element::attenuation_tag::constant_isotropic>
material(rho, vs, vp, Qmu, Qkappa, static_cast<type_real>(0.0));
const int index = materials.add_material(material);
mapping.push_back(
{ specfem::element::medium_tag::elastic,
specfem::element::property_tag::isotropic,
specfem::element::attenuation_tag::constant_isotropic, index,
imat });
}

} else {
Expand Down
110 changes: 25 additions & 85 deletions core/specfem/macros/macros_impl/utils.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,15 +34,20 @@

#define _TRANSFORM_INSTANTIATE(s, data, elem) (elem, )

#define _OP_OR(s, state, elem) BOOST_PP_OR(state, elem)
#define _OP_AND(s, state, elem) BOOST_PP_AND(state, elem)

#define _SEQ_FOR_TAGS_2 MEDIUM_TAGS
#define _OP_OR(s, state, elem) BOOST_PP_OR(state, elem)

#define _SEQ_FOR_TAGS_3
#define _ALL_SEQS \
(MEDIUM_TAGS)(ELEMENT_TYPES)(MATERIAL_SYSTEMS)(EDGES)(EDGES_AND_FLUX_SCHEME)

#define _SEQ_FOR_TAGS_4 ELEMENT_TYPES EDGES MATERIAL_SYSTEMS
#define _IS_VALID_SEQ(r, size, seq) \
BOOST_PP_IF( \
BOOST_PP_EQUAL(BOOST_PP_TUPLE_SIZE(BOOST_PP_SEQ_ELEM(0, seq)), size), \
seq, _EMPTY_MACRO())

#define _SEQ_FOR_TAGS_5 EDGES_AND_FLUX_SCHEME
#define _GET_VALID_SEQS(size) \
BOOST_PP_SEQ_FOR_EACH(_IS_VALID_SEQ, size, _ALL_SEQS)

/**
* @brief Declare a variable or instantiante a template based on the type
Expand Down Expand Up @@ -209,93 +214,28 @@
BOOST_PP_EQUAL(_GET_ENUM_ID(enum1), _GET_ENUM_ID(enum2)), \
BOOST_PP_IF(BOOST_PP_EQUAL(_GET_ID(enum1), _GET_ID(enum2)), 1, 0), 0)

/**
* @brief Compare each item in the sequence for a sequence pair of length 2, 3
* and 4.
*/
#define _TYPE_MATCH_2(elem, seq) \
BOOST_PP_IF( \
BOOST_PP_EQUAL(BOOST_PP_TUPLE_SIZE(seq), 2), \
BOOST_PP_IF(_TAG_IN_SEQ(BOOST_PP_TUPLE_ELEM(0, elem), \
BOOST_PP_TUPLE_ELEM(0, seq)), \
BOOST_PP_IF(_TAG_IN_SEQ(BOOST_PP_TUPLE_ELEM(1, elem), \
BOOST_PP_TUPLE_ELEM(1, seq)), \
1, 0), \
0), \
0)

#define _TYPE_MATCH_3(elem, seq) \
BOOST_PP_IF( \
BOOST_PP_EQUAL(BOOST_PP_TUPLE_SIZE(seq), 3), \
BOOST_PP_IF( \
_TAG_IN_SEQ(BOOST_PP_TUPLE_ELEM(0, elem), \
BOOST_PP_TUPLE_ELEM(0, seq)), \
BOOST_PP_IF(_TAG_IN_SEQ(BOOST_PP_TUPLE_ELEM(1, elem), \
BOOST_PP_TUPLE_ELEM(1, seq)), \
BOOST_PP_IF(_TAG_IN_SEQ(BOOST_PP_TUPLE_ELEM(2, elem), \
BOOST_PP_TUPLE_ELEM(2, seq)), \
1, 0), \
0), \
0), \
0)

#define _TYPE_MATCH_4(elem, seq) \
BOOST_PP_IF( \
BOOST_PP_EQUAL(BOOST_PP_TUPLE_SIZE(seq), 4), \
BOOST_PP_IF( \
_TAG_IN_SEQ(BOOST_PP_TUPLE_ELEM(0, elem), \
BOOST_PP_TUPLE_ELEM(0, seq)), \
BOOST_PP_IF( \
_TAG_IN_SEQ(BOOST_PP_TUPLE_ELEM(1, elem), \
BOOST_PP_TUPLE_ELEM(1, seq)), \
BOOST_PP_IF( \
_TAG_IN_SEQ(BOOST_PP_TUPLE_ELEM(2, elem), \
BOOST_PP_TUPLE_ELEM(2, seq)), \
BOOST_PP_IF(_TAG_IN_SEQ(BOOST_PP_TUPLE_ELEM(3, elem), \
BOOST_PP_TUPLE_ELEM(3, seq)), \
1, 0), \
0), \
0), \
0), \
0)

#define _TYPE_MATCH_5(elem, seq) \
BOOST_PP_IF( \
BOOST_PP_EQUAL(BOOST_PP_TUPLE_SIZE(seq), 5), \
BOOST_PP_IF( \
_TAG_IN_SEQ(BOOST_PP_TUPLE_ELEM(0, elem), \
BOOST_PP_TUPLE_ELEM(0, seq)), \
BOOST_PP_IF( \
_TAG_IN_SEQ(BOOST_PP_TUPLE_ELEM(1, elem), \
BOOST_PP_TUPLE_ELEM(1, seq)), \
BOOST_PP_IF( \
_TAG_IN_SEQ(BOOST_PP_TUPLE_ELEM(2, elem), \
BOOST_PP_TUPLE_ELEM(2, seq)), \
BOOST_PP_IF( \
_TAG_IN_SEQ(BOOST_PP_TUPLE_ELEM(3, elem), \
BOOST_PP_TUPLE_ELEM(3, seq)), \
BOOST_PP_IF(_TAG_IN_SEQ(BOOST_PP_TUPLE_ELEM(4, elem), \
BOOST_PP_TUPLE_ELEM(4, seq)), \
1, 0), \
0), \
0), \
0), \
0), \
0)
#define _TYPE_MATCH(elem, seq) \
BOOST_PP_SEQ_FOLD_LEFT( \
_OP_AND, 1, \
BOOST_PP_SEQ_FOR_EACH_I(_TAG_IN_SEQ, elem, BOOST_PP_TUPLE_TO_SEQ(seq)))

/**
* @brief Check if a given tag sequence is in the list of available tag
* sequences.
*/
#define _TAG_EQ(s, tag1, tag2) _CHECK_ENUM(tag1, tag2)
#define _TAG_IN_SEQ(elem, seq) \
BOOST_PP_SEQ_FOLD_LEFT(_OP_OR, 0, BOOST_PP_SEQ_TRANSFORM(_TAG_EQ, elem, seq))

#define _TAG_IN_SEQ(r, elem, i, seq) \
(BOOST_PP_SEQ_FOLD_LEFT( \
_OP_OR, 0, \
BOOST_PP_SEQ_TRANSFORM(_TAG_EQ, BOOST_PP_TUPLE_ELEM(i, elem), seq)))

/**
* Check if a given tag sequence is in the list of available tag sequences,
* write declaration and code block for the sequence if it is in the list.
* Check if a given tag sequence is in the list of available tag
* sequences, write declaration and code block for the sequence if it is
* in the list.
*/
#define _FOR_ONE_TAG_SEQ(s, code, elem) \
BOOST_PP_IF(BOOST_PP_CAT(_TYPE_MATCH_, BOOST_PP_TUPLE_SIZE(elem))( \
elem, BOOST_PP_SEQ_HEAD(code)), \
_CHECK_DECLARE, _EMPTY_MACRO) \
BOOST_PP_IF(_TYPE_MATCH(elem, BOOST_PP_SEQ_HEAD(code)), _CHECK_DECLARE, \
_EMPTY_MACRO) \
(elem, BOOST_PP_SEQ_TAIL(code))
12 changes: 4 additions & 8 deletions core/specfem/macros/material_iterators.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -198,14 +198,10 @@
ATTENUATION_TAG_CONSTANT_ISOTROPIC))( \
(DIMENSION_TAG_DIM2, MEDIUM_TAG_ELASTIC_SH, PROPERTY_TAG_ANISOTROPIC, \
ATTENUATION_TAG_CONSTANT_ISOTROPIC))( \
(DIMENSION_TAG_DIM2, MEDIUM_TAG_ELASTIC_PSV_T, \
PROPERTY_TAG_ISOTROPIC_COSSERAT, ATTENUATION_TAG_CONSTANT_ISOTROPIC))( \
(DIMENSION_TAG_DIM2, MEDIUM_TAG_ACOUSTIC, PROPERTY_TAG_ISOTROPIC, \
ATTENUATION_TAG_CONSTANT_ISOTROPIC))( \
(DIMENSION_TAG_DIM2, MEDIUM_TAG_POROELASTIC, PROPERTY_TAG_ISOTROPIC, \
ATTENUATION_TAG_CONSTANT_ISOTROPIC))( \
(DIMENSION_TAG_DIM2, MEDIUM_TAG_ELECTROMAGNETIC_TE, \
PROPERTY_TAG_ISOTROPIC, ATTENUATION_TAG_CONSTANT_ISOTROPIC))
ATTENUATION_TAG_CONSTANT_ISOTROPIC))

#define MATERIAL_SYSTEMS_DIM3 \
((DIMENSION_TAG_DIM3, MEDIUM_TAG_ELASTIC, PROPERTY_TAG_ISOTROPIC, \
Expand Down Expand Up @@ -319,8 +315,8 @@
BOOST_PP_VARIADIC_TO_SEQ(__VA_ARGS__))

#define FOR_EACH_IN_PRODUCT(seq, ...) \
BOOST_PP_SEQ_FOR_EACH( \
_FOR_ONE_TAG_SEQ, (seq)BOOST_PP_VARIADIC_TO_SEQ(__VA_ARGS__), \
BOOST_PP_CAT(_SEQ_FOR_TAGS_, BOOST_PP_TUPLE_SIZE(seq)))
BOOST_PP_SEQ_FOR_EACH(_FOR_ONE_TAG_SEQ, \
(seq)BOOST_PP_VARIADIC_TO_SEQ(__VA_ARGS__), \
_GET_VALID_SEQS(BOOST_PP_TUPLE_SIZE(seq)))

#include "interface_iterators.hpp"
2 changes: 2 additions & 0 deletions core/specfem/medium/dim2/acoustic/isotropic/material.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,8 @@ struct AttenuationValues<
public:
type_real Qkappa; ///< Attenuation factor for bulk modulus

AttenuationValues() = default;

AttenuationValues(const type_real &Qkappa) : Qkappa(Qkappa) {
if (this->Qkappa <= 0.0) {
throw std::runtime_error(
Expand Down
2 changes: 2 additions & 0 deletions core/specfem/medium/dim2/elastic/isotropic/material.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,8 @@ struct AttenuationValues<
type_real Qkappa; ///< Attenuation factor for bulk modulus
type_real Qmu; ///< Attenuation factor for shear modulus

AttenuationValues() = default;

AttenuationValues(const type_real &Qkappa, const type_real &Qmu)
: Qkappa(Qkappa), Qmu(Qmu) {
if (this->Qkappa <= 0.0 || this->Qmu <= 0.0) {
Expand Down
2 changes: 2 additions & 0 deletions core/specfem/medium/dim2/poroelastic/isotropic/material.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,8 @@ struct AttenuationValues<
}
};

AttenuationValues() = default;

bool operator==(const AttenuationValues &other) const {
return (std::abs(this->Qmu - other.Qmu) < 1e-6);
}
Expand Down
41 changes: 41 additions & 0 deletions core/specfem/mesh/dim2/materials/materials.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,41 @@
#include "materials.hpp"
#include "enumerations/interface.hpp"
#include "specfem/logger.hpp"
#include "specfem/macros.hpp"

void specfem::mesh::materials<specfem::element::dimension_tag::dim2>::print()
const {
std::ostringstream message;
message << "Total number of materials: " << this->n_materials << "\n";
FOR_EACH_IN_PRODUCT(
(DIMENSION_TAG(DIM2),
MEDIUM_TAG(ELASTIC_PSV, ELASTIC_SH, ACOUSTIC, POROELASTIC, ELASTIC_PSV_T,
ELECTROMAGNETIC_TE),
PROPERTY_TAG(ISOTROPIC, ANISOTROPIC, ISOTROPIC_COSSERAT),
ATTENUATION_TAG(NONE, CONSTANT_ISOTROPIC)),
CAPTURE() {
const auto &material_container =
this->get_container<_medium_tag_, _property_tag_,
_attenuation_tag_>();

if (material_container.n_materials > 0) {

message << "Material Type: \n"
<< "\t Medium Tag: "
<< specfem::element::to_string(_medium_tag_) << "\n"
<< "\tProperty Tag: "
<< specfem::element::to_string(_property_tag_) << "\n"
<< "\tAttenuation Tag: "
<< specfem::element::to_string(_attenuation_tag_) << "\n";

for (int i = 0; i < material_container.n_materials; ++i) {
message << "Material Index: " << i << "\n";
message << material_container.element_materials[i].print();
}
}
})

specfem::Logger::info(message.str());

return;
}
Loading