Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://www.mooseframework.org
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #pragma once
11 :
12 : #include "KokkosTypes.h"
13 : #include "KokkosVariableValue.h"
14 : #include "KokkosMaterialProperty.h"
15 : #include "KokkosMaterialPropertyValue.h"
16 : #include "KokkosFunction.h"
17 :
18 : #include "ConsoleStreamInterface.h"
19 :
20 : #include "peglib.h"
21 :
22 : namespace Moose::Kokkos
23 : {
24 :
25 : class FunctionBase;
26 :
27 : /**
28 : * Parsing Expression Grammar (PEG)
29 : */
30 : class PEGParser
31 : {
32 : public:
33 : /**
34 : * Constructor
35 : * @param expression The function expression
36 : * @param console The console object
37 : */
38 : PEGParser(const std::string & expression, const ConsoleStream * console = nullptr);
39 :
40 : /**
41 : * Get AST
42 : * @returns The AST
43 : */
44 445 : auto ast() const { return _ast; }
45 : /**
46 : * Get input expression
47 : * @returns The input expression
48 : */
49 8 : const std::string & expression() const { return _expression; }
50 :
51 : private:
52 : /**
53 : * Parser object
54 : */
55 : peg::parser _parser;
56 : /**
57 : * Abstract Syntax Tree (AST)
58 : */
59 : std::shared_ptr<peg::Ast> _ast;
60 : /**
61 : * Input expression
62 : */
63 : const std::string _expression;
64 : };
65 :
66 : /**
67 : * Reverse Polish Notation (RPN) builder
68 : */
69 : class RPNBuilder
70 : {
71 : public:
72 : /**
73 : * Constructor
74 : * @param expression The function expression
75 : * @param console The console object
76 : */
77 471 : RPNBuilder(const std::string & expression, const ConsoleStream * console = nullptr)
78 260 : : _parser(expression, console)
79 : {
80 445 : }
81 :
82 : /**
83 : * RPN opcode
84 : */
85 : enum class Opcode
86 : {
87 : NUM,
88 : VAR,
89 : NEG,
90 : NOT,
91 : // Binary operators
92 : ADD,
93 : SUB,
94 : MUL,
95 : DIV,
96 : AND,
97 : OR,
98 : EQ,
99 : NEQ,
100 : LT,
101 : LEQ,
102 : GT,
103 : GEQ,
104 : // Functions
105 : ABS,
106 : ACOS,
107 : ACOSH,
108 : ASIN,
109 : ASINH,
110 : ATAN,
111 : ATAN2,
112 : ATANH,
113 : CBRT,
114 : CEIL,
115 : COS,
116 : COSH,
117 : COT,
118 : CSC,
119 : EXP,
120 : EXP2,
121 : FLOOR,
122 : HYPOT,
123 : IF,
124 : INT,
125 : LOG,
126 : LOG2,
127 : LOG10,
128 : MAX,
129 : MIN,
130 : POW,
131 : SEC,
132 : SIN,
133 : SINH,
134 : SQRT,
135 : TAN,
136 : TANH,
137 : TRUNC
138 : };
139 :
140 : /**
141 : * RPN instruction
142 : */
143 : struct Instruction
144 : {
145 : // Opcode
146 : Opcode op;
147 : // Original string
148 : std::string text;
149 : // Constant of variable index
150 : unsigned int arg = libMesh::invalid_uint;
151 : };
152 :
153 : /**
154 : * Variable identifier
155 : */
156 : struct Variable
157 : {
158 : // Variable index
159 : unsigned int idx;
160 : // Pointer to a scalar value
161 : const Real * scalar = nullptr;
162 : // Pointer to a field variable
163 : const VariableValue * field = nullptr;
164 : // Pointer to a material property
165 : const MaterialProperty<Real> * property = nullptr;
166 : // Pointer to a function
167 : const Function * function = nullptr;
168 : // Whether the variable was associated
169 2514 : bool associated() const { return scalar || field || property || function; }
170 : };
171 :
172 : /**
173 : * Build RPN from AST
174 : * @param ast The current node
175 : */
176 : ///@{
177 : void build(const peg::Ast & ast);
178 445 : void build() { build(*_parser.ast()); }
179 : ///@}
180 :
181 : /**
182 : * Print RPN sequence for debugging
183 : * @param console The console object
184 : */
185 : void printRPN(const ConsoleStream & console) const;
186 : /**
187 : * Get RPN sequence
188 : * @returns The RPN sequence
189 : */
190 437 : const std::vector<Instruction> & getRPN() const { return _rpn; }
191 :
192 : /**
193 : * Add default variables
194 : */
195 : void addDefaultVariables();
196 : /**
197 : * Get whether default variables were added
198 : * @returns Whether default variables were added
199 : */
200 437 : bool hasDefaultVariables() const { return _has_default_variables; }
201 :
202 : /**
203 : * Add a parsed function constant
204 : * @param number The constant
205 : * @returns The constant index
206 : */
207 : unsigned int addNumber(Real number);
208 : /**
209 : * Add a parsed function variable
210 : * @param name The variable name
211 : * @returns The variable index
212 : */
213 : unsigned int addVariable(const std::string & name);
214 : /**
215 : * Associate a variable with a scalar value
216 : * @param name The variable name
217 : * @param value The pointer to the scalar value
218 : */
219 : void associateScalar(const std::string & name, const Real * scalar);
220 : /**
221 : * Associate a variable with a field variable
222 : * @param name The variable name
223 : * @param variable The pointer to the field variable
224 : */
225 : void associateField(const std::string & name, const VariableValue * field);
226 : /**
227 : * Associate a variable with a material property
228 : * @param name The variable name
229 : * @param variable The pointer to the material property
230 : */
231 : void associateProperty(const std::string & name, const MaterialProperty<Real> * property);
232 : /**
233 : * Associate a variable with a function
234 : * @param name The variable name
235 : * @param function The pointer to the function
236 : */
237 : void associateFunction(const std::string & name, const Function * function);
238 :
239 : /**
240 : * Get numbers used in the expression
241 : * @returns The numbers
242 : */
243 437 : const std::vector<Real> & getNumbers() const { return _numbers; }
244 : /**
245 : * Get variables used in the expression
246 : * @returns The variables
247 : */
248 437 : const std::unordered_map<std::string, Variable> & getVariables() const { return _variables; }
249 :
250 : /**
251 : * Finalize the builder and prevent further changes
252 : */
253 437 : void finalize() { _finalized = true; }
254 : /**
255 : * Get whether the builder was finalized
256 : * @returns Whether the builder was finalized
257 : */
258 437 : bool finalized() const { return _finalized; }
259 :
260 : private:
261 : /**
262 : * Map from unary operators to opcodes
263 : */
264 : static inline const std::map<std::string, Opcode> _unary_opcode_map = {{"-", Opcode::NEG},
265 : {"!", Opcode::NOT}};
266 :
267 : /**
268 : * Map from binary operators to opcodes
269 : */
270 : static inline const std::map<std::string, Opcode> _binary_opcode_map = {{"+", Opcode::ADD},
271 : {"-", Opcode::SUB},
272 : {"*", Opcode::MUL},
273 : {"/", Opcode::DIV},
274 : {"^", Opcode::POW},
275 : {"&", Opcode::AND},
276 : {"|", Opcode::OR},
277 : {"=", Opcode::EQ},
278 : {"!=", Opcode::NEQ},
279 : {"<", Opcode::LT},
280 : {"<=", Opcode::LEQ},
281 : {">", Opcode::GT},
282 : {">=", Opcode::GEQ}};
283 :
284 : /**
285 : * Map from functions to opcodes and the expected number of arguments
286 : */
287 : static inline const std::map<std::string, std::pair<Opcode, unsigned int>> _function_opcode_map =
288 : {{"abs", {Opcode::ABS, 1}}, {"acos", {Opcode::ACOS, 1}}, {"acosh", {Opcode::ACOSH, 1}},
289 : {"asin", {Opcode::ASIN, 1}}, {"asinh", {Opcode::ASINH, 1}}, {"atan", {Opcode::ATAN, 1}},
290 : {"atan2", {Opcode::ATAN2, 2}}, {"atanh", {Opcode::ATANH, 1}}, {"cbrt", {Opcode::CBRT, 1}},
291 : {"ceil", {Opcode::CEIL, 1}}, {"cos", {Opcode::COS, 1}}, {"cosh", {Opcode::COSH, 1}},
292 : {"cot", {Opcode::COT, 1}}, {"csc", {Opcode::CSC, 1}}, {"exp", {Opcode::EXP, 1}},
293 : {"exp2", {Opcode::EXP2, 1}}, {"floor", {Opcode::FLOOR, 1}}, {"hypot", {Opcode::HYPOT, 2}},
294 : {"if", {Opcode::IF, 3}}, {"int", {Opcode::INT, 1}}, {"log", {Opcode::LOG, 1}},
295 : {"log2", {Opcode::LOG2, 1}}, {"log10", {Opcode::LOG10, 1}}, {"max", {Opcode::MAX, 2}},
296 : {"min", {Opcode::MIN, 2}}, {"pow", {Opcode::POW, 2}}, {"sec", {Opcode::SEC, 1}},
297 : {"sin", {Opcode::SIN, 1}}, {"sinh", {Opcode::SINH, 1}}, {"sqrt", {Opcode::SQRT, 1}},
298 : {"tan", {Opcode::TAN, 1}}, {"tanh", {Opcode::TANH, 1}}, {"trunc", {Opcode::TRUNC, 1}}};
299 :
300 : /**
301 : * PEG parser
302 : */
303 : PEGParser _parser;
304 : /**
305 : * RPN sequence
306 : */
307 : std::vector<Instruction> _rpn;
308 : /**
309 : * Numbers used in the function
310 : */
311 : std::vector<Real> _numbers;
312 : /**
313 : * Variables used in the function
314 : */
315 : std::unordered_map<std::string, Variable> _variables;
316 : /**
317 : * Whether default variables were added
318 : */
319 : bool _has_default_variables = false;
320 : /**
321 : * Whether builder was finalized
322 : */
323 : bool _finalized = false;
324 :
325 : /**
326 : * Print a pretty error showing the position of error
327 : * @param ast The erroneous AST node
328 : * @param message The error message
329 : */
330 : [[noreturn]] void builderError(const peg::Ast & ast, const std::string & message) const;
331 :
332 : /**
333 : * Error on attempts to update the builder after finalization
334 : */
335 : void checkFinalized();
336 : };
337 :
338 : /**
339 : * Reverse Polish Notation (RPN) evaluator
340 : */
341 : class RPNEvaluator
342 : {
343 : using Opcode = RPNBuilder::Opcode;
344 :
345 : public:
346 : /**
347 : * Default constructor
348 : */
349 445 : RPNEvaluator() = default;
350 : /**
351 : * Copy constructor for parallel dispatch
352 : */
353 : RPNEvaluator(const RPNEvaluator & evaluator);
354 :
355 : /**
356 : * Initialize RPN evaluator from an RPN builder
357 : * @param builder The RPN builder
358 : */
359 : void init(const RPNBuilder & builder);
360 :
361 : /**
362 : * Evaluate RPN at point (t,x,y,z)
363 : * @param t The time
364 : * @param p The location in space (x,y,z)
365 : * @param qp The local quadrature point index
366 : * @param datum The Datum object of the current thread
367 : * @returns The evaluated value
368 : */
369 : KOKKOS_FUNCTION Real eval(const Real t,
370 : const Real3 p,
371 : const unsigned int qp = 0,
372 : Datum * datum = nullptr) const;
373 :
374 : private:
375 : /**
376 : * RPN instruction
377 : */
378 : struct Instruction
379 : {
380 : // Opcode
381 : Opcode op;
382 : // Constant of variable index
383 : unsigned int arg = libMesh::invalid_uint;
384 : };
385 : /**
386 : * Types of variables
387 : */
388 : enum class VariableType
389 : {
390 : SCALAR,
391 : FIELD,
392 : MATERIAL,
393 : FUNCTION
394 : };
395 :
396 : /**
397 : * RPN sequence
398 : */
399 : Array<Instruction> _rpn;
400 : /**
401 : * Numbers used in the function
402 : */
403 : Array<Real> _numbers;
404 : /**
405 : * Scalar variables used in the function
406 : */
407 : Array<Real> _scalars;
408 : /**
409 : * Field variables used in the function
410 : */
411 : Array<VariableValue> _fields;
412 : /**
413 : * Material properties used in the function
414 : */
415 : Array<MaterialProperty<Real>> _properties;
416 : /**
417 : * Functions used in the function
418 : */
419 : Array<Function> _functions;
420 : /**
421 : * Types of variables and indices
422 : */
423 : Array<::Kokkos::pair<VariableType, unsigned int>> _variables;
424 : /**
425 : * Pointers to the associated quantities of variables
426 : */
427 : std::vector<const void *> _pointers;
428 :
429 : /**
430 : * Number of scalar variables
431 : */
432 : unsigned int _num_scalars = 0;
433 : /**
434 : * Number of field variables
435 : */
436 : unsigned int _num_fields = 0;
437 : /**
438 : * Number of material properties
439 : */
440 : unsigned int _num_properties = 0;
441 : /**
442 : * Number of functions
443 : */
444 : unsigned int _num_functions = 0;
445 :
446 : /**
447 : * Fixed stack size
448 : */
449 : static constexpr unsigned int _stack_size = 10;
450 : /**
451 : * Epsilon for equality comparison
452 : */
453 : static constexpr double _epsilon = 1.0e-12;
454 : };
455 :
456 : #define KOKKOS_FPARSER_COMPARE(code, op, epsilon) \
457 : case Opcode::code: \
458 : stack[head - 2] = stack[head - 2] - stack[head - 1] op epsilon; \
459 : --head; \
460 : break
461 :
462 : #define KOKKOS_FPARSER_BINARY(code, op) \
463 : case Opcode::code: \
464 : stack[head - 2] = stack[head - 2] op stack[head - 1]; \
465 : --head; \
466 : break
467 :
468 : #define KOKKOS_FPARSER_FUNCTION_1(code, func) \
469 : case Opcode::code: \
470 : stack[head - 1] = ::Kokkos::func(stack[head - 1]); \
471 : break
472 :
473 : #define KOKKOS_FPARSER_FUNCTION_2(code, func) \
474 : case Opcode::code: \
475 : stack[head - 2] = ::Kokkos::func(stack[head - 2], stack[head - 1]); \
476 : --head; \
477 : break
478 :
479 : #define KOKKOS_FPARSER_FUNCTION_INV_1(code, func) \
480 : case Opcode::code: \
481 : stack[head - 1] = 1.0 / ::Kokkos::func(stack[head - 1]); \
482 : break
483 :
484 : KOKKOS_FUNCTION inline Real
485 3214432 : RPNEvaluator::eval(const Real t, const Real3 p, const unsigned int qp, Datum * datum) const
486 : {
487 : Real stack[_stack_size];
488 :
489 : // Stack head position
490 3214432 : unsigned int head = 0;
491 :
492 44347731 : for (unsigned int pos = 0; pos < _rpn.size(); ++pos)
493 : {
494 : KOKKOS_ASSERT(head < _stack_size);
495 :
496 41133299 : const auto inst = _rpn[pos];
497 :
498 41133299 : switch (_rpn[pos].op)
499 : {
500 6583916 : case Opcode::NUM:
501 6583916 : stack[head] = _numbers[inst.arg];
502 6583916 : ++head;
503 6583916 : break;
504 16366108 : case Opcode::VAR:
505 16366108 : if (inst.arg < 3)
506 8573322 : stack[head] = p(inst.arg);
507 7792786 : else if (inst.arg == 3)
508 2022200 : stack[head] = t;
509 : else
510 : {
511 5770586 : const auto type = _variables[inst.arg].first;
512 5770586 : const auto idx = _variables[inst.arg].second;
513 :
514 5770586 : switch (type)
515 : {
516 4038506 : case VariableType::SCALAR:
517 4038506 : stack[head] = _scalars[idx];
518 4038506 : break;
519 1121240 : case VariableType::FIELD:
520 : {
521 : KOKKOS_ASSERT(datum);
522 :
523 1121240 : stack[head] = _fields[idx](*datum, qp);
524 1121240 : break;
525 : }
526 144000 : case VariableType::MATERIAL:
527 : {
528 : KOKKOS_ASSERT(datum);
529 :
530 144000 : stack[head] = _properties[idx](*datum, qp);
531 144000 : break;
532 : }
533 466840 : case VariableType::FUNCTION:
534 466840 : stack[head] = _functions[idx].value(t, p);
535 466840 : break;
536 : }
537 : }
538 16366108 : ++head;
539 16366108 : break;
540 7 : case Opcode::NEG:
541 7 : stack[head - 1] = -stack[head - 1];
542 7 : break;
543 4 : case Opcode::NOT:
544 4 : stack[head - 1] = !::Kokkos::round(stack[head - 1]);
545 4 : break;
546 1 : case Opcode::EQ:
547 1 : stack[head - 2] = ::Kokkos::abs(stack[head - 2] - stack[head - 1]) <= _epsilon;
548 1 : --head;
549 1 : break;
550 1 : case Opcode::NEQ:
551 1 : stack[head - 2] = ::Kokkos::abs(stack[head - 2] - stack[head - 1]) > _epsilon;
552 1 : --head;
553 1 : break;
554 8 : case Opcode::AND:
555 8 : stack[head - 2] = ::Kokkos::round(stack[head - 2]) && ::Kokkos::round(stack[head - 1]);
556 8 : --head;
557 8 : break;
558 5 : case Opcode::OR:
559 5 : stack[head - 2] = ::Kokkos::round(stack[head - 2]) || ::Kokkos::round(stack[head - 1]);
560 5 : --head;
561 5 : break;
562 2019203 : case Opcode::IF:
563 2019203 : stack[head - 3] = ::Kokkos::round(stack[head - 3]) ? stack[head - 2] : stack[head - 1];
564 2019203 : head -= 2;
565 2019203 : break;
566 10 : KOKKOS_FPARSER_COMPARE(LT, <, -_epsilon);
567 2 : KOKKOS_FPARSER_COMPARE(LEQ, <=, _epsilon);
568 2019202 : KOKKOS_FPARSER_COMPARE(GT, >, _epsilon);
569 2 : KOKKOS_FPARSER_COMPARE(GEQ, >=, -_epsilon);
570 7020983 : KOKKOS_FPARSER_BINARY(ADD, +);
571 560622 : KOKKOS_FPARSER_BINARY(SUB, -);
572 5068869 : KOKKOS_FPARSER_BINARY(MUL, *);
573 6 : KOKKOS_FPARSER_BINARY(DIV, /);
574 2 : KOKKOS_FPARSER_FUNCTION_1(ABS, abs);
575 1 : KOKKOS_FPARSER_FUNCTION_1(ACOS, acos);
576 1 : KOKKOS_FPARSER_FUNCTION_1(ACOSH, acosh);
577 1 : KOKKOS_FPARSER_FUNCTION_1(ASIN, asin);
578 1 : KOKKOS_FPARSER_FUNCTION_1(ASINH, asinh);
579 1 : KOKKOS_FPARSER_FUNCTION_1(ATAN, atan);
580 1 : KOKKOS_FPARSER_FUNCTION_2(ATAN2, atan2);
581 1 : KOKKOS_FPARSER_FUNCTION_1(ATANH, atanh);
582 1 : KOKKOS_FPARSER_FUNCTION_1(CBRT, cbrt);
583 1 : KOKKOS_FPARSER_FUNCTION_1(CEIL, ceil);
584 3 : KOKKOS_FPARSER_FUNCTION_1(COS, cos);
585 1 : KOKKOS_FPARSER_FUNCTION_1(COSH, cosh);
586 1 : KOKKOS_FPARSER_FUNCTION_INV_1(COT, tan);
587 1 : KOKKOS_FPARSER_FUNCTION_INV_1(CSC, sin);
588 1 : KOKKOS_FPARSER_FUNCTION_1(EXP, exp);
589 1 : KOKKOS_FPARSER_FUNCTION_1(EXP2, exp2);
590 1 : KOKKOS_FPARSER_FUNCTION_1(FLOOR, floor);
591 1 : KOKKOS_FPARSER_FUNCTION_2(HYPOT, hypot);
592 2 : KOKKOS_FPARSER_FUNCTION_1(INT, round);
593 1 : KOKKOS_FPARSER_FUNCTION_1(LOG, log);
594 1 : KOKKOS_FPARSER_FUNCTION_1(LOG2, log2);
595 1 : KOKKOS_FPARSER_FUNCTION_1(LOG10, log10);
596 2 : KOKKOS_FPARSER_FUNCTION_2(MAX, max);
597 1 : KOKKOS_FPARSER_FUNCTION_2(MIN, min);
598 1027470 : KOKKOS_FPARSER_FUNCTION_2(POW, pow);
599 1 : KOKKOS_FPARSER_FUNCTION_INV_1(SEC, cos);
600 466844 : KOKKOS_FPARSER_FUNCTION_1(SIN, sin);
601 1 : KOKKOS_FPARSER_FUNCTION_1(SINH, sinh);
602 1 : KOKKOS_FPARSER_FUNCTION_1(SQRT, sqrt);
603 2 : KOKKOS_FPARSER_FUNCTION_1(TAN, tan);
604 1 : KOKKOS_FPARSER_FUNCTION_1(TANH, tanh);
605 1 : KOKKOS_FPARSER_FUNCTION_1(TRUNC, trunc);
606 : }
607 : }
608 :
609 : KOKKOS_ASSERT(head == 1);
610 :
611 3214432 : return stack[head - 1];
612 : }
613 :
614 : } // namespace Moose::Kokkos
|
