v0.14.0
ScalingMethod.cpp
Go to the documentation of this file.
1 /**
2  * @file ScalingMethod.cpp
3  * @brief
4  * @version 0.1
5  * @date 2022-08-12
6  *
7  * @copyright Copyright (c) 2022
8  *
9  */
10 
11 #include <MoFEM.hpp>
12 
13 namespace MoFEM {
14 
15 double ScalingMethod::getScale(const double time) {
16  THROW_MESSAGE("getScale not implemented");
17 }
18 
19 const std::string TimeScale::defaultDelimiter =
20  "(\\s*,\\s*|\\s+)"; ///< comma or space
21 
22 TimeScale::TimeScale(std::string file_name, bool error_if_file_not_given)
23  : TimeScale(file_name, defaultDelimiter, error_if_file_not_given) {}
24 
25 TimeScale::TimeScale(std::string file_name, std::string delimiter,
26  bool error_if_file_not_given)
27  : fileName(file_name), errorIfFileNotGiven(error_if_file_not_given) {
28  CHK_THROW_MESSAGE(timeData(file_name, delimiter),
29  "Error in reading time data");
30 }
31 
32 MoFEMErrorCode TimeScale::timeData(std::string fileName,
33  std::string delimiter) {
35  MOFEM_LOG_CHANNEL("WORLD");
36  // Set the argument found flag to false as default
37  PetscBool arg_found = PETSC_FALSE;
38  char time_file_name[255] = {'\0'};
39  // Check to see if a filename has been provided
40  if (fileName.empty()) {
41  // If no filename, look for command line argument
42  CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, fileNameFlag.c_str(),
43  time_file_name, 255, &arg_found);
44  if (arg_found) {
45  fileName = std::string(time_file_name);
46  }
47  } else {
48  // Set the command line flag to true for correct flow control using provided
49  // filename
50  arg_found = PETSC_TRUE;
51  }
52  if (!arg_found && fileName.empty() && errorIfFileNotGiven) {
53  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
54  "*** ERROR %s (time_data FILE NEEDED)", fileName.c_str());
55  } else if (!arg_found && fileName.empty() && !errorIfFileNotGiven) {
56  MOFEM_LOG_C("WORLD", Sev::warning,
57  "The %s file not provided. Loading scaled with time.",
58  fileName.c_str());
59  scalingMethod = [this](double time) { return this->getLinearScale(time); };
61  }
62 
63  std::ifstream in_file_stream(fileName);
64  if (!in_file_stream.is_open() && errorIfFileNotGiven) {
65  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
66  "*** ERROR data file < %s > open unsuccessful", fileName.c_str());
67  } else if (!in_file_stream.is_open() && !errorIfFileNotGiven) {
68  MOFEM_LOG("WORLD", Sev::warning)
69  << "*** Warning data file " << fileName
70  << " open unsuccessful. Using linear time scaling.";
71  scalingMethod = [this](double time) { return this->getLinearScale(time); };
73  }
74 
75  in_file_stream.seekg(0);
76  std::string line;
77  double time = 0.0, value = 0.0;
78  tSeries[time] = value;
79 
80  std::regex rgx(delimiter.c_str());
81  std::sregex_token_iterator end;
82  while (std::getline(in_file_stream, line)) {
83  std::sregex_token_iterator iter(line.begin(), line.end(), rgx, -1);
84  auto time_str = iter;
85  auto value_str = ++iter;
86  if (time_str == end || value_str == end) {
87  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
88  "*** ERROR read data file error (check input time data file) ");
89  }
90  time = std::stod(time_str->str());
91  value = std::stod(value_str->str());
92  tSeries[time] = value;
93  }
94  in_file_stream.close();
95  scalingMethod = [this](double time) { return this->getScaleFromData(time); };
96  if (in_file_stream.is_open()) {
97  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
98  "*** ERROR file close unsuccessful");
99  }
101 }
102 
103 double TimeScale::getLinearScale(const double time) { return time; }
104 
105 double TimeScale::getScaleFromData(const double time) {
106  if (tSeries.empty())
107  return 0.;
108 
109  auto it = tSeries.find(time);
110  if (it != tSeries.end()) {
111  return it->second;
112  } else {
113  auto it = tSeries.lower_bound(time);
114  if (it == tSeries.end()) {
115  return (--it)->second;
116  }
117  if (it == tSeries.begin()) {
118  return it->second;
119  }
120  auto upper = *(it);
121  it--;
122  if (it == tSeries.end()) {
123  return upper.second;
124  }
125  auto lower = *it;
126  double t = (time - lower.first) / (upper.first - lower.first);
127  double scale1 = upper.second;
128  double scale0 = lower.second;
129  return scale0 + t * (scale1 - scale0);
130  }
131 }
132 
133 double TimeScale::getScale(const double time) { return scalingMethod(time); }
134 
135 // Vector scale
136 
137 // template <int SPACE_DIM> TimeScaleVector<SPACE_DIM>::TimeScaleVector() {}
138 template <int SPACE_DIM>
140  bool error_if_file_not_given)
141  : readFile(0), debug(0), nAme(name),
142  errorIfFileNotGiven(error_if_file_not_given) {
143  CHK_THROW_MESSAGE(timeData(), "Error in reading time data");
144 }
145 
146 template <int SPACE_DIM>
148  bool error_if_file_not_given)
149  : readFile(0), debug(0), errorIfFileNotGiven(error_if_file_not_given) {
150  nAme = name + std::to_string(ms_id);
151  CHK_THROW_MESSAGE(timeData(), "Error in reading time data");
152 }
153 
156 
157  char time_file_name[255];
158  PetscBool flg = PETSC_FALSE;
159  if (!nAme.empty())
160  CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, nAme.c_str(),
161  time_file_name, 255, &flg);
162 
163  if (!flg) {
164  if (errorIfFileNotGiven)
165  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
166  "*** ERROR %s (time_data FILE NEEDED)", nAme.c_str());
168  }
169 
170  FILE *time_data = fopen(time_file_name, "r");
171  if (time_data == NULL) {
172  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
173  "*** ERROR data file < %s > open unsuccessful", time_file_name);
174  }
175  double no1 = 0.0;
177  (tSeries[no1])(i) = 0.;
178  while (!feof(time_data)) {
179  FTensor::Tensor1<double, 3> no2{0., 0., 0.};
180  int n =
181  fscanf(time_data, "%lf %lf %lf %lf", &no1, &no2(0), &no2(1), &no2(2));
182  if (n < 0) {
183  fgetc(time_data);
184  continue;
185  }
186  if (n != 4) {
187  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
188  "*** ERROR read data file error (check input time data file) "
189  "{ n = %d }",
190  n);
191  }
192  (tSeries[no1])(i) = no2(i);
193  }
194  int r = fclose(time_data);
195 
196  MOFEM_LOG_CHANNEL("WORLD");
197  for (auto &[ts, vec] : tSeries) {
198  MOFEM_TAG_AND_LOG_C("WORLD", Sev::verbose, "TimeScaleVector",
199  "** read vector %3.2e time %3.2e %3.2e %3.2e", ts,
200  vec(0), vec(1), vec(2));
201  }
202  MOFEM_LOG_CHANNEL("WORLD");
203 
204  if (r != 0) {
205  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
206  "*** ERROR file close unsuccessful");
207  }
208  readFile = 1;
209 
210  if (readFile == 1)
211  scalingMethod = [this](double time) {
212  return this->getVectorFromData(time);
213  };
214 
216 }
217 
218 template <int SPACE_DIM>
221  return scalingMethod(time);
222 }
223 
224 template <int SPACE_DIM>
227 
230  FTensor::Tensor1<double, 3> acc0 = tSeries.begin()->second;
231 
234 
235  double t0 = 0, t1, dt;
236  for (auto &[ts, vec] : tSeries) {
237  if (ts > time) {
238  t1 = ts;
239  dt = time - t0;
240  acc(I) = acc0(I) + ((vec(I) - acc0(I)) / (t1 - t0)) * dt;
241  break;
242  }
243  t0 = ts;
244  acc0 = vec;
245  acc = acc0;
246  }
247  Nf(i) = acc(i);
248  return Nf;
249 }
250 
251 template class TimeScaleVector<3>;
252 template class TimeScaleVector<2>;
253 
254 } // namespace MoFEM
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
FTensor::Tensor1< double, 3 >
MOFEM_LOG_CHANNEL
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:609
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::TimeScaleVector::getVectorFromData
virtual FTensor::Tensor1< double, SPACE_DIM > getVectorFromData(const double time)
Definition: ScalingMethod.cpp:226
MoFEM.hpp
THROW_MESSAGE
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
Definition: definitions.h:574
sdf.r
int r
Definition: sdf.py:8
I
constexpr IntegrationType I
Definition: operators_tests.cpp:31
MoFEM::TimeScaleVector::timeData
MoFEMErrorCode timeData()
Definition: ScalingMethod.cpp:154
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::TimeScale::TimeScale
TimeScale(std::string file_name="", bool error_if_file_not_given=false)
TimeScale constructor.
Definition: ScalingMethod.cpp:22
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
MoFEM::TimeScale
Force scale operator for reading two columns.
Definition: ScalingMethod.hpp:32
MoFEM::TimeScale::tSeries
std::map< double, double > tSeries
Definition: ScalingMethod.hpp:82
MoFEM::TimeScale::fileName
std::string fileName
Definition: ScalingMethod.hpp:83
MoFEM::TimeScale::timeData
MoFEMErrorCode timeData(std::string fileName, std::string delimiter)
Definition: ScalingMethod.cpp:32
t
constexpr double t
plate stiffness
Definition: plate.cpp:58
MoFEM::TimeScaleVector
Force scale operator for reading four columns (time and vector)
Definition: ScalingMethod.hpp:97
debug
static const bool debug
Definition: dm_create_subdm.cpp:12
MoFEM::TimeScaleVector::TimeScaleVector
TimeScaleVector(std::string name="-time_vector_file", bool error_if_file_not_given=false)
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
MoFEM::TimeScale::getScale
double getScale(const double time)
Get scaling at a given time.
Definition: ScalingMethod.cpp:133
FTensor::Index< 'i', 3 >
convert.n
n
Definition: convert.py:82
MoFEM::TimeScale::errorIfFileNotGiven
bool errorIfFileNotGiven
Definition: ScalingMethod.hpp:89
MoFEM::TimeScale::defaultDelimiter
static const std::string defaultDelimiter
comma or space
Definition: ScalingMethod.hpp:87
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
MoFEM::TimeScale::fileNameFlag
std::string fileNameFlag
Definition: ScalingMethod.hpp:84
MOFEM_TAG_AND_LOG_C
#define MOFEM_TAG_AND_LOG_C(channel, severity, tag, format,...)
Tag and log in channel.
Definition: LogManager.hpp:370
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::PetscOptionsGetString
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
Definition: DeprecatedPetsc.hpp:172
MoFEM::TimeScaleVector::nAme
string nAme
Definition: ScalingMethod.hpp:114
MoFEM::TimeScale::scalingMethod
std::function< double(double)> scalingMethod
Definition: ScalingMethod.hpp:90
MoFEM::ScalingMethod::getScale
virtual double getScale(const double time)
Get scaling at given time.
Definition: ScalingMethod.cpp:15
dt
double dt
Definition: heat_method.cpp:26
MoFEM::TimeScale::getScaleFromData
double getScaleFromData(const double time)
Get scaling at a given time when the scalar values have been provided. Uses linear interpolation on t...
Definition: ScalingMethod.cpp:105
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
MoFEM::TimeScale::getLinearScale
double getLinearScale(const double time)
Returns the value of time.
Definition: ScalingMethod.cpp:103
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
MOFEM_INVALID_DATA
@ MOFEM_INVALID_DATA
Definition: definitions.h:36
MoFEM::TimeScaleVector::getVector
virtual FTensor::Tensor1< double, SPACE_DIM > getVector(const double time)
Definition: ScalingMethod.cpp:220