Line data Source code
1 : //
2 : // This is a general-purpose routine for specifying time-interpolated
3 : // quantities in an input file.
4 : // Interpolators obey the following syntax:
5 : //
6 : // .. code-block:: bash
7 : //
8 : // (t1,t2,..,tn;v1,v2,...,vn)
9 : //
10 : // where :code:`tn` are times and :code:`vn` are values at those times.
11 : //
12 : // :bdg-danger-line:`Note: do not include whitespace in an interpolator, this will mess with the parser.`
13 : //
14 : // Interpolators can usually be used in place of regular numbers, but only
15 : // when supported by the calling parse function.
16 : //
17 : // For instance, the `BC::Operator::Elastic::TensionTest` allows the user to specify
18 : // fixed displacement; for instance, to specify a value of 0.1:
19 : //
20 : // .. code-block:: bash
21 : //
22 : // bc.tenstion_test.disp = 0.1
23 : //
24 : // However, it may be desirable to change the value over time.
25 : // In this case, one may specify an interpolator string instead:
26 : //
27 : // .. code-block:: bash
28 : //
29 : // bc.tenstion_test.disp = (1.0,1.5,2.0;0.0,0.1,0.0)
30 : //
31 : // This will cause the displacement value to be zero for :code:`t<1.0`;
32 : // ramp to :code:`0.1` as :code:`1.0<t<2.0`;
33 : // ramp back to :code:`0.0` as :code:`2.0<t<3.0`;
34 : // and then remain at the constant value of :code:`3.0` for :code:`t>3.0`.
35 : //
36 :
37 : #ifndef NUMERIC_INTERPOLATOR_LINEAR_H_
38 : #define NUMERIC_INTERPOLATOR_LINEAR_H_
39 :
40 : #include <AMReX.H>
41 : #include <AMReX_MultiFab.H>
42 : #include "Set/Set.H"
43 : #include "Util/Util.H"
44 : #include "Numeric/Interpolator/Interpolator.H"
45 : #include "IO/ParmParse.H"
46 :
47 : #include "Unit/Unit.H"
48 :
49 : namespace Numeric
50 : {
51 : namespace Interpolator
52 : {
53 : template <class T>
54 : class Linear : public Interpolator<T>
55 : {
56 : public:
57 466 : Linear() {};
58 1546 : ~Linear() {};
59 530 : Linear(const Set::Scalar &a_val)
60 530 : {
61 530 : interval_points.push_back(0.0);
62 530 : data_points.push_back(a_val);
63 530 : }
64 222 : Linear(const std::string a_str, Unit unit_time = Unit::Less(), Unit unit_data = Unit::Less())
65 222 : {
66 222 : define(a_str, unit_time, unit_data);
67 222 : }
68 2 : Linear(const std::vector<T> _data_points, const std::vector<Set::Scalar> _interval_points)
69 2 : {
70 2 : define(_data_points,_interval_points);
71 2 : };
72 :
73 224 : void define(const std::string a_str, Unit a_unit_time, Unit a_unit_val)
74 : {
75 224 : std::string str = a_str;
76 :
77 2588 : if (Util::String::Contains(str,"(") || Util::String::Contains(str,":") || Util::String::Contains(str,")") || Util::String::Contains(str,","))
78 : {
79 108 : if (!Util::String::Contains(str,"(")) Util::Abort(INFO,"Mismatched parentheses while trying to parse ",str);
80 108 : if (!Util::String::Contains(str,")")) Util::Abort(INFO,"Mismatched parentheses while trying to parse ",str);
81 108 : if (!Util::String::Contains(str,":")) Util::Abort(INFO,"Missing : in interpolator while trying to parse ",str);
82 180 : Util::String::ReplaceAll(str,"(","");
83 180 : Util::String::ReplaceAll(str,")","");
84 :
85 : std::vector<std::string> splitstr = Util::String::Split(str,':');
86 36 : std::string str_time = splitstr[0];
87 36 : std::string str_val = splitstr[1];
88 :
89 : std::vector<std::string> str_time_arr = Util::String::Split(str_time,',');
90 : std::vector<std::string> str_val_arr = Util::String::Split(str_val,',');
91 :
92 36 : if (str_time_arr.size() != str_val_arr.size()) Util::Abort(INFO,"Mismatched number of time values vs array values while trying to parse ", str);
93 :
94 108 : for (unsigned int i = 0; i < str_time_arr.size(); i++)
95 : {
96 72 : Unit unit_time = Unit::Parse(str_time_arr[i]);
97 72 : if (!unit_time.isType(a_unit_time) && !unit_time.isType(Unit::Less()))
98 0 : Util::Exception(INFO,"Invalid unit, ", str_time_arr[i]);
99 72 : Unit unit_val = Unit::Parse(str_val_arr[i]);
100 72 : if (!unit_val.isType(a_unit_val) && !unit_val.isType(Unit::Less()))
101 0 : Util::Exception(INFO,"Invalid unit, ", str_val_arr[i]);
102 72 : interval_points.push_back(unit_time.normalized_value());
103 72 : if (i > 0) if (interval_points[i] < interval_points[i-1]) Util::Abort(INFO,"Time series must monotonically increase - error while trying to parse ", str);
104 72 : data_points.push_back(unit_val.normalized_value());
105 : }
106 36 : }
107 : else
108 : {
109 188 : interval_points.push_back(0.0);
110 188 : data_points.push_back(std::stod(a_str));
111 : }
112 224 : }
113 :
114 2 : void define(const std::vector<T> _data_points, const std::vector<Set::Scalar> _interval_points)
115 : {
116 2 : data_points = _data_points;
117 2 : interval_points = _interval_points;
118 2 : if(data_points.size() != interval_points.size())
119 0 : Util::Abort(INFO,"Data points and interval points have different sizes");
120 2 : };
121 :
122 94796850 : T operator() (const Set::Scalar point) const
123 : {
124 : // If there is only one point stored, return it.
125 94796850 : if (data_points.size() == 1) return data_points[0];
126 :
127 : // Do this if point is below interpolation range
128 9348 : if(point < interval_points[0])
129 : {
130 200 : return data_points[0];
131 : }
132 :
133 9148 : Set::Scalar interval = interval_points[1]-interval_points[0];
134 9148 : int start = 0;
135 9148 : if (uniform)
136 : {
137 9148 : start = (int)( (point - interval_points[0]) / interval);
138 9148 : start -= 1;
139 9148 : start = std::max(start,0);
140 9148 : start = std::min(start,(int)(interval_points.size()-2));
141 : }
142 :
143 : // Do this if point is within interpolation range
144 9544 : for (unsigned int i = start; i<interval_points.size()-1; i++)
145 : {
146 9346 : if (interval_points[i] <= point && point <= interval_points[i+1])
147 : {
148 8950 : return data_points[i] +
149 8950 : (point - interval_points[i]) * (data_points[i+1] - data_points[i]) /
150 8950 : (interval_points[i+1] - interval_points[i]);
151 :
152 : }
153 : }
154 :
155 : // Otherwise point must be outside interpolation range, so
156 : // return the value for the highest point
157 198 : return data_points[interval_points.size()-1];
158 : }
159 :
160 2 : static void Parse(Linear<T> & value, IO::ParmParse & pp, Unit unit_time = Unit::Less(), Unit unit_value = Unit::Less())
161 : {
162 2 : std::string in;
163 2 : pp_query("str",in); // Interpolator string used when Parsed from queryclass.
164 2 : value.define(in, unit_time, unit_value);
165 2 : }
166 :
167 :
168 : static Linear<T> Read(std::string filename, int derivative = 0);
169 :
170 : protected:
171 : std::vector<T> data_points;
172 : std::vector<Set::Scalar> interval_points;
173 : bool uniform = true;
174 :
175 : public:
176 : template<class U> friend Linear<U> operator * (const Set::Scalar alpha, const Linear<U> &b);
177 : };
178 :
179 : template<class T>
180 : AMREX_FORCE_INLINE
181 : Linear<T> operator * (const Set::Scalar alpha, const Linear<T> &b)
182 : {
183 : Linear<T> ret;
184 : for (unsigned int i = 0; i < b.data_points.size(); i++)
185 : ret.data_points.push_back(alpha * b.data_points[i]);
186 : ret.interval_points = b.interval_points;
187 : ret.uniform = b.uniform;
188 : return ret;
189 : }
190 :
191 : template<>
192 : ALAMO_SINGLE_DEFINITION
193 0 : Linear<Set::Scalar> Linear<Set::Scalar>::Read(std::string filename, int derivative)
194 : {
195 0 : std::ifstream input;
196 0 : input.open(filename);
197 0 : std::string line;
198 0 : std::vector<Set::Scalar> theta, w;
199 : // parse input
200 0 : while (std::getline(input, line))
201 : { std::vector<std::string> dat = Util::String::Split(line);
202 0 : theta.push_back(std::stof(dat[0]));
203 0 : w.push_back(std::stof(dat[1]));
204 0 : }
205 0 : if (derivative == 0)
206 : {
207 0 : return Linear<Set::Scalar>(w,theta);
208 : }
209 0 : else if (derivative == 1)
210 : {
211 0 : std::vector<Set::Scalar> thetasmall, dw;
212 0 : for (unsigned int i = 1; i < theta.size() - 1; i++)
213 : {
214 0 : thetasmall.push_back(theta[i]);
215 0 : dw.push_back((w[i + 1] - w[i - 1]) / (theta[i + 1] - theta[i - 1]));
216 : }
217 0 : return Linear<Set::Scalar>(dw,thetasmall);
218 0 : }
219 0 : else if (derivative == 2)
220 : {
221 0 : std::vector<Set::Scalar> thetasmall, ddw;
222 0 : for (unsigned int i = 1; i < theta.size() - 1; i++)
223 : {
224 0 : thetasmall.push_back(theta[i]);
225 0 : ddw.push_back((w[i + 1] - 2.0 * w[i] + w[i - 1]) / ((theta[i + 1] - theta[i]) * (theta[i] - theta[i - 1])));
226 : }
227 0 : return Linear<Set::Scalar>(ddw,thetasmall);
228 0 : }
229 0 : else Util::Abort(INFO,"derivative must be 0,1,2 but got ",derivative);
230 0 : exit(-1);
231 0 : }
232 :
233 :
234 : }
235 : }
236 :
237 : #endif
|
