Skip to content

Commit 8ce97e7

Browse files
lsawadelsawade
authored
Merge pull request #1628 from PrincetonUniversity/issue-1593-part-3
Issue 1593 part 3 - Reads attenuation database
2 parents f0c75ab + 8bffee0 commit 8ce97e7

File tree

28 files changed

+1154
-460
lines changed

28 files changed

+1154
-460
lines changed

‎core/specfem/io/mesh/impl/fortran/dim2/mesh.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -195,7 +195,7 @@ specfem::io::read_2d_mesh(
195195
MEDIUM_TAG(ELASTIC_PSV, ELASTIC_SH, ACOUSTIC, POROELASTIC, ELASTIC_PSV_T,
196196
ELECTROMAGNETIC_TE),
197197
PROPERTY_TAG(ISOTROPIC, ANISOTROPIC, ISOTROPIC_COSSERAT),
198-
ATTENUATION_TAG(NONE)),
198+
ATTENUATION_TAG(NONE, CONSTANT_ISOTROPIC)),
199199
{
200200
total_materials_read += mesh.materials
201201
.get_container<_medium_tag_, _property_tag_,

‎core/specfem/io/mesh/impl/fortran/dim2/read_material_properties.cpp‎

Lines changed: 184 additions & 332 deletions
Large diffs are not rendered by default.

‎core/specfem/io/mesh/impl/fortran/dim3/read_materials.cpp‎

Lines changed: 34 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -52,7 +52,16 @@ specfem::io::mesh::impl::fortran::dim3::read_materials(std::ifstream &stream,
5252
"materials.");
5353
}
5454

55-
if ((std::abs(Qmu - 9999.0) < 1e-6) || (std::abs(Qmu) < 1e-6)) {
55+
if (!((std::abs(Qmu - 9999.0) < 1e-6) || (std::abs(Qmu) < 1e-6))) {
56+
std::ostringstream error_message;
57+
error_message
58+
<< "Qmu should be set to 9999 or 0 for acoustic materials. "
59+
<< "Found Qmu = " << Qmu << " for material index " << imat
60+
<< "." << "[" << __FILE__ << ":" << __LINE__ << "]\n";
61+
throw std::runtime_error(error_message.str());
62+
}
63+
64+
if ((std::abs(Qkappa - 9999.0) < 1e-6) || (std::abs(Qkappa) < 1e-6)) {
5665

5766
specfem::medium_container::material<
5867
specfem::element::dimension_tag::dim3,
@@ -66,8 +75,18 @@ specfem::io::mesh::impl::fortran::dim3::read_materials(std::ifstream &stream,
6675
specfem::element::attenuation_tag::none, index,
6776
imat });
6877
} else {
69-
throw std::runtime_error(
70-
"Attenuation not yet supported for acoustic materials in 3D");
78+
specfem::medium_container::material<
79+
specfem::element::dimension_tag::dim3,
80+
specfem::element::medium_tag::acoustic,
81+
specfem::element::property_tag::isotropic,
82+
specfem::element::attenuation_tag::constant_isotropic>
83+
material(rho, vp, Qkappa, static_cast<type_real>(0.0));
84+
const int index = materials.add_material(material);
85+
mapping.push_back(
86+
{ specfem::element::medium_tag::acoustic,
87+
specfem::element::property_tag::isotropic,
88+
specfem::element::attenuation_tag::constant_isotropic, index,
89+
imat });
7190
}
7291
} else if (vs > 0.0) {
7392
// Isotropic elastic material
@@ -92,8 +111,18 @@ specfem::io::mesh::impl::fortran::dim3::read_materials(std::ifstream &stream,
92111
specfem::element::attenuation_tag::none, index,
93112
imat });
94113
} else {
95-
throw std::runtime_error(
96-
"Attenuation not yet supported for elastic materials in 3D");
114+
specfem::medium_container::material<
115+
specfem::element::dimension_tag::dim3,
116+
specfem::element::medium_tag::elastic,
117+
specfem::element::property_tag::isotropic,
118+
specfem::element::attenuation_tag::constant_isotropic>
119+
material(rho, vs, vp, Qmu, Qkappa, static_cast<type_real>(0.0));
120+
const int index = materials.add_material(material);
121+
mapping.push_back(
122+
{ specfem::element::medium_tag::elastic,
123+
specfem::element::property_tag::isotropic,
124+
specfem::element::attenuation_tag::constant_isotropic, index,
125+
imat });
97126
}
98127

99128
} else {

‎core/specfem/macros/macros_impl/utils.hpp‎

Lines changed: 25 additions & 85 deletions
Original file line numberDiff line numberDiff line change
@@ -34,15 +34,20 @@
3434

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

37-
#define _OP_OR(s, state, elem) BOOST_PP_OR(state, elem)
37+
#define _OP_AND(s, state, elem) BOOST_PP_AND(state, elem)
3838

39-
#define _SEQ_FOR_TAGS_2 MEDIUM_TAGS
39+
#define _OP_OR(s, state, elem) BOOST_PP_OR(state, elem)
4040

41-
#define _SEQ_FOR_TAGS_3
41+
#define _ALL_SEQS \
42+
(MEDIUM_TAGS)(ELEMENT_TYPES)(MATERIAL_SYSTEMS)(EDGES)(EDGES_AND_FLUX_SCHEME)
4243

43-
#define _SEQ_FOR_TAGS_4 ELEMENT_TYPES EDGES MATERIAL_SYSTEMS
44+
#define _IS_VALID_SEQ(r, size, seq) \
45+
BOOST_PP_IF( \
46+
BOOST_PP_EQUAL(BOOST_PP_TUPLE_SIZE(BOOST_PP_SEQ_ELEM(0, seq)), size), \
47+
seq, _EMPTY_MACRO())
4448

45-
#define _SEQ_FOR_TAGS_5 EDGES_AND_FLUX_SCHEME
49+
#define _GET_VALID_SEQS(size) \
50+
BOOST_PP_SEQ_FOR_EACH(_IS_VALID_SEQ, size, _ALL_SEQS)
4651

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

212-
/**
213-
* @brief Compare each item in the sequence for a sequence pair of length 2, 3
214-
* and 4.
215-
*/
216-
#define _TYPE_MATCH_2(elem, seq) \
217-
BOOST_PP_IF( \
218-
BOOST_PP_EQUAL(BOOST_PP_TUPLE_SIZE(seq), 2), \
219-
BOOST_PP_IF(_TAG_IN_SEQ(BOOST_PP_TUPLE_ELEM(0, elem), \
220-
BOOST_PP_TUPLE_ELEM(0, seq)), \
221-
BOOST_PP_IF(_TAG_IN_SEQ(BOOST_PP_TUPLE_ELEM(1, elem), \
222-
BOOST_PP_TUPLE_ELEM(1, seq)), \
223-
1, 0), \
224-
0), \
225-
0)
226-
227-
#define _TYPE_MATCH_3(elem, seq) \
228-
BOOST_PP_IF( \
229-
BOOST_PP_EQUAL(BOOST_PP_TUPLE_SIZE(seq), 3), \
230-
BOOST_PP_IF( \
231-
_TAG_IN_SEQ(BOOST_PP_TUPLE_ELEM(0, elem), \
232-
BOOST_PP_TUPLE_ELEM(0, seq)), \
233-
BOOST_PP_IF(_TAG_IN_SEQ(BOOST_PP_TUPLE_ELEM(1, elem), \
234-
BOOST_PP_TUPLE_ELEM(1, seq)), \
235-
BOOST_PP_IF(_TAG_IN_SEQ(BOOST_PP_TUPLE_ELEM(2, elem), \
236-
BOOST_PP_TUPLE_ELEM(2, seq)), \
237-
1, 0), \
238-
0), \
239-
0), \
240-
0)
241-
242-
#define _TYPE_MATCH_4(elem, seq) \
243-
BOOST_PP_IF( \
244-
BOOST_PP_EQUAL(BOOST_PP_TUPLE_SIZE(seq), 4), \
245-
BOOST_PP_IF( \
246-
_TAG_IN_SEQ(BOOST_PP_TUPLE_ELEM(0, elem), \
247-
BOOST_PP_TUPLE_ELEM(0, seq)), \
248-
BOOST_PP_IF( \
249-
_TAG_IN_SEQ(BOOST_PP_TUPLE_ELEM(1, elem), \
250-
BOOST_PP_TUPLE_ELEM(1, seq)), \
251-
BOOST_PP_IF( \
252-
_TAG_IN_SEQ(BOOST_PP_TUPLE_ELEM(2, elem), \
253-
BOOST_PP_TUPLE_ELEM(2, seq)), \
254-
BOOST_PP_IF(_TAG_IN_SEQ(BOOST_PP_TUPLE_ELEM(3, elem), \
255-
BOOST_PP_TUPLE_ELEM(3, seq)), \
256-
1, 0), \
257-
0), \
258-
0), \
259-
0), \
260-
0)
261-
262-
#define _TYPE_MATCH_5(elem, seq) \
263-
BOOST_PP_IF( \
264-
BOOST_PP_EQUAL(BOOST_PP_TUPLE_SIZE(seq), 5), \
265-
BOOST_PP_IF( \
266-
_TAG_IN_SEQ(BOOST_PP_TUPLE_ELEM(0, elem), \
267-
BOOST_PP_TUPLE_ELEM(0, seq)), \
268-
BOOST_PP_IF( \
269-
_TAG_IN_SEQ(BOOST_PP_TUPLE_ELEM(1, elem), \
270-
BOOST_PP_TUPLE_ELEM(1, seq)), \
271-
BOOST_PP_IF( \
272-
_TAG_IN_SEQ(BOOST_PP_TUPLE_ELEM(2, elem), \
273-
BOOST_PP_TUPLE_ELEM(2, seq)), \
274-
BOOST_PP_IF( \
275-
_TAG_IN_SEQ(BOOST_PP_TUPLE_ELEM(3, elem), \
276-
BOOST_PP_TUPLE_ELEM(3, seq)), \
277-
BOOST_PP_IF(_TAG_IN_SEQ(BOOST_PP_TUPLE_ELEM(4, elem), \
278-
BOOST_PP_TUPLE_ELEM(4, seq)), \
279-
1, 0), \
280-
0), \
281-
0), \
282-
0), \
283-
0), \
284-
0)
217+
#define _TYPE_MATCH(elem, seq) \
218+
BOOST_PP_SEQ_FOLD_LEFT( \
219+
_OP_AND, 1, \
220+
BOOST_PP_SEQ_FOR_EACH_I(_TAG_IN_SEQ, elem, BOOST_PP_TUPLE_TO_SEQ(seq)))
221+
285222
/**
286223
* @brief Check if a given tag sequence is in the list of available tag
287224
* sequences.
288225
*/
289226
#define _TAG_EQ(s, tag1, tag2) _CHECK_ENUM(tag1, tag2)
290-
#define _TAG_IN_SEQ(elem, seq) \
291-
BOOST_PP_SEQ_FOLD_LEFT(_OP_OR, 0, BOOST_PP_SEQ_TRANSFORM(_TAG_EQ, elem, seq))
227+
228+
#define _TAG_IN_SEQ(r, elem, i, seq) \
229+
(BOOST_PP_SEQ_FOLD_LEFT( \
230+
_OP_OR, 0, \
231+
BOOST_PP_SEQ_TRANSFORM(_TAG_EQ, BOOST_PP_TUPLE_ELEM(i, elem), seq)))
292232

293233
/**
294-
* Check if a given tag sequence is in the list of available tag sequences,
295-
* write declaration and code block for the sequence if it is in the list.
234+
* Check if a given tag sequence is in the list of available tag
235+
* sequences, write declaration and code block for the sequence if it is
236+
* in the list.
296237
*/
297238
#define _FOR_ONE_TAG_SEQ(s, code, elem) \
298-
BOOST_PP_IF(BOOST_PP_CAT(_TYPE_MATCH_, BOOST_PP_TUPLE_SIZE(elem))( \
299-
elem, BOOST_PP_SEQ_HEAD(code)), \
300-
_CHECK_DECLARE, _EMPTY_MACRO) \
239+
BOOST_PP_IF(_TYPE_MATCH(elem, BOOST_PP_SEQ_HEAD(code)), _CHECK_DECLARE, \
240+
_EMPTY_MACRO) \
301241
(elem, BOOST_PP_SEQ_TAIL(code))

‎core/specfem/macros/material_iterators.hpp‎

Lines changed: 4 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -198,14 +198,10 @@
198198
ATTENUATION_TAG_CONSTANT_ISOTROPIC))( \
199199
(DIMENSION_TAG_DIM2, MEDIUM_TAG_ELASTIC_SH, PROPERTY_TAG_ANISOTROPIC, \
200200
ATTENUATION_TAG_CONSTANT_ISOTROPIC))( \
201-
(DIMENSION_TAG_DIM2, MEDIUM_TAG_ELASTIC_PSV_T, \
202-
PROPERTY_TAG_ISOTROPIC_COSSERAT, ATTENUATION_TAG_CONSTANT_ISOTROPIC))( \
203201
(DIMENSION_TAG_DIM2, MEDIUM_TAG_ACOUSTIC, PROPERTY_TAG_ISOTROPIC, \
204202
ATTENUATION_TAG_CONSTANT_ISOTROPIC))( \
205203
(DIMENSION_TAG_DIM2, MEDIUM_TAG_POROELASTIC, PROPERTY_TAG_ISOTROPIC, \
206-
ATTENUATION_TAG_CONSTANT_ISOTROPIC))( \
207-
(DIMENSION_TAG_DIM2, MEDIUM_TAG_ELECTROMAGNETIC_TE, \
208-
PROPERTY_TAG_ISOTROPIC, ATTENUATION_TAG_CONSTANT_ISOTROPIC))
204+
ATTENUATION_TAG_CONSTANT_ISOTROPIC))
209205

210206
#define MATERIAL_SYSTEMS_DIM3 \
211207
((DIMENSION_TAG_DIM3, MEDIUM_TAG_ELASTIC, PROPERTY_TAG_ISOTROPIC, \
@@ -319,8 +315,8 @@
319315
BOOST_PP_VARIADIC_TO_SEQ(__VA_ARGS__))
320316

321317
#define FOR_EACH_IN_PRODUCT(seq, ...) \
322-
BOOST_PP_SEQ_FOR_EACH( \
323-
_FOR_ONE_TAG_SEQ, (seq)BOOST_PP_VARIADIC_TO_SEQ(__VA_ARGS__), \
324-
BOOST_PP_CAT(_SEQ_FOR_TAGS_, BOOST_PP_TUPLE_SIZE(seq)))
318+
BOOST_PP_SEQ_FOR_EACH(_FOR_ONE_TAG_SEQ, \
319+
(seq)BOOST_PP_VARIADIC_TO_SEQ(__VA_ARGS__), \
320+
_GET_VALID_SEQS(BOOST_PP_TUPLE_SIZE(seq)))
325321

326322
#include "interface_iterators.hpp"

‎core/specfem/medium/dim2/acoustic/isotropic/material.hpp‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -19,6 +19,8 @@ struct AttenuationValues<
1919
public:
2020
type_real Qkappa; ///< Attenuation factor for bulk modulus
2121

22+
AttenuationValues() = default;
23+
2224
AttenuationValues(const type_real &Qkappa) : Qkappa(Qkappa) {
2325
if (this->Qkappa <= 0.0) {
2426
throw std::runtime_error(

‎core/specfem/medium/dim2/elastic/isotropic/material.hpp‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -22,6 +22,8 @@ struct AttenuationValues<
2222
type_real Qkappa; ///< Attenuation factor for bulk modulus
2323
type_real Qmu; ///< Attenuation factor for shear modulus
2424

25+
AttenuationValues() = default;
26+
2527
AttenuationValues(const type_real &Qkappa, const type_real &Qmu)
2628
: Qkappa(Qkappa), Qmu(Qmu) {
2729
if (this->Qkappa <= 0.0 || this->Qmu <= 0.0) {

‎core/specfem/medium/dim2/poroelastic/isotropic/material.hpp‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -27,6 +27,8 @@ struct AttenuationValues<
2727
}
2828
};
2929

30+
AttenuationValues() = default;
31+
3032
bool operator==(const AttenuationValues &other) const {
3133
return (std::abs(this->Qmu - other.Qmu) < 1e-6);
3234
}

‎core/specfem/mesh/dim2/materials/materials.cpp‎

Lines changed: 41 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,41 @@
1+
#include "materials.hpp"
2+
#include "enumerations/interface.hpp"
3+
#include "specfem/logger.hpp"
4+
#include "specfem/macros.hpp"
5+
6+
void specfem::mesh::materials<specfem::element::dimension_tag::dim2>::print()
7+
const {
8+
std::ostringstream message;
9+
message << "Total number of materials: " << this->n_materials << "\n";
10+
FOR_EACH_IN_PRODUCT(
11+
(DIMENSION_TAG(DIM2),
12+
MEDIUM_TAG(ELASTIC_PSV, ELASTIC_SH, ACOUSTIC, POROELASTIC, ELASTIC_PSV_T,
13+
ELECTROMAGNETIC_TE),
14+
PROPERTY_TAG(ISOTROPIC, ANISOTROPIC, ISOTROPIC_COSSERAT),
15+
ATTENUATION_TAG(NONE, CONSTANT_ISOTROPIC)),
16+
CAPTURE() {
17+
const auto &material_container =
18+
this->get_container<_medium_tag_, _property_tag_,
19+
_attenuation_tag_>();
20+
21+
if (material_container.n_materials > 0) {
22+
23+
message << "Material Type: \n"
24+
<< "\t Medium Tag: "
25+
<< specfem::element::to_string(_medium_tag_) << "\n"
26+
<< "\tProperty Tag: "
27+
<< specfem::element::to_string(_property_tag_) << "\n"
28+
<< "\tAttenuation Tag: "
29+
<< specfem::element::to_string(_attenuation_tag_) << "\n";
30+
31+
for (int i = 0; i < material_container.n_materials; ++i) {
32+
message << "Material Index: " << i << "\n";
33+
message << material_container.element_materials[i].print();
34+
}
35+
}
36+
})
37+
38+
specfem::Logger::info(message.str());
39+
40+
return;
41+
}

0 commit comments

Comments
 (0)