diff --git a/README.md b/README.md index 6b49fc9..a52ae79 100644 --- a/README.md +++ b/README.md @@ -1,7 +1,12 @@ -# Finite Element models 4 Intrinsic Nonlinear Aeroelastics in JAX [FENIAX] +# Finite Element models 4 Nonlinear Intrinsic Aeroelastics in JAX [FENIAX] FENIAX is an aeroelastic toolbox written and parallelized in Python, which acts as a post-processor of commercial software such as MSC Nastran. Arbitrary FE models built for linear aeroelastic analysis are enhanced with geometric nonlinear effects, flight dynamics and linearized state-space solutions about nonlinear equilibrium. +Some of the key features of the software are: +- Leveraging on the numerical library JAX and optimised algorithms, a high performance is achieved that leads to simulation times comparable to the linear counterparts on conventional platforms. +- The software runs on modern hardware architectures such as GPUs in a addition to standard CPUs. +- Algorithm differentiation (AD) of the aeroelastic response is available via JAX primitives. +- Concurrent simulations for multiple load cases are being developed. ## Installation @@ -17,30 +22,56 @@ pip install . pip install -e .[all] ``` -- see setup.py file for the options available. Python 3.9+ is required but 3.11+ is recommended. +- see pyproject.toml file for the options available. Python 3.10+ is required. + +- To install with GPU support install jax first: +``` +pip install -U "jax[cuda12]" +pip install -e ".[all]" +``` + + +### Python environment +Although it is not necessary, it is recommended that the software is installed in its own environment. Options that have been tested out follow. + +- _Conda_: + Although it is not necessary, If conda is being used as package manager, one can make a specific environment as, ``` conda create -n feniax python=3.11 -conda activate fem4inas +conda activate feniax ``` -- If pytest has been installed, check everything is OK by running the tests: +If pytest has been installed, check everything is OK by running the tests: ``` -pytest tests +pytest ``` - Thus a typical installation would comprise of these 4 steps: ``` conda create -n feniax.python=3.11 -conda activate fem4inas +conda activate feniax pip install -e .[all] -pytest tests +pytest ``` +- _pyenv_: Navigate to the root directory and run the following: + +``` + pyenv install 3.11.10 + pyenv virtualenv 3.11.10 feniax + pyenv local feniax + pip install -e .[all] + pytest +``` +By setting pyenv local to feniax, every time one moves to feniax directory the environment is automatically activated + + ## Documentation Available at https://acea15.github.io/FENIAX/ + ## Simulation Examples The most relevant examples in the code base are shown here, these and more can be found in the folder `/examples` @@ -77,17 +108,4 @@ They are also part of a large test suite that is integrated into the development #### 3D dynamics ![Free flying structure 3D](./docs/media/SimoFFB3D_optimized.gif) -### Industrial Aircraft model -!!! success - Linear response validated with MSC Nastran linear aeroelastic solution (sol 146) - -#### Gust clamped model - -[Notebook](./docs/documentation/examples/XRF1/xrf1_nb.md) - - -![XRF1-gustclamped](./docs/media/xrf1_gust_optimized.gif) - -#### Gust trimmed flight -![XRF1-Trim+gust](./docs/media/xrf1_trimgust_optimized.gif) diff --git a/docs/reports/scitech25/figs_ext/BUGmodel.png b/docs/reports/scitech25/figs_ext/BUGmodel.png new file mode 100644 index 0000000..97ed366 Binary files /dev/null and b/docs/reports/scitech25/figs_ext/BUGmodel.png differ diff --git a/docs/reports/scitech25/figs_ext/BUGmodel2.png b/docs/reports/scitech25/figs_ext/BUGmodel2.png new file mode 100644 index 0000000..2fdac39 Binary files /dev/null and b/docs/reports/scitech25/figs_ext/BUGmodel2.png differ diff --git a/docs/reports/scitech25/figs_ext/DiscreteL0.png b/docs/reports/scitech25/figs_ext/DiscreteL0.png new file mode 100644 index 0000000..49b376d Binary files /dev/null and b/docs/reports/scitech25/figs_ext/DiscreteL0.png differ diff --git a/docs/reports/scitech25/figs_ext/DiscreteL2.png b/docs/reports/scitech25/figs_ext/DiscreteL2.png new file mode 100644 index 0000000..7e4d0a1 Binary files /dev/null and b/docs/reports/scitech25/figs_ext/DiscreteL2.png differ diff --git a/docs/reports/scitech25/figs_ext/DiscreteL4.pdf b/docs/reports/scitech25/figs_ext/DiscreteL4.pdf new file mode 100644 index 0000000..cec4471 Binary files /dev/null and b/docs/reports/scitech25/figs_ext/DiscreteL4.pdf differ diff --git a/docs/reports/scitech25/figs_ext/MC1.png b/docs/reports/scitech25/figs_ext/MC1.png new file mode 100644 index 0000000..89bc940 Binary files /dev/null and b/docs/reports/scitech25/figs_ext/MC1.png differ diff --git a/docs/reports/scitech25/figs_ext/bug_gust3d.pdf b/docs/reports/scitech25/figs_ext/bug_gust3d.pdf new file mode 100644 index 0000000..f52a67e Binary files /dev/null and b/docs/reports/scitech25/figs_ext/bug_gust3d.pdf differ diff --git a/docs/reports/scitech25/figs_ext/bug_model2.pdf b/docs/reports/scitech25/figs_ext/bug_model2.pdf new file mode 100644 index 0000000..3c87e8c Binary files /dev/null and b/docs/reports/scitech25/figs_ext/bug_model2.pdf differ diff --git a/docs/reports/scitech25/figs_ext/bug_model3.pdf b/docs/reports/scitech25/figs_ext/bug_model3.pdf new file mode 100644 index 0000000..feda018 Binary files /dev/null and b/docs/reports/scitech25/figs_ext/bug_model3.pdf differ diff --git a/docs/reports/scitech25/figs_ext/bug_model7.pdf b/docs/reports/scitech25/figs_ext/bug_model7.pdf new file mode 100644 index 0000000..d185478 Binary files /dev/null and b/docs/reports/scitech25/figs_ext/bug_model7.pdf differ diff --git a/docs/reports/scitech25/figs_ext/monoeuvre3D.pdf b/docs/reports/scitech25/figs_ext/monoeuvre3D.pdf new file mode 100644 index 0000000..10d08ab Binary files /dev/null and b/docs/reports/scitech25/figs_ext/monoeuvre3D.pdf differ diff --git a/docs/reports/scitech25/main.org b/docs/reports/scitech25/main.org index 918b68b..cc6a7e1 100644 --- a/docs/reports/scitech25/main.org +++ b/docs/reports/scitech25/main.org @@ -1,11 +1,14 @@ -#+TITLE: Aircraft Nonlinear Dynamic Loads at Large Scale Using an Accelerator-Based Distributed Solution -#+TITLE: Massively parallel Aircraft nonlinear dynamic loads using an Accelerator-Based Distributed Solution +#+TITLE: Geometrically Nonlinear Aircraft Loads at Large Scale Using an Accelerator-Based Concurrent Solution +#+AUTHOR: Alvaro Cea\footnote{Research Associate, CAGB 308, South Kensington Campus. (alvaro.cea-esteban15@imperial.ac.uk)}, Rafael Palacios\footnote{Professor in Computational Aeroelasticity, CAGB 338, South Kensington Campus. AIAA Associate Fellow (r.palacios@imperial.ac.uk)} and Lucian Iorga\footnote{Wing Airframe Integrator} + +# \author{Alvaro Cea\footnote{Research Associate, Department of Aeronautics, CAGB 308, South Kensington Campus. (alvaro.cea-esteban15@imperial.ac.uk)}} +# \author{Rafael Palacios\footnote{Professor in Computational Aeroelasticity, Department of Aeronautics and Brahmal Vasudevan Institute for Sustainable Aviation, CAGB 310, South Kensington Campus. AIAA Associate Fellow (r.palacios@imperial.ac.uk)}} +# \affil{Imperial College London, SW7 2AZ, United Kingdom} + +# \author{Lucian Iorga\footnote{Wing Airframe Integrator}} +# \affil{Airbus Operations Ltd., Filton, BS99 7AR, United Kingdom} -# #+TITLE: Parallelized Aeroelastic Solution for Large Scale Simulation of Nonlinear Dynamic Loads on Accelerators -# #+AUTHOR: Alvaro Cea and Rafael Palacios -#+AUTHOR: Alvaro Cea\footnote{Research Associate, CAGB 308, South Kensington Campus. (alvaro.cea-esteban15@imperial.ac.uk)}, and Rafael Palacios\footnote{Professor in Computational Aeroelasticity, CAGB 338, South Kensington Campus. AIAA Associate Fellow (r.palacios@imperial.ac.uk)} #+DATE: -#+BIBLIOGRAPHY:/home/acea/Documents/Engineering.bib :LATEX_PROPERTIES: #+OPTIONS: toc:nil #+OPTIONS: broken-links:mark @@ -26,30 +29,18 @@ #+LATEX_HEADER: \usepackage{amsfonts} #+LATEX_HEADER: \usepackage{enumitem} #+LATEX_HEADER: \usepackage{mathalpha} -#+LATEX_HEADER: \renewcommand{\figurename}{\bf \small Figure} -#+LATEX_HEADER: \renewcommand{\tablename}{\bf \small Table} -#+LATEX_HEADER: \newcommand{\de}{\delta} -#+LATEX_HEADER: \newcommand{\ve}{\text{v}} -#+LATEX_HEADER: \newcommand{\lo}{\mathcal{L}} -#+LATEX_HEADER: \newcommand{\vt}{\overline{\delta\bm{\theta}}} -#+LATEX_HEADER: \newcommand{\vu}{\overline{\delta\bm{u}}} -#+LATEX_HEADER: \newcommand{\e}{\bm{\mathfrak{e}}} -#+LATEX_HEADER: \newcommand{\E}{\bm{\mathbb{E}}} -#+LATEX_HEADER: \newcommand{\T}{\bm{\mathcal{T}}} -#+LATEX_HEADER: \newcommand{\fra}{(\mathtt{1})} -#+LATEX_HEADER: \newcommand{\frb}{(\mathtt{2})} -#+LATEX_HEADER: \newcommand{\fri}{(\mathfrak{i})} -#+LATEX_HEADER: \newcommand{\bs}[1]{\boldsymbol{#1}} -#+LATEX_HEADER: \newcommand{\rhoinf}{\rho} -#+LATEX_HEADER: \newcommand{\Vinf}{U} -#+LATEX_HEADER: \newcommand{\Cl}[1]{c_{l_{#1}}} -#+LATEX_HEADER: \newcommand{\barCl}[1]{\bar{c}_{l_{#1}}} -#+LATEX_HEADER: \newcommand{\Cm}[1]{c_{m_{#1}}} -#+LATEX_HEADER: \newcommand{\barCm}[1]{\bar{c}_{m_{#1}}} -#+LATEX_HEADER: \newcommand{\AIC}{\bs{\mathcal{A}}} +#+LATEX_HEADER: \usepackage{setspace} +#+LATEX_HEADER: \onehalfspacing +# % or: +# \doublespacing :END: +#+begin_abstract +We have extended our newly developed framework for geometrically nonlinear aeroelastic analysis with parallelisation capabilities of multiple load cases. The software can be deployed in modern hardware architectures and the parallelisation has been tested both on CPUs and GPUs. This allows us to perform uncertainty quantification studies and build load envelopes of static and dynamic cases in computational times of the order of seconds. +We explore a high-aspect-ratio-wing aircraft configuration, showing the ability of the solvers to capture nonlinear effects for loadings inducing large deformations and also the importance of these effects when studying modern aircraft concepts. Discrete, manoeuvre and gust loads are assessed for multiple flow and parametric conditions running concurrently. +#+end_abstract + * House keeping :noexport: #+begin_src elisp :results none :tangle no :exports none @@ -59,10 +50,14 @@ '("se" . "src elisp")) (setq org-confirm-babel-evaluate nil) (define-key org-mode-map (kbd "C-c ]") 'org-ref-insert-link) + ;(setq org-latex-pdf-process + ; '("latexmk -pdflatex='pdflatex --syntex=1 -interaction nonstopmode' -pdf -bibtex -f %f")) + ; (setq org-latex-pdf-process (list "latexmk -f -pdf -interaction=nonstopmode -output-directory=%o %f")) (setq org-latex-pdf-process '("latexmk -pdflatex='pdflatex --syntex=1 -interaction nonstopmode' -pdf -bibtex -f %f")) - ;; (setq org-latex-pdf-process (list "latexmk -f -pdf -interaction=nonstopmode -output-directory=%o %f")) - (pyvenv-workon "fem4inas") + ;; (setq org-latex-pdf-process (list "latexmk -f -pdf -interaction=nonstopmode output-directory=%o %f")) + + (pyvenv-workon "feniax") (require 'org-tempo) ;; Veval_blocks -> eval blocks of latex ;; Veval_blocks_run -> eval blocks to obtain results @@ -166,67 +161,34 @@ ** Helper functions * Introduction -Aeroelastic analysis are expected to become critical in the very early phases of the wing design process: while the field was more important in post-design stages to ensure in-flight integrity, it now becomes paramount to capture the cross-couplings between disciplines. +Aeroelastic analysis are expected to become critical in the very early phases of the wing design process: while the field was more important in post-design stages to ensure in-flight integrity, it now becomes paramount to capture the cross-couplings between disciplines. As highlighted in cite:&LIVNE2018, formulations that include nonlinear effects should be developed that not only enhance current modelling techniques but that also allow rapid data turnaround for the industry. Real-time, hardware-in-the-loop flight simulators would also benefit of actively controlled, deformable airplane models. This leads to a more nonlinear landscape, where the overall aerodynamic performance needs to be calculated around a flight shape with large deformations cite:&GRAY2021; the input for efficient control laws account for the steady state and nonlinear couplings cite:&Artola2021; and the loads ultimately sizing the wings are atmospheric disturbances computed in the time-domain cite:&CESNIK2014a. This is also the case for more radical configurations that may or may not exhibit high flexibility but whose aeroelastic behaviour is more uncertain. A more holistic approach to the design also increases the complexity of the processes exponentially, and the trade-offs and cost-benefit analysis may not be possible until robust computational tools are in-place to simulate the different assumptions. Certification of new air vehicles is another important aspect that requires 100,000s of load cases simulations cite:&Kier2017, as it considers manoeuvres and gust loads at different velocities and altitudes, and for a range of mass cases and configurations. This poses another challenge for new methods that aim to include new physics since they normally incur in prohibitively expensive computational times. -Lastly, the mathematical representation of the airframe, embodied in the complex Finite-Element Models (FEMs) built by organizations, encompasses a level of knowledge that is to be preserved when including the new physics mentioned above. -Those previous considerations set the goals for our previous work [[cite:&CEA2023;&CEA2024]]: 1) to be able to perform geometrically nonlinear aeroelastic analysis, 2) to work with generic FEMs in a non-intrusive manner, and 3) to achieve a computational efficiency that is equivalent to present linear methods (if not faster). +Lastly, the mathematical representation of the airframe, embodied in the complex Finite-Element Models (FEMs) built by organizations, encompasses a level of knowledge that is to be preserved when including the new physics mentioned above. \\ +Those previous considerations set the goals for our previous work [[cite:&CEA2023;&CEA2024]]: 1) to be able to perform geometrically nonlinear aeroelastic analysis, 2) to work with generic FEMs in a non-intrusive manner, and 3) to achieve a computational efficiency that is equivalent to present linear methods (if not faster). In this work we explore the latest advances on accelerator's parallelisation and how to integrate them into our solution process to enable large scale aeroelastic simulations under geometrically nonlinear assumptions. Specifically, we set out to characterise the dynamics of highly flexible aircraft in response to the very large envelopes of in-flight loads encountered in the certification process. -In our latest developments we have leveraged the numerical library JAX cite:&jax2018github to build a new simulation environment for time-domain nonlinear aeroelastic analysis that achieves two orders of magnitude speed-ups with respect to standard implementations [[cite:&CEA2024]], is suitable for modern hardware architectures such as GPUs [[cite:&ALVAROCEA2024]], and that is also capable of computing derivatives via algorithmic differentiation [[cite:&ALVAROCEA2024a]]. The power of JAX for scientific computation has been proved recently in fluid dynamics cite:&BEZGIN2023 and solid mechanics cite:&XUE2023 applications. +In our latest developments we have leveraged the numerical library JAX cite:&jax2018github to build a new simulation environment for time-domain nonlinear aeroelastic analysis that achieves two orders of magnitude speed-ups with respect to standard implementations [[cite:&CEA2024]], is suitable for modern hardware architectures such as GPUs [[cite:&ALVAROCEA2024]], and that is also capable of computing derivatives via algorithmic differentiation [[cite:&ALVAROCEA2024a]]. The strength and suitability of JAX for scientific computation has been proved recently in fluid dynamics cite:&BEZGIN2023 and solid mechanics cite:&XUE2023 applications. We want to go one step further by adding parallelisation and distributed computational capabilities to the codebase to tackle the very demanding task of calculating load envelopes while introducing new physics to account for the large displacements and rotations ultra-high-aspect-ratio wings undergo. In this multi-process environment, a Single Program Multiple Data (SPMD) paradigm is employed with the main computation spanning as many devices as available in the cluster and performing collective operations to communicate between devices. -By addressing in one program a substantial part of scenarios during flight (manoeuvres and gust responses at different velocities and altitudes, and for a range of mass cases and configurations), we will be able to produce the critical loading characteristics of the aircraft at a fraction of time. Moreover, the boundaries of these critical cases will be differentiated using the AD, providing gradients for optimization studies as well as additional insights to the designer. -A representative configuration of a high-aspect-ratio aircraft is going to be studied and verification against MSC Nastran will be shown for the linear cases. Afterwards an assessment of geometrically nonlinear effects for trim, manoeuvre and gust loads will be carried out. -This application of modern hardware architectures to aircraft nonlinear load analysis is novel and could potentially be introduced inside current industrial processes. - -* Theory and implementation -In this section we briefly describe the backbone theory of the proposed methods for nonlinear aeroelastic modelling as continuation of the work in [[cite:&CEA2021b;&CEA2023]]. A summary of the main formulation and its integration into an aeroelastic framework are presented next, along with some implementation details. -We start with a global FE model of the airframe as illustrated in Fig. [[workflow]]. - -#+NAME: workflow -#+CAPTION: Workflow of the solution process -#+ATTR_LATEX: :width 1.\textwidth -[[./figs_ext/workflowAIAA3.pdf]] -It is common practice for large-scale aeroelastic models to feature lumped masses along a load path axis that are attached to their corresponding cross-sectional nodes via interpolation elements. -With those characteristics a reduced model can be obtained from a static or dynamic condensation that captures well the stiffness and inertia properties in the condensed matrices, $\pmb{K}_a$ and $\pmb{M}_a$ (Step 1 in Fig. [[workflow]]). The eigenvalue solution of the FEM yields the modal shapes, $\pmb \Phi_0$, and frequencies $\pmb \omega$ (Step 2, however, $\pmb \Phi_0$ is defined on the master nodes and the figure shows the full reconstructed modal shapes). The dynamics of this reduced model are described by a system on nonlinear equations [[cite:&HODGES2003]] written in material velocities, $\bm x_1$, and stresses, $\bm x_2$, as state variables. A modal expansion of those is a key step in seamlessly mapping the global FEM into the nonlinear description. The intrinsic modes are introduced and the projection of the state variables is such $\pmb{x}_1 = \pmb{\Phi}_1\pmb{q}_1$ and $\pmb{x}_2 = \pmb{\Phi}_2\pmb{q}_2$. -A resulting set of four intrinsic modal shapes are directly linked to the displacement modal shapes coming from the global FEM: - -1. Velocity modes, $\bm \Phi_1 = \bm \Phi_0$, which follow after the linear relation with displacements: $\bm x_1 = \dot{\bm x}_0$, $\bm \Phi_1 \bm q_1 = \bm \Phi_0 \dot{\bm q}_0$. +By addressing in one program a substantial part of scenarios during flight (manoeuvres and gust responses at different velocities and altitudes, and for a range of mass cases and configurations), we will be able to produce the critical loading characteristics of the aircraft at a fraction of time. Moreover, as future work we aim to differentiate the boundaries of these critical cases using the already demonstrated capabilities AD, thereby providing gradients for optimization studies as well as additional insights to the designer. +\\ +The paper is organised as follows: Sec. [[Theoretical and computational background]] gives and overview of the theoretical and computational developments that underpin this work with a focus on the new parallelisation capabilities. In sec. [[Results]], a representative configuration of an ultra-high-aspect-ratio aircraft is studied under various loading scenarios that have been parallelised; namely structural static loads, manoeuvre cases for varying flow conditions and dynamic loads with multiple gusts running concurrently. This application of modern hardware architectures to aircraft nonlinear load analysis is novel and could potentially be introduced inside current industrial processes. We conclude in Sec. [[Conclusions]] with a summary of the main advances and the future work that is needed to finalise a formulation that may run in parallel on modern hardware architectures as well as being differentiated. +* Theoretical and computational background +The main aspects of the aeroelastic framework we have developed are presented in this section. +The approach is built on a non-intrusive reduction order process combined with a nonlinear description of the dominant dimension for slender structures. It achieves a nonlinear representation of aeroelastic models of arbitrary complexity in a very efficient manner and without losing the characteristics of the linear model. We target the calculation of flight loads herein, but it can also be applied to the computation of aeroelastic stability phenomena such as flutter or divergence [[cite:&CEA2023]] and to broader multidisciplinary design optimisation problems, which are currently being explored. +The key features of the formulation are: -2. Momentum modes, $\bm \Psi_1 = \bm M_a \bm \Phi_0$. Note from this definition that, for arbitrary distributed mass models, the dynamic condensation technique will produce a fully-populated mass matrix, and the various couplings will be captured after the matrix multiplication. +- Geometrically nonlinear aeroelastic analysis using complex GFEMs: achieved via a three step process in which a condensed model is first produced, the dynamics of this reduced model are described by a system on nonlinear equations [[cite:&HODGES2003]] written in material velocities and stresses, and a modal expansion of those variables is the final key step in seamlessly mapping the global FEM into the nonlinear description [[cite:&PALACIOS2011]]. The overall process can be found in [[cite:&CEA2021a]]. +- Maximum performance: as a combination of a highly optimised and vectorised codebase, numerical library JAX with its JIT compiler and accelerator capabilities driving the calculations, and the newly added added parallelisation of load cases. +- Differentiation and sensitivity analysis: using JAX algorithmic differentiation toolbox, the entire process, from inputs to aeroelastic outputs can be differentiated [[cite:&CEA2024a]]. -3. Force/moment modes, $\bm \Phi_2 = \mathcal{S}(\bm K_a \bm \Phi_0)$, represent the internal stress resultants in the structure as the sum, $\mathcal{S}$, along the main load-paths of equilibrium forces and moments produced by the modal deformations. - # Note that if $\bm{\mathfrak{f}} = \bm K_a \bm \Phi_0|^{1-3}$ are the internal forces and $\bm{\mathfrak{m}} = \bm K_a \bm \Phi_0|^{3-6}$ the internal moments, the moments produced by the internal forces also need to be taken into account: $\bm \Phi_2|^{3-6} = \mathcal{S}(\bm{\mathfrak{m}} + \bm{r}_{\frac{1}{2}} \times \bm{\mathfrak{f}})$. - Results are presented in the mid-point between nodes because more information cannot be extracted in terms of linear stresses from one node to the other. - -4. Strain modes, $\bm \Psi_2 = -\bm \Phi_{0d} + \pmb{E}^{\top}\bm \Phi_{0m}$, with $\bm \Phi_{0d}$ the approximate derivative along $s$: $\bm \Phi_{0d}^i = \frac{\bm \Phi_0^{i+1} - \bm \Phi_0^{i}}{\Delta s_i}$; and $\bm \Phi_{0m} = \frac{\bm \Phi_0^{i+1} + \bm \Phi_0^{i}}{2}$, the displacement modal shape in between nodes. $\pmb{E}^{\top}$ is a constant matrix. - After the intrinsic modes have been computed, a dynamic system is obtained after a Galerkin projection of the equations of motion \cite[Ch. 8]{PALACIOS2023}: - -\begin{equation} -\label{eq2:sol_qs} -\begin{split} -\dot{\pmb{q}}_{1} &= \pmb{\omega} \odot \pmb{q}_{2} - \pmb{\Gamma}_{1} \pmb{:} \left(\pmb{q}_{1} \otimes \pmb{q}_{1} \right) - \pmb{\Gamma}_{2} \pmb{:} \left( \pmb{q}_{2} \otimes \pmb{q}_{2} \right) + \bm{\eta} \\ -\dot{\pmb{q}}_{2} &= -\pmb{\omega} \odot \pmb{q}_{1} + \pmb{\Gamma}_{2}^{\top} \pmb{:} \left( \pmb{q}_{2} \otimes \pmb{q}_{1} \right) -\end{split} -\end{equation} -where $\odot$ is the Hadamard product (element-wise multiplication), $\otimes$ is the tensor product operation and $\pmb{:}$ is the double dot product[fn:1: The double dot product represents a contraction of the last two indexes of the first tensor with the first two indexes of the second one; it however needs further specification as two alternative definitions can be adopted and here we opt for the following: \(\pmb{a} \pmb{:} \pmb{b} = a_{..ij} b_{ij..} \). This has implications on the definition of the transpose of \(\bm{\Gamma}_2 \) in the second equation since for high order tensors multiple transpose operators can be defined. Consistency is achieved by ensuring the dot product operation satisfies the following: \( \pmb{x} \cdot \left(\bm{\Gamma} \pmb{:} \left( \pmb{y} \otimes \pmb{z} \right) \right) = \pmb{y} \cdot \left(\bm{\Gamma}^{\top} \pmb{:} \left(\pmb{z} \otimes \pmb{x} \right) \right) \), which leads to the transpose of the third order tensor, \( \bm{\Gamma} = \Gamma^{ijk} \), as \( \bm{\Gamma}^{\top} = \Gamma^{jki} \).]. -The form of the equations in compact tensorial notation is in fact the way they have been implemented and vectorised. This description is geometrically-exact, with nonlinearities encapsulated in the modal couplings of the third-order tensors $\pmb{\Gamma}_{1}$ and $\pmb{\Gamma}_{2}$ (the former introduces the gyroscopic terms in the dynamics and the latter introduces the strain-force nonlinear relation). $\pmb{\eta}$ is the modal projection of the external forcing terms. They are computed as integrals along the load-paths as an inner product: $\langle \pmb{u},\pmb{v} \rangle = \int_\Gamma \pmb{u}^\top \pmb{v} ds$, for any $\pmb{u}\in\mathbb{R}^6$ and $\pmb{v}\in\mathbb{R}^6$: -# ?? messy in tensorial notation (computation is via vmap) -\begin{align}\label{eq2:gammas12} -\Gamma_{1}^{ijk} & = \langle \pmb{\Phi}_{1i}, \lo_1(\pmb{\Phi}_{1j})\pmb{\Psi}_{1k}\rangle, \nonumber \\ -\Gamma_{2}^{ijk} & = \langle \pmb{\Phi}_{1i}, \lo_2(\pmb{\Phi}_{2j})\pmb{\Psi}_{2k}\rangle, \\ -\eta_{i} & = \langle \pmb{\Phi}_{1i}, \pmb{f}_1\rangle \nonumber -\end{align} -with $\lo_1$ and $\lo_2$ linear operators. The solution of Eqs. \ref{eq2:sol_qs} correspond to Step 3 in Fig. [[workflow]], and can be extended to form the full aeroelastic system with gravity forces, $\bm{\eta}_g$, aerodynamic forces and gust disturbances, $\bm{v}_g$. Control states can also be included [[cite:&CEA2021a]], but they are not necessary for this work. For a set of reduced frequencies and a given Mach number, the DLM (or a higher fidelity aerodynamic method) yields the Generalised Aerodynamic Forces (GAFs). The current implementation uses Roger's rational function approximation to those GAFs [[cite:&Roger1977]], which results in the follower modal forces: - -# \begin{equation} -# \lo_1 (\pmb{x}_1) = \begin{bmatrix} \tilde{\pmb{\omega}} & \pmb{0} \\ \tilde{\pmb{\ve}} & \tilde{\pmb{\omega}} \end{bmatrix} \hspace{0.5cm} ; \hspace{0.5cm} -# \lo_2 (\pmb{x}_2)= \begin{bmatrix} \pmb{0} & \tilde{\pmb{f}} \\ \tilde{\pmb{f}} & \tilde{\pmb{m}} \end{bmatrix} \hspace{0.5cm} ; \hspace{0.5cm} \pmb{\mathsf{E}}= \lo_1 \begin{pmatrix} \begin{bmatrix} 1 \\ \bm 0_5 \end{bmatrix} \end{pmatrix} -# \end{equation} +** Nonlinear aeroelastic system +Given a general GFEM, a reduced model is obtained from a static or dynamic condensation that captures well the stiffness and inertia properties in the condensed matrices, $\pmb{K}_a$ and $\pmb{M}_a$. The eigenvalue solution of the FEM yields the modal shapes, $\pmb \Phi_0$, and frequencies $\pmb \omega$. A projection of the state variables, velocities $\pmb{x}_1 = \pmb{\Phi}_1\pmb{q}_1$ and stresses $\pmb{x}_2 = \pmb{\Phi}_2\pmb{q}_2$, and a Galerkin projection of the equations of motion leads to the system of ODEs that is solved in time domain. +Aerodynamic forces are obtained via Generalised Aerodynamic Forces (GAFs) using a panel-based DLM solver and Roger's rational function approximation[[cite:&Roger1977]] to bring the forces to the time domain, resulting in a modal force component given as: \begin{equation}\label{eq3:eta_full} \begin{split} @@ -234,9 +196,8 @@ with $\lo_1$ and $\lo_2$ linear operators. The solution of Eqs. \ref{eq2:sol_qs} & \left. + \pmb{\mathcal{A}}_{g0}\bm{v}_g +\frac{c}{2U_\infty}\pmb{\mathcal{A}}_{g1} \dot{\bm{v}}_g +\left(\frac{c}{2U_\infty}\right)^2 \pmb{\mathcal{A}}_{g2}\ddot{\bm{v}}_g + \sum_{p=1}^{N_p} \pmb{\lambda}_p \right) \end{split} \end{equation} -Where the $\pmb{\mathcal{A}}_is$ are real matrices, $c$ is the reference chord, $\tfrac12\rho_\infty U_\infty^2$, $\pmb{\lambda}_p$ the aerodynamic states and $N_p$ the number of lags. -The coupling of the structure and aerodynamic equations combined with the aerodynamic lags gives the final ODE system: - +Where the $\pmb{\mathcal{A}}_is$ are real matrices, $c$ is the reference chord, $\tfrac12\rho_\infty U_\infty^2$, $\pmb{\lambda}_p$ the aerodynamic states and $N_p$ the number of lags. Note these forces naturally follow the structure given the formulation is written in the material frame of reference. +The coupling of the structure and aerodynamic equations combined with the aerodynamic lags, gravity forces, $\bm{\eta}_g$, and gust disturbances, $\bm{v}_g$, gives the final ODE system: \begin{equation} \label{eq2:sol_qs} \begin{split} @@ -247,166 +208,278 @@ The coupling of the structure and aerodynamic equations combined with the aerody -\frac{2U_\infty\gamma_p}{c}\bm{\lambda}_{p} \end{split} \end{equation} -in this system the aerodynamic added-mass effect has been moved to the left hand side such that $\bm{\mathrm{A}}_2 = (\pmb{I} - \frac{\rho c^2}{8}\pmb{\mathcal{A}}_2)^{-1}$, and it couples all DoF in $\pmb q_1$. Thus the natural frequency terms become $\hat{\pmb{\Omega}} = \bm{\mathrm{A}}_2 \textup{diag}(\pmb{\omega})$ and the nonlinear terms $\hat{\pmb{\Gamma}} = \bm{\mathrm{A}}_2 \bm{\Gamma}$. The effect of all external forces, aero, $\bm{\eta}_a$, gravity, $\bm{\eta}_g$, and others, $\bm{\eta}_f$, are combined in such that $\hat{\bm{\eta}} = \bm{\mathrm{A}}_2 \left( \left( \bm{\eta}_a - \frac{\rho c^2}{8} \pmb{\mathcal{A}}_2\dot{\bm{q}}_1 \right) + \bm{\eta}_g + \bm{\eta}_f \right)$. +where $\odot$ is the Hadamard product (element-wise multiplication), $\otimes$ is the tensor product operation and $\pmb{:}$ is the double dot product. +In this system the aerodynamic added-mass effect has been moved to the left hand side such that $\bm{\mathrm{A}}_2 = (\pmb{I} - \frac{\rho c^2}{8}\pmb{\mathcal{A}}_2)^{-1}$, and it couples all DoF in $\pmb q_1$. Thus the natural frequency terms become $\hat{\pmb{\Omega}} = \bm{\mathrm{A}}_2 \textup{diag}(\pmb{\omega})$ and the nonlinear terms $\hat{\pmb{\Gamma}} = \bm{\mathrm{A}}_2 \bm{\Gamma}$. The effect of all external forces, aero, $\bm{\eta}_a$, gravity, $\bm{\eta}_g$, and others, $\bm{\eta}_f$, are combined in such that $\hat{\bm{\eta}} = \bm{\mathrm{A}}_2 \left( \left( \bm{\eta}_a - \frac{\rho c^2}{8} \pmb{\mathcal{A}}_2\dot{\bm{q}}_1 \right) + \bm{\eta}_g + \bm{\eta}_f \right)$. The aerodynamic matrices $\hat{\bm{\mathcal{A}}}_{p+2}$ have also been scaled accordingly. + The nonlinearities in the system are encapsulated in the modal couplings of the third-order tensors $\pmb{\Gamma}_1$ and $\pmb{\Gamma}_2$ (the former introduces the gyroscopic terms in the dynamics and the latter introduces the strain-force nonlinear relation). \\ - -In the second instance, the rotation and position in the inertial reference system are calculated by integration of strains along the domain, as in the Frenet-Serret formulas of differential geometry. Following definition of strains and curvatures, -\begin{equation}\label{eq2:urecover_s} -\begin{split} -\pmb{R}_{ab}^{\prime} &= \pmb{R}_{ab}\tilde{\pmb{k}} \\ -\pmb{r}_a'&=\pmb{R}_{ab}(\pmb{\gamma} + \pmb{e}_x) -\end{split} -\end{equation} -Analytical solutions to Eq. \eqref{eq2:urecover_s} can be obtained when the strain is assumed constant between nodes and a piecewise constant integration is carried out, as is the case in the current implementation. If the beam path is discretized in n+1 points, strain and curvatures are defined in the mid-points of the spatial discretization (n in total). $\gamma_n$ and $\kappa_n$ are constant within the segment $s_{n-1} \leq s \leq s_n$, and the position and rotation matrix after integration are -\begin{equation}\label{eq:strain_integration} -\begin{split} -\bm{R}_{ab}(s) &= \bm{R}_{ab}(s_{n-1})\pmb{\mathcal{H}}^0(\bm{k},s) \\ -\bm{r}_a(s) &= \bm{r}_a(s_{n-1}) + \bm{R}_{ab}(s_{n-1})\pmb{\mathcal{H}}^1(\bm{k}, s)\left(\bm{e}_x+\bm{\gamma}_n\right) -\end{split} -\end{equation} -Where the operators $\pmb{\mathcal{H}}^0(\bm{k}, s)$ and $\pmb{\mathcal{H}}^1(\bm{k}, s)$ are obtained from integration of the exponential function as defined in \cite{Palacios2010}. -# \begin{equation} -# \begin{split} -# \pmb{\mathcal{H}}^0(\bm{k},s) &= e^{\Delta\tilde{\bm{k}}} =\pmb I + \frac{\sin(\Delta \phi)}{\Delta{ \phi}}\Delta \tilde{\pmb\Psi}+ \frac{1-\cos(\Delta \phi)}{\Delta \phi^2}\Delta \tilde{\pmb\Psi}\Delta \tilde{\pmb\Psi} \\ -# \pmb{\mathcal{H}}^1(\bm{k},s) &= \Delta s\left(\pmb I + \frac{1-\cos(\Delta \phi)}{\Delta \phi^2}\Delta \tilde{\pmb\Psi} + \frac{\Delta \phi -\sin(\Delta \phi)}{\Delta \phi^3}\Delta\tilde{\pmb\Psi}\Delta\tilde{\pmb\Psi} \right) -# \end{split} -# \end{equation} -# with $\Delta s = s- s_{n-1}$, $\Delta \pmb{\Psi} = \bm{k} \Delta s$ and $\Delta \phi = ||\Delta \pmb{\Psi}||$. -Note that when position and rotations are recovered from strain integration, there is still one point that is either clamped or needs to be tracked from integration of its local velocity -In the next section an optimized implementation of this algorithm is shown in JAX. -\\ -Once the nonlinear solution of the condensed model is computed, the corresponding full 3D state is calculated via a two postprocessing steps: firstly the displacements of the cross-sectional nodes linked to the reduced model via the interpolation elements are computed using the positions and rotations of the latter; secondly, Radial Basis Functions (RBFs) kernels are placed on those cross-sections, thus building an intermediate model that is utilised to extrapolate the positions of the remaining nodes in the full model. +Once the nonlinear solution of the condensed model is computed, the corresponding full 3D state is calculated via two postprocessing steps: firstly the displacements of the cross-sectional nodes linked to the reduced model via the interpolation elements are computed using the positions and rotations of the latter; secondly, Radial Basis Functions (RBFs) kernels are placed on those cross-sections, thus building an intermediate model that is utilised to extrapolate the positions of the remaining nodes in the full model. This paves the way for a broader multidisciplinary analysis where CFD-based aerodynamic loading could be used for the calculation of the nonlinear static equilibrium, and also with the transfer of the full deformed state back to the original FE solver to study other phenomena such as local buckling. -* Results -The following results show the strength of the approach in problems with large geometric nonlinearities, the ability to run on CPUs and GPUs and to automatically compute gradients of dynamic problems, and the edge in performance when compare with commercial toolbox Nastran. -** Automatic differentiation of the nonlinear dynamics of a representative wing -This test case demonstrates the accuracy of the NMROM approach for dynamic geometrically-nonlinear calculations and was first introduced in [[cite:&CEA2021b]]. The right wing of Fig. [[fig:SailPlane2]] is considered and dynamic nonlinear simulations are carried out and compared to MSC Nastran linear and nonlinear analysis (SOL 109 and 400, respectively) on the full FE model. +** Computational implementation +# Bring parallelisation on distributed accelerators into our solution process, thus applying the state-of-the-art techniques in used for the large problems in deep-learning. Combined with our already fast simulations times, this could allow the prediction of those sizing aeroelastic loads that include thousands of cases in commercial aircraft, the computation of their gradients with respect to design variables, with geometrically nonlinear effects accounted for, and at such performance that the framework could be integrated into a larger multidisciplinary optimization. + +# One of the main contribution of this work is a new computational implementation that achieves accelerations of over 2 orders of magnitude with respect to its predecessor \footnote{Both the new implementation and the examples of this paper can be found at \url{https://github.com/ACea15/FENIAX}}. In addition, a highly modular, flexible architecture based on software design patterns has been put in place, which was further described in \cite{CEA2024}. Moreover, the resulting nonlinear aeroelastic framework is suitable for modern hardware architectures and able to compute sensitivities via algorithmic differentiation (AD), as will be demonstrated herein. +# The key enabler was moving from standard Python to a highly vectorised, JAX-based numerical implementation. JAX is a Python library designed for high-performance numerical computing with focus on machine learning activities \cite{jax2018github}. It combines XLA (accelerated linear algebra) and Autograd, the former being a compiler that optimises models for different hardware platforms, the latter is an Automatic Differentiation (AD) tool in Python. +# The XLA compiler orchestrates the conversion of high-level Python code into efficient, low-level machine-specific instructions. When it comes to leveraging the computational power of GPUs, the link between XLA and CUDA kernels is critical. Here's how these components interact: -#+NAME: fig:SailPlane2 -#+CAPTION: Representative plane structural and aerodynamic models -#+ATTR_LATEX: :width 0.7\textwidth -[[file:figs_ext/SailPlaneRef.png]] +# 1. **JAX and XLA**: JAX is a numerical computing library that provides high-level interfaces for differentiable programming. It incorporates automatic differentiation and just-in-time (JIT) compilation to optimize performance. The JIT compilation is powered by XLA. -A force is applied at the wing tip with a triangular loading profile, followed by a sudden release of the applied force to heavily excite the wing. The force profile is given in Fig. [[fig:ramping_load]]. The applied force is then \(f_{tip} = \alpha \times \textup{\pmb{f}}_{max} f(0.05, 4) = [-2\times 10^5, 0., 6\times 10^5]f(0.05, 4)\) where $\alpha$ has been set to $1$. +# 2. **XLA (Accelerated Linear Algebra)**: XLA is a domain-specific compiler that takes the computation expressed in JAX, optimizes it, and targets it for execution on specific hardware, such as CPUs, GPUs, or TPUs. It performs optimizations such as operation fusion, constant folding, and reducing memory transfers, which are crucial for high-performance computing. -#+NAME: fig:ramping_load -#+CAPTION: Ramping load profile for dynamic simulation of representative wing -#+ATTR_LATEX: :width 0.6\textwidth -[[./figs_ext/ramping_load.pdf]] -The dynamic response is presented in Fig. [[fig:wsp_3d]], where results have been normalised with the wing semi-span (28.8 m.). As expected, linear analysis over-predicts vertical displacements and does not capture displacements in the $x$ and $y$ directions. NMROMs were built with 5, 15, 30, 50 and 100 modes. A Runge-Kutta four is used to march the equation in time with time steps corresponding to the inverse of the largest eigenvalue in the NMROM, i.e. $dt = [27.34, 6.62, 2.49, 1.27, 0.575] \times 10^{-3}$ s. +# 3. **Targeting GPUs and CUDA Kernels**: +# - **Translation to HLO (High-Level Optimizations)**: XLA first converts the JAX computational graphs into an intermediate representation called HLO. This provides an opportunity to perform high-level optimizations on the computations before they are translated into device-specific code. +# - **Mapping to CUDA Kernels**: For GPU execution, XLA translates HLO into GPU-specific instructions. This involves generating CUDA kernels that can be executed by NVIDIA GPUs. CUDA kernels are the fundamental units of execution in the CUDA programming model provided by NVIDIA, allowing parallel execution of computations on the GPU. +# - **Execution on GPUs**: Once XLA generates the optimized CUDA kernels, they are executed on the GPU. This process leverages the parallel processing capabilities of the GPU, allowingJAX to efficiently handle large-scale computations. -#+NAME: fig:wsp_3d -#+ATTR_LATEX: :width 1\textwidth -#+CAPTION: Span-normalised wing-tip displacements -#+RESULTS: WSPsubplots -[[file:figs/WSPsubplots.png]] +# In summary, the link between XLA and CUDA in the context of JAX involves the transformation of high-level computation graphs into optimized CUDA kernels via XLA's compilation pipeline. This is what enables JAX to achieve high performance when running on NVIDIA GPUs, harnessing their parallel computing architecture. + +# The required abstractions are by constraining the python code a functional programming paradigm + +# Concurrent simulations in multi-GPU environments involve distributing and executing computations across multiple GPU devices to leverage their collective computational power. In JAX, achieving this involves several strategies and technologies. Here's how concurrent simulations are facilitated: + +# 1. **JAX's Support for Parallelism**: +# - **pmap (Parallel Map)**: JAX offers a `pmap` function designed for parallel execution across multiple devices, including GPUs. `pmap` allows you to automatically batch computations across different devices in parallel. It effectively maps a function across multiple input sets, distributing the workload across available GPUs. +# - **Replication**: Typically with `pmap`, JAX replicates the computational workload across multiple GPUs, each performing the computation on a subset of the data. + +# 2. **XLA's Role in Multi-GPU Execution**: +# - XLA splits the high-level operations into device-specific computations, managing inter-device communication and synchronization. This includes optimization and coordination of data transfer between devices to minimize overhead. +# - XLA's compiler applies device-specific optimizations to ensure efficient workload distribution, synchronization, and communication between GPUs. + +# 3. **Data Sharding**: +# - Data can be explicitly sharded across GPUs. Each GPU processes a portion of the data, which is crucial for reducing memory overhead and increasing throughput. +# - Efficient data sharding and collection strategies are critical, especially to ensure that no GPU becomes a bottleneck due to data transfer overheads. + +# 4. **Custom Collective Operations**: +# - JAX provides APIs for defining custom collective operations. These operations include complex data communications, such as all-reduce operations necessary for aggregating mirrored data across GPUs. +# - This allows for synchronization and aggregation of results from each device, essential for consistent updates in simulations or model training. + +# 5. **Inter-Device Communication**: +# - JAX manages low-level inter-device communication, often using NVIDIA's NCCL (NVIDIA Collective Communications Library) for efficient peer-to-peer data transfers and collective operations like broadcasts and reductions. + +# 6. **Resource Management**: +# - Proper allocation and management of GPU resources, such as memory allocation and computational tasks, are crucial to maximizing throughput and minimizing potential conflicts or resource contention. -The 3D shape of the model is retrieved and compared against the full nonlinear dynamic solution as illustrated in Fig. [[wsp_3d]] (Nastran solution in yellow and NMROM with 50 modes in blue). The times at positive and negative peaks are displayed. Even though a wing of such characteristics would never undergo this level of deformations, these results further support the viability of the methodology to solve highly geometrically nonlinear dynamics, on complex models and with minimal computational effort. -#+NAME: wsp_3d -#+CAPTION: Snapshots of wing 3D dynamic response comparing NMROM (yellow) and NLFEM3D (blue) -#+ATTR_LATEX: :width 1\textwidth -[[./figs_ext/WSP_3D-front.png]] -Next we look at the differences of the dynamic simulations with the same metric employed above that now evolves in time. Integrator errors accumulate and discrepancies grow with time but still remain small. In fact the differences between Nastran and our dynamic solvers are comparable to the static example with the highest load (around the $5\times 10^{-5}$ mark), both cases inducing over 25\% percent deformations of the wing semi-span. -#+NAME: WSP_error -#+CAPTION: L-2 norm per node differences between Nastran full FE solution and NMROM with 50 modes +# In JAX, you need to only specify how you want the input and output of your code to be partitioned, and the compiler will figure out how to: 1) partition everything inside; and 2) compile inter-device communications. + +# The XLA compiler behind jit includes heuristics for optimizing computations across multiple devices. In the simplest of cases, those heuristics boil down to computation follows data. + +# To demonstrate how auto-parallelization works in JAX, below is an example that uses a jax.jit()-decorated staged-out function: it’s a simple element-wise function, where the computation for each shard will be performed on the device associated with that shard, and the output is sharded in the same way: + +# NVIDIA Collective Communications Library (NCCL) + + +# In JAX, the process of utilizing the GPU involves both generating custom CUDA kernels and leveraging existing libraries for standard operations. Here's how that typically works: + +# 1. **Custom CUDA Kernel Generation**: +# - JAX, through XLA, can generate custom CUDA kernels for specific operations that are not efficiently covered by existing libraries. This is particularly true for operations that benefit from domain-specific optimizations. +# - The XLA compiler analyzes the computational graph and determines whether generating a custom kernel is beneficial, considering factors like operation fusion, data locality, and computation complexity. + +# 2. **Using Existing Optimized Libraries**: +# - For many standard and well-understood operations (e.g., matrix multiplications, convolutions, basic arithmetic operations), JAX relies heavily on highly optimized libraries rather than generating new kernels from scratch. +# - **cuBLAS**: For linear algebra operations (e.g., matrix multiplication), JAX calls functions from the NVIDIA cuBLAS library, which is optimized for performance on NVIDIA GPUs. +# - **cuDNN**: For deep learning operations like convolutions, JAX uses the cuDNN library, which provides highly optimized implementations of commonly used deep learning primitives. +# - **CUB**: JAX may use NVIDIA CUB (CUDA Unbound) for efficient operations on parallel processors, such as reduces, scans, and histogram operations. + +# 3. **Fusion and Optimization**: +# - JAX, via XLA, aims to minimize the overhead of launching kernels by fusing operations. Instead of launching multiple kernels for a sequence of operations, XLA attempts to combine them into a single operation when possible, reducing memory transfers and improving execution efficiency. + +# 4. **Hybrid Approach**: +# - Often, JAX will employ a hybrid approach, using acombination of custom kernels for unique, complex operations while leveraging optimized library calls for standard computations. This balance helps capitalize on existing optimization work while providing flexibility for custom performance tuning when needed. + +# In summary, JAX effectively combines the use of CUDA kernel generation for specific needs and the efficient calls to established libraries for common operations, aiming to provide both flexibility and high performance when running computations on GPUs. + +# device parallelism for Single-Program Multi-Data (SPMD) code in JAX. SPMD is a parallelism technique where the same computation, such as the forward pass of a neural network, can be run on different input data (for example, different inputs in a batch) in parallel on different devices, such as several GPUs + +Algorithm [[alg:process]] shows the main components in the solution process, highlighting the time and space complexities, $O(time, space)$, of the data structures being generated. We assume a single analysis is being run, for instance a dynamic simulation computing the response to multiple gusts that will be run in parallel for a total number of $N_c$ cases. $N_t$ time-steps are used in the integration scheme with a resolution of $N_m$ modal shapes. The FE model has been condensed to $N_N$ number of nodes. + + +#+NAME: alg:process +\begin{algorithm}[h!] +\DontPrintSemicolon +\SetKwInOut{Input}{input} +\SetKwInOut{Output}{output} +\Input{Input file: settings.yaml; FE model: $\bm{K}_a$, $\bm{M}_a$, $\bm{X}_a$; Aerodynamic matrices: $\bm{\mathcal{A}}$} +\Output{Nonlinear aeroealastic solutioxn} +\Begin{ + \BlankLine +$\bm{\phi}$, $\bm{\psi}$ $\longleftarrow$ modes($\bm{K}_a$, $\bm{M}_a$, $\bm{X}_a$) \Comment{Intrinsic modes: O($N_n^2 \times N_m$; $N_n \times N_m$)} \; +$\bm{\Gamma}$ $\longleftarrow$ couplings($\bm{\phi}$, $\bm{\psi}$) \Comment{Nonlinear couplings O($N_n \times N_m^3$; $N_m^3$)} \; +$\bm{q}$ $\longleftarrow$ system($\bm{\Gamma}$, $\bm{\mathcal{A}}$, $\bm{\phi}$, $\bm{X}_a$) \Comment{Modal coordinates: O($\frac{N_c}{N_d} \times N_t \times N_m^3$; $N_c \times N_t \times N_m$)} \; +$\bm{X}_1$, $\bm{X}_{2}$, $\bm{X}_{3}$ $\longleftarrow$ ivars($\bm{q}$, $\bm{\phi}$, $\bm{\psi}$) \Comment{velocity/strain fields: O($\frac{N_c}{N_d} \times N_t \times N_n \times N_m$; $N_c \times N_t \times N_n$)} \; +$\bm{r}_a$, $\bm{R}_{a}$ $\longleftarrow$ integration($\bm{X}_{3}$, $\bm{X}_a$) \Comment{Positional/rotational fields: O($\frac{N_c}{N_d} \times N_t \times N_n \times N_m$; $N_c \times N_t \times N_n$)} \; +\BlankLine +} +\caption{Main components in solution process} +\end{algorithm} + +*** Two-level parallelisation + +* Results +In this section we show the main strengths of our solvers to: a) run a representative aircraft model undergoing very large nonlinear displacements; b) leverage on modern hardware architectures and a parallelisation across devices to unlock problems such as quantifying the uncertainties in the nonlinear response given a loading field that is not fully determinate; c) build load envelopes of static and dynamic aeroelastic simulations. +The University of Bristol Ultra-Green (BUG) aircraft model [[cite:&STODIECK2018]] is the chosen platform to demonstrate these capabilities as it showcases high-aspect ratio wings that are built using a representative GFEM of current industrial models and it is available open-source. The main components of the aeroelastic model are shown in Fig. [[fig:BUG]]. The GFEM is formed of over 16000 grid nodes, 17000 Quad and Triangular shell elements, 16000 beam elements, 400 discrete mass elements, 30 rigid elements and 100 interpolation elements. + +#+NAME: fig:BUG +#+CAPTION: BUG model GFEM and DLM models +#+ATTR_LATEX: :width 1\textwidth +[[file:figs_ext/bug_model7.pdf]] + + +# A modal analysis is first showing + +Structural and aeroelastic static simulations follow, all solved via a Newton-Raphson solver with tolerance of $10^{-6}$, as well as an assessment of the aircraft dynamics in response to a gust. +Calculations are carried out on a CPU Intel Xeon Silver 4108 with 1.80GHz speed, 6 cores and a total 12 threads, as well as on an Nvidia GPU A100 80GB SXM. +** Structural static analysis +Two exercises are studied to assess two levels of parallisation in the current implementation. That corresponding to the operations within a single solution, and that of multiple load cases. +Firstly, a tip load is prescribed as to induce very large deformations such that a big modal basis is needed to accurately capture the response. These extreme cases are very good to be solved in modern hardware architectures that can run many of the operations involving tensors in parallel. Secondly, uncertainty quantification is performed to the nonlinear response to a loading field that is non-deterministic. Hundreds to thousands of simulations are employed in Monte Carlo type of analysis to resolve for the statistics, for which parallelisation of the independent simulations become critical. +*** Extremely large deformations under discrete loads +A total of eight different loading cases are computed with tip loads as forces and moments in the $x, y, z$ directions, and a combination of both. +The number of modes used in the solution is 300. +Even though a wing would never undergo such a level of deformations before breaking, we are showing them to demonstrate to both the efficiency and robustness of the solvers. +Fig. [[fig:BUG_tipL2]] shows the static equilibrium for follower tip forces in the $z-$ direction of 2.5 \times 10^4, 4.5 \times 10^4, 1.5 \times 10^5, 2.125 \times 10^5 Newtons; similarly in Fig. [[fig:BUG_tipL0]] $x-$ direction forces of 1.25 \times 10^5, 2.25 \times 10^5, 7.5 \times 10^5, 1.125 \times 10^6 are applied; and Fig. [[fig:BUG_tipL4]] results from the application of torsional moments from 1.75 \times 10^5 to 1.4 \times 10^6 N m. +# 1, 3, 7, 10 -> 2.5e4, 4.5e4, 1.5e5, 2.125e5 +# 1, 3, 7, 11 -> 1.25e5, 2.25e5, 7.5e5, 1.125 e6 +# 1,3,5,6,7,8,9 -> 1.75e5 - 1.4e6 N m +#+NAME: fig:BUG_tipL2 +#+CAPTION: Static equilibrium for out-of-plane tip loads +#+ATTR_LATEX: :width 0.8\textwidth +[[file:figs_ext/DiscreteL2.png]] + +#+NAME: fig:BUG_tipL0 +#+CAPTION: Static equilibrium for in-plane tip loads #+ATTR_LATEX: :width 0.7\textwidth -#+RESULTS: WSP_error -[[file:figs/WSP_error.png]] +[[file:figs_ext/DiscreteL0.png]] -An impressive reduction of computational time is achieved by our solvers as highlighted in Table [[table:WSP_times]]. The nonlinear response of the full model in Nastran took 1 hour 22 minutes, which is over two orders of magnitude slower than the NMROM with 50 modes resolution, which proved very accurate. The significant increase in computational effort when moving from a solution with 50 modes to 100 modes is due to various factors: vectorised operations are limited and the quadratic nonlinearities ultimately lead to O($N_m^3$) algorithms; the time-step needs to be decreased for the Runge-Kutta integration to remain stable; the additional overheads that come with saving and moving larger tensors, from the modal shapes, the cubic modal couplings, to the system states (note times shown account for all the steps from start to end of the simulation, including saving all the data for postprocessing). +#+NAME: fig:BUG_tipL4 +#+CAPTION: Static equilibrium for torsional tip loads +#+ATTR_LATEX: :width 0.8\textwidth +[[file:figs_ext/DiscreteL4.pdf]] -#+NAME: table:WSP_times -#+CAPTION: Computational times representative wing dynamic solution -#+ATTR_LATEX: :center t -#+RESULTS: WSP_times -| | NMROM (modes: 5, 15, 30, 50, 100) | NASTRAN 400 | NASTRAN 109 | -|----------+-----------------------------------+-------------+-------------| -| Time [s] | 2.79, 2.92, 4.85, 7.14, 155.3 | 4920 | 33.6 | +Table [[table:times_discrete]] presents the computational times for a single load case in the CPU and therefore the time it would take to run the 8 cases sequentially; for the CPU to run the 8 loads within the parallelised algorithm; and for the GPU to run in parallel. +It worth noting the parallelisation on the CPU does not attain a very good speed up, the reason being the modal base is so large and some of the tensor calculations in JAX, which already run in concurrently by default, will be taking 6 cores and 12 threads available in the CPU. The A100, one of the latest GPUs of for high performance computations, has 6912 CUDA cores that deliver a remarkable increase in performance. -We move now to one of the main highlights of this work, i.e. the ability to compute gradients via automatic differentiation in geometrically nonlinear dynamic problems. The maximum root loads occurring in a wing subjected to dynamic loads is a good test case as it can be a critical metric in sizing the aircraft wings, especially high-aspect ratio ones. Thus we look at the variation of the maximum z-component of the vertical internal forces as a function of \(\alpha\) in the loading profile of Fig. [[fig:ramping_load]]. Effectively, the slope of the loading increases with \(\alpha\). Table [[table:AD_WSP]] shows the derivatives computed using FD with $\epsilon=10^{-4}$ and AD in reverse-mode on the example with 50 modes resolution. In this case the FD needed tweaking of $\epsilon$ while application of AD was straight forward with no need for checkpoints and took around three times the speed of a single calculation. +#+CAPTION: Computational times tip-forces +#+ATTR_LATEX: :center t +#+NAME: table:times_discrete +| Device | Time (sec.) | +|----------------+-------------------| +| CPU (single) | 31 \times 8 = 248 | +| CPU (parallel) | 188.1 | +| GPU | 24.2 | + +*** Uncertainty quantification of nonlinear response +This example is manufactured to resemble the workflow of flight loads and wing stress analysis in an industrial setup: the flight physics department would compute the in-flight loads for various conditions and pass the maximum of these loads to the stress engineers who would check the integrity of the airframe in their more detailed models. There will always be an element of uncertainty around the loads computed by the flight department, and what we show here is how for nonlinear assumptions the statistics need to be computed for every distinct loading. And for this, having a parallisation strategy as the one presented could potentially allow the computation of complex correlations and averages that are more easily calculated under linear assumptions. +A constant loading force is prescribed along the wings consisting of follower forces in the $z-$ direction as well as torsional moments, with the characteristic that the force follows a normal distribution with $N(\mu=1.5 \times 10^4 \times \mu_0, \sigma=0.15 \times \mu)$ for the vertical forces and $N(\mu=3 \times 10^4 \times \mu_0, \sigma=0.15 \times \mu)$ for the moments. Three scenarios are studied: one in which very large nonlinear deformations are induced with $\mu_0 = 1$, and two small loading with $\mu_0 = 10^{-2}$ and $\mu_0 = 10^{-3}$. +The distribution of displacements is characterised by means of Montecarlo simulations that run in parallel for a total of 1600 simulations. The modal resolutions consists of 100 modes. +Fig. [[fig:BUG_mc]] shows the equilibrium for the high loading calculations for two cases out of the 1600. +#+NAME: fig:BUG_mc +#+CAPTION: Static equilibrium for two cases of the random excitation +#+ATTR_LATEX: :width 0.8\textwidth +[[file:figs_ext/MC1.png]] -#+NAME: table:AD_WSP -#+CAPTION: Automatic differentiation in dynamic problem +Table [[table:BUG_mc]] shows the statistics gathered from the response +#+CAPTION: Tip displacement statistics #+ATTR_LATEX: :center t -| $\alpha$ | $f(\alpha)$ [KN/m] | $f'(\alpha)$ (AD) | $f'(\alpha)$ (FD) | -|----------+----------------------+-----------------------+-----------------------| -| 0.5 | $1723.2 \times 10^3$ | $3587.71 \times 10^3$ | $3587.77 \times 10^3$ | -| 1.0 | $3624.4 \times 10^3$ | $3735.26 \times 10^3$ | $3735.11 \times 10^3$ | -| 1.5 | $5608.3 \times 10^3$ | $3957.81 \times 10^3$ | $3958.31 \times 10^3$ | - -** Accelerator benchmark on a very flexible free flying structure -This example exemplifies the ability of our solvers to turn a generic linear free-free finite-element model into a fully nonlinear solution that accounts for the rigid-body dynamics coupled with large elastic deformations, which has already been presented in [[cite:&PALACIOS2019]]. The novelties introduced herein are the new optimised implementation that can run on accelerators and the approach to recover the full 3D state from the reduced model. -The beam version of this structure was first studied by Simo and Vu-Quoc [[cite:&SIMO1988]] and has served to verify several implementations of nonlinear beam dynamics with rigid body motions [[cite:&HESSE2014]]. -A straight structure of constant square cross section (side = 3, wall thickness = 3/10) is built consisting of 784 shell elements linked to 50 spanwise nodes via interpolation elements as depicted in Fig. [[fig:FFS]] together with the material properties and two types of loading: firstly, a dead-force in the x-direction and dead-moment in the z-direction that yield a planar motion in the x-y plane; and secondly, the addition of a moment in the y-direction which produces a three dimensional motion. - -#+NAME: fig:FFS -#+CAPTION: FFS geometry, material properties and load cases -#+ATTR_LATEX: :width 0.7\textwidth -[[file:figs_ext/ffbw10.pdf]] +#+NAME: table:BUG_mc +| Case | Tip displacement mean (m) | Tip displacement std | +|-------------------------------+---------------------------+----------------------| +| Nonlinear ($\mu_0 = 1$) | 11.57 | 1.35 | +| Linear ($\mu_0 = 0.01$) | 0.148 | 0.024 | +| Very Linear ($\mu_0 = 0.001$) | 0.0149 | 0.0023 | + + +# Mean displacement node 35: 11.566769265603666 +# std displacement node 35: 1.3448662385231276 +# Ratio displacement node 35: 8.600683796111781 +# *************** +# Mean displacement node 35: 0.14768956221710616 +# std displacement node 35: 0.024150658437415644 +# ratio displacement node 35: 6.115343090948471 +# *************** +# Mean displacement node 35: 0.01485757200729988 +# std displacement node 35: 0.002342569483498701 +# ratio displacement node 35: 6.342425320554263 +# *************** + +We can see the statistics of the linear response are fully captured by one example, whereas for a nonlinear response such as $\mu_0 = 2$, the 1600 simulations would need to be computed again. Table [[table:times_MC]] shows the times taken for the nonlinear case. The computation of 1600 independent simulations of Fig. [[fig:BUG_mc]], which presents deformations of over 40% the wing semi-span, in just over a minute, highlights the potential of this methodology in more complex uncertainty quantification problems. + +#+CAPTION: Computational times uncertainty quantification +#+ATTR_LATEX: :center t +#+NAME: table:times_MC +| Device | Time (sec.) | +|----------------+--------------------------| +| CPU (single) | 16.8 \times 1600 = 26880 | +| CPU (parallel) | 317.4 | +| GPU | 67.6 | + -The free-flying evolution of the 3D model is shown in Fig. [[fig:FFB_2D]] for the planar motion and Fig. [[fig:FFB_3D]] for the loads giving rise to the full 3D deformations. It worth remarking the later motion also exhibits large torsional deformations which are combined with the also large bending displacements and rigid-body modes. +** Steady manoeuvre loads +We extend the analysis to an static aeroelastic case for varying angles of attack that represent a manoeuvre scenario. The number of modes used was 100, more than necessary for this type of response, which indicates even faster calculations are possible on this type of analysis. We test the parallelisation by varying the flow density ($\pm 20 \%$ of the reference density 0.41 Kg/ m$^3$) as well and the flow velocity ($\pm 20 \%$ of the reference velocity 209.6 m/s). 16 different points for both density and velocity make a total number of 256 simulations. -#+NAME: fig:FFB_2D -#+CAPTION: Free-flying structure in the 2D plane +Fig. [[fig:BUG_manoeuvre3D]] illustrates the 3D equilibrium of the airframe at the reference flow values. +#+NAME: fig:BUG_manoeuvre3D +#+CAPTION: Aeroelastic steady equilibrium for increasing angle of attack #+ATTR_LATEX: :width 0.8\textwidth -[[file:figs_ext/FFB_2D2.pdf]] +[[file:figs_ext/monoeuvre3D.pdf]] +# Fig. [[fig:BUG_manoeuvreRa]] shows the tip displacement with angle of attack in both the parallel and a single simulation run normally for validation. +# #+NAME: fig:BUG_manoeuvreRa +# #+CAPTION: Aeroelastic steady equilibrium for increasing angle of attack +# #+ATTR_LATEX: :width 0.8\textwidth +# [[file:figs_ext/monoeuvre3D.pdf]] -#+NAME: fig:FFB_3D -#+CAPTION: Free-flying structure in the 3D plane -#+ATTR_LATEX: :width 1.\textwidth -[[file:figs_ext/FFB_3Dt.pdf]] +The simulations show displacements differ if compare to a linear projection for the higher angles, which highlight the potential need for geometrically nonlinear aeroelastic tools in future aircraft configurations under high loading scenarios. +# As the angle of attack, AoA, is increased, the tip displacement falls down the linear projection between the 0 and 0.5 AoA as expected. +Table [[table:times_manoeuvre]] shows the computational times of this case, which shows near no overhead in adding a few hundred of static calculations when moving from the single load case in the CPU to the GPU (nearly 8 seconds to 14 seconds, which amounts for 6 seconds cost when adding an extra 255 cases). + +#+CAPTION: Computational times multiple manoeuvre problem +#+ATTR_LATEX: :center t +#+NAME: table:times_manoeuvre +| Device | Time (sec.) | +|----------------+--------------------------| +| CPU (single) | 7.71 \times 256 = 1973.8 | +| CPU (parallel) | 52.8 | +| GPU | 14.4 | -Because the applied load is a dead force we can track the position of the center-of-gravity (Cg) analytically as a verification exercise. Furthermore, the highly nonlinear nature of this problem makes it a good example to showcase the strength of accelerators for large problems and to gain insights as to when it might be better to deploy the codes in standard CPUs instead. Therefore we perform a sweep with the number of modes kept in the solution from 50 to 300, which determines the size of the system to be solved. The full modal basis is employed at 300 modes and due to the nonlinear cubic term this entails operations of the order of $O(10^7)$ at every time step of the solution, making it a good case for accelerators. The increase in the number of modes also restricts the incremental time step used in the explicit solver to preserve stability. Table [[table:FFB_times]] shows both computational time and Cg error for the planar case and in two scenarios: linking the integration time-step to the largest eigenvalue $\lambda$ in the solution $dt=\lambda^{-0.5}$; and fixing it to $dt=10^{-3}$. -Computations have been carried out in AMD EPYC 7742 CPU processors and Nvidia GPU RTX 6000 at the Imperial College cluster. -# time steps = 0.001, 0.0028, 0.0061 +** Dynamic loads at large scale +In this final example we perform a dynamic aeroelastic analysis to study the response of the aircraft to multiple 1-cos gusts for varying length, intensity and the density of the airflow. The mach number is kept constant at 0.7. A Rungue-Kutta solver is employed to march in time the equations with a time step of $10^{-6}$ and the total number of modes used was 100. Note the large size of the aeroelastic ODE system: 2 \times 100 nonlinear equations plus 5 \times 100 linear equations for the aerodynamic states with 5 poles, plus 4 equations for the quaternion tracking the rigid-body motion, for a combined ODE system of 704 equations. +In addition, a total of 512 gusts cases are run concurrently for all possible combinations of 8 gust lengths between 25 and 265 meters, 8 gust intensities between 1 and 30 m/s, and 8 flow densities between 0.34 and 0.48 Kg/m$^3$. This means that $512 \times 704 = 360448$ equations are being marched in time, in this case for 2 seconds which is enough to capture peak loads and with a time step of 0.001 seconds. +Table [[table:times_gust]] contains the simulation times of the calculation, which show one order of magnitude increase in performance when running in parallel in the CPU versus a complete single simulation running sequentially, and another order of magnitude when moving from the CPU to a modern GPU. -#+NAME: table:FFB_times -#+CAPTION: FFB computational times and Cg error +#+CAPTION: Computational times multiple gust problem #+ATTR_LATEX: :center t -| Arch/Nmodes | 50M | 100M | 150M | 200M | 250M | 300M | -|---------------------+------------+-------------+------------+------------+------------+--------------| -| CPU HPC (time/err) | 7/1.3e-1 | 9.3/5.7e-2 | 34/2.2e-2 | 79/2e-3 | 474/5.3e-4 | 1869/2.54e-5 | -| GPU HPC (var. dt) | 9.9/1.3e-1 | 10.4/5.7e-2 | 14/2.2e-2 | 22/2e-3 | 38/5.3e-4 | 111/2.54e-5 | -|---------------------+------------+-------------+------------+------------+------------+--------------| -| CPU HPC (time/err) | 42/2.1e-2 | 184/1.2e-2 | 287/5.6e-3 | 421/7.2e-4 | 893/2.7e-4 | 1869/2.54e-5 | -| GPU HPC (const. dt) | 58/2.1e-2 | 65/1.2e-2 | 67/5.6e-3 | 76/7.2e-4 | 94/2.7e-4 | 111/2.54e-5 | -|---------------------+------------+-------------+------------+------------+------------+--------------| - -Fig. [[fig:FFBtimes2]] and [[fig:FFBerror2]] illustrate the times and error results in the table for the second case with fixed time step. The gain in performance from the GPU is is more impressive the larger the system to solve, and for the full modal basis the CPU takes more than 31 minutes versus the less than 2 minutes in the GPU. Computational times in the 3D problem are similar and the error on the Cg position is slightly higher: for the 300 modes, the error is $6.9e-5$ versus the $2.54e-5$ of the planar case. - -#+NAME: fig:FFBtimes2 -#+CAPTION: Performance CPU vs GPU comparison in free-flying structure -#+ATTR_LATEX: :width 0.5\textwidth -#+RESULTS: FFBtimes2 -[[file:figs/FFBtimes2.png]] - -#+NAME: fig:FFBerror2 -#+CAPTION: Error metric Cg position for planar case -#+ATTR_LATEX: :width 0.5\textwidth -#+RESULTS: FFBerror2 -[[file:figs/FFBerror2.png]] - - -* Further work -The next step is to characterise the geometrically nonlinear aeroelastic response of the BUG High-Aspect-Ratio aircraft [[cite:&STODIECK2018a]] shown in Fig. [[fig:bug]]. We aim to achieve the following points: - -- Fine-tune the aeroelastic solvers to introduce an updating mechanism of the normal of the aerodynamic panels to account for the nonlinear effect around deformed configurations such as trimmed flight. -- Build manoeuvre and dynamic load envelopes that can also be differentiated via AD. -- Bring parallelisation on distributed accelerators into our solution process, thus applying the state-of-the-art techniques in used for the large problems in deep-learning. Combined with our already fast simulations times, this could allow the prediction of those sizing aeroelastic loads that include thousands of cases in commercial aircraft, the computation of their gradients with respect to design variables, with geometrically nonlinear effects accounted for, and at such performance that the framework could be integrated into a larger multidisciplinary optimization. - -#+NAME: fig:bug -#+CAPTION: BUG static equilibrium -#+ATTR_LATEX: :width 0.8\textwidth -[[file:figs_ext/bug_model.pdf]] +#+NAME: table:times_gust +| Device | Time (sec.) | +|----------------+---------------------------| +| CPU (single) | 27.8 \times 512 = 14233.6 | +| CPU (parallel) | 922.6 | +| GPU | 38.2 | + + +Fig. [[fig:BUG_Gust3D]] shows the 3D flight shape of the airframe for a gust of 150 m length, intensity of 20 m/s and flow density of 0.41 Kg/m$^3$. +#+NAME: fig:BUG_Gust3D +#+CAPTION: Dynamic response to 1-cos gust excitation +#+ATTR_LATEX: :width 1\textwidth +[[file:figs_ext/bug_gust3d.pdf]] + +# #+NAME: fig:BUG_Gust_envelope +# #+CAPTION: Load envelope for maximum loads in the dynamic response to 1-cos gust excitation with +# #+ATTR_LATEX: :width 1\textwidth +# [[file:figs_ext/]] + +* Conclusions +A modal-based, geometrically nonlinear formulation of the aircraft dynamics has been enhanced with multiple load-cases parallelisation in modern hardware architectures. We have applied state-of-the-art techniques and tools employed for large problems in Deep Machine Learning to the computation and prediction of the sizing aeroelastic loads in commercial aircraft, which can expand thousands of simulations. +Remarkable computational times of under a minute are achieved for 256 manoeuvres varying flow conditions and for 512 dynamic gust responses, including geometrically nonlinear effects in the simulations. +Such a performance potentially unlocks two different applications: uncertainty quantification of the nonlinear aircraft response to a non-deterministic loading and integration of the software in larger multidisciplinary frameworks. +The former has been demonstrated on a problem where a field of forces with an stochastic component induces very large deformation; it has been shown that while the statistics in the linear response can be easily forecast from one complete experiment, in the nonlinear case a Montecarlo simulation needs to be carried out for each new set of loading scenario. +For the latter, differentiation of the load envelopes via the AD capabilities within JAX will be the next step. Since we are already in the memory boundaries of a single GPU or CPU, this will require the use of multiple devices, for which we have already built the implementation. +Scaling up the process to include various mass cases, as it is done in industrial scenarios, is also a feasible target. Thus combining prediction of sizing aeroelastic loads that include thousands of cases in commercial aircraft with the computation of their gradients with respect to design variables in a framework for multidisciplinary design optimization. bibliographystyle:unsrt # bibliography:/home/acea/Documents/Engineering.bib diff --git a/docs/schemas/draft1.org b/docs/schemas/draft1.org index e8958dc..8da8c48 100644 --- a/docs/schemas/draft1.org +++ b/docs/schemas/draft1.org @@ -1,4 +1,4 @@ -#+TITLE: FENIAX.initial architecture +#+TITLE: FENIAX architecture #+AUTHOR: Alvaro Cea #+Date: 6/06/2023 #+STARTUP: inlineimages @@ -7,7 +7,7 @@ (setq org-image-actual-width nil) (setq org-confirm-babel-evaluate nil) (require 'org-tempo) - (pyvenv-workon "feniax.) + (pyvenv-workon "feniax") #+end_src diff --git a/docs/schemas/draft2.org b/docs/schemas/draft2.org new file mode 100644 index 0000000..c90c6dd --- /dev/null +++ b/docs/schemas/draft2.org @@ -0,0 +1,410 @@ +#+TITLE: FENIAX architecture +#+AUTHOR: Alvaro Cea +#+Date: 11/29/2024 +#+STARTUP: inlineimages +#+LATEX_HEADER: \usepackage{algpseudocode} +#+LATEX_HEADER: \usepackage[ruled,vlined]{algorithm2e} + +* HouseKeeping +#+begin_src elisp :results none :tangle no + (setq org-image-actual-width nil) + (setq org-confirm-babel-evaluate nil) + ;; (setq org-latex-pdf-process + ;; '("latexmk -pdflatex='pdflatex --syntex=1 -interaction nonstopmode' -pdf -bibtex -f %f")) + (setq org-latex-pdf-process + '("latexmk -pdflatex='pdflatex --syntex=1 -shell-escape -interaction nonstopmode' -pdf -bibtex -f %f" + "latexmk -pdflatex='pdflatex --syntex=1 -shell-escape -interaction nonstopmode' -pdf -bibtex -f %f")) + + (require 'org-tempo) + (pyvenv-workon "feniax") + +#+end_src + +* Challenging problems driving the research :noexport: +** Geometrically nonlinear aeroelastic analysis using complex GFEMs +- Flight loads +- Airframe integrity: flutter or divergence assessment +- Multidisciplinary design optimisation +** Max performance: +- Critical for time-domain simulations +- JAX Just in Time Compilation (JIT) +- Code vectorisation and parallelisation. +- Run on accelerators (GPUs, TPUs, etc.) +** Derivatives of aeroelastic for design optimisation. +- Using JAX algorithmic differentiation. +*** Efficient derivatives +JAX AD package using pure functions from functional programming style. +*** Compute and arbitrary number of loads cases +Via subcase option of driver class that modifies the xloads option in each case. For example to get the max. loads envelope without having to manage an array of simulations which might also overlap computations. +*** Compute the derivative of function applied to the previous point +If for instance one wants to set an optimisation with boundaries on the maximum loads, it would not make sense to do it for only one atmospheric load case. +*** Perform calculations on different models in a single run +Via the supercase option of driver. This might be used to compute finite differences; or if the response of if the response of a fractured component is to be compared. + +* Code design and software architecture :noexport: +- Config object with input settings. +- Driver class: initialises all relevant objects such as the simulation, and the systems that will solved in the solution process. +- Simulation class: responsible for running the various systems appropriately, including setting the initial conditions from one system to another. +- System class: set the computations to solve the corresponding system of equations, including the solver library that should be called, the system of equations and the arguments to the solvers. + +** Config +config.engine = intrinsic +config.supercase.fems[dict] +config.subcase.system[].xloads +config.simulation.typeof[Serial, parallel, Single] +config.simulation.build_grads +config.simulation.optimize +config. +config.systems[] +config.system.name +config.system.solver.library +config.system.solver.settings +** Driver +(only driver gets to modify object) +run_cases + -- set_case + -- integration.pre_simulation() + -- simulation.trigger() + -- +set_case +(modifies config object) + -- Supercase + -- Subcase +set_integration +set_simulation -> simulation +** Integration +run + -- calculate_modalshapes + -- calculate_modalcouplings +derivatives + +** Simulation +- trigger +- _prerun +- _run + -- system.set_init +- _pull_solution +- _postrun +*** SerialSimulation +-_run +*** ParallelSimulation +-_run +*** SingleSimulation +- _run +** System +- set_init -> q0 +- set_name +- set_generator -> dq +- set_solver +- solve -> q + self.solver(self.dq) +- save + + +Static and dynamic systems for static and dynamic simulations + +Systems with labels: + +[[file:~/projects/FENIAX/feniax/systems/intrinsic_system.py::label = f"dq_{self.settings.label}"][intrinsic_system]] +[[file:~/projects/FENIAX/feniax/systems/intrinsicAD.py::label = f"main_{label_sys}_{label_ad}"][AD_system]] +[[file:~/projects/FENIAX/feniax/systems/intrinsicShard.py::self.label = f"main_{label_sys}_{label_shard}"][shard_system]] + + +** XForces +*** prescribed_follower +*** prescribed_dead +*** gravity +*** modal_aero +** inputs +*** container +*** fields +- value +- description +- default +- options + + +** UML + +#+Name: classes_architecture +#+begin_src plantuml :file UML_software1.png + abstract Driver { + +pre_simulation() + +run_cases() + } + + class IntrinsicDriver { + #integration: IntrinsicIntegration + #simulation: Simulation + #opt: Optimisation + #systems: [System] + -__init__(config: Config) + #_set_case() + #_set_integration() + #_set_simulation() + #_set_systems() + } + + class XLoads { + +q: [jnp.ndarray] + +Rab: [jnp.ndarray] + +GAFs: [jnp.ndarray] + -__init__(config.systems.loads, + q, Rab, GAFs) + +followerF() + +deadF() + +gravityF() + +modalAero() + } + + /' + ' abstract class Integration { + ' +run() + ' } + '/ + + class IntrinsicIntegration { + + phi_1, phi_2, psi_1, psi_2 + + Gamma_1, Gamma_2 + -__init__(X, Ka, Ma) + +run() + #compute_modalshapes() + #compute_modalcouplings() + } + + abstract class Simulation { + +systems: [System] + #workflow: dict[str:str] + #opt: Optimisation + -__init__(config.simulation, + systems, opt, config.simulation) + +trigger() + #run_systems() + #post_run() + } + + /' + ' package Simulations { + ' class SerialSimulation { + ' } + ' class ParallelSimulation { + ' } + ' class SingleSimulation { + ' } + ' class CoupledSimulation { + ' } + ' } + '/ + + class SerialSimulation { + } + class ParallelSimulation { + } + class SingleSimulation { + } + class CoupledSimulation { + } + + abstract class System { + +set_ic(q0) + +solve() -> sol + +pull_solution() -> qs + } + + class IntrinsicSystem { + -__init__(name[str], settings:config.Dsystem, + fem: config.Dfem, + sol: solution.IntrinsicSolution) + -dq: callable + -solver: callable + +sol: obj + #set_generator() -> dq + #set_solver() -> solver + + } + + class ControlSystem { + } + + class MultibodySystem { + } + + /' + ' Simulation <|-- SerialSimulation + ' Simulation <|-- ParallelSimulation + ' Simulation <|-- SingleSimulation + ' Simulation <|-- CoupledSimulation + '/ + abstract class Optimisation { + +save_grads() + +assemble() + } + + abstract class Sollibs { + +name() + +pull_name() + } + enum dq { + - sol_dict + - dq_label + } + + enum loads { + - eta_dict[] + - eta_{label} + } + + 'Simulation <|-- Simulations + Simulation <|-- SingleSimulation + SingleSimulation -- SerialSimulation + SerialSimulation -- ParallelSimulation + ParallelSimulation -- CoupledSimulation + 'Driver "1" -- "1" Integration : composition + 'Driver "1" -- "1" Simulation : composition' + IntrinsicIntegration -* IntrinsicDriver + Driver <|-- IntrinsicDriver + IntrinsicDriver *-- Optimisation + IntrinsicDriver *-- Simulation + IntrinsicDriver *- System + System ..> Simulation + 'Integration <|-- IntrinsicIntegration + System <|-- IntrinsicSystem + System *- Sollibs + IntrinsicSystem -- ControlSystem + IntrinsicSystem o- XLoads + ControlSystem -- MultibodySystem +#+end_src + +#+RESULTS: classes_architecture +[[file:UML_software1.png]] + +* Algorithms + + +$N_N$ Number of discretisation condensed nodes in the model +$N_m$ Number of modal shapes in the solution +$N_t$ Number of time-steps in the integration scheme. +$N_c$ Number of cases to be run in parallel + + +#+NAME: alg:process +\begin{algorithm}[h!] +\DontPrintSemicolon +\SetKwInOut{Input}{input} +\SetKwInOut{Output}{output} +\Input{Input file: settings.yaml; FE model: $\bm{K}_a$, $\bm{M}_a$, $\bm{X}_a$; Aerodynamic matrices: $\bm{\mathcal{A}}$} +\Output{Nonlinear aeroealastic solutioxn} +\Begin{ + \BlankLine +$\bm{\phi}$, $\bm{\psi}$ $\longleftarrow$ modes($\bm{K}_a$, $\bm{M}_a$, $\bm{X}_a$) \Comment{Intrinsic modes: O($N_n^2 \times N_m$; $N_n \times N_m$)} \; +$\bm{\Gamma}$ $\longleftarrow$ couplings($\bm{\phi}$, $\bm{\psi}$) \Comment{Nonlinear couplings O($N_n \times N_m^3$; $N_m^3$)} \; +$\bm{q}$ $\longleftarrow$ system($\bm{\Gamma}$, $\bm{\mathcal{A}}$, $\bm{\phi}$, $\bm{X}_a$) \Comment{Modal coordinates: O($\frac{N_l}{N_d} \times N_t \times N_m^3$; $N_l \times N_t \times N_m$)} \; +$\bm{X}_1$, $\bm{X}_{2}$, $\bm{X}_{3}$ $\longleftarrow$ ivars($\bm{q}$, $\bm{\phi}$, $\bm{\psi}$) \Comment{velocity/strain fields: O($\frac{N_l}{N_d} \times N_t \times N_n \times N_m$; $N_l \times N_t \times N_n$)} \; +$\bm{r}_a$, $\bm{R}_{a}$ $\longleftarrow$ integration($\bm{X}_{3}$, $\bm{X}_a$) \Comment{Positional/rotational fields: O($\frac{N_l}{N_d} \times N_t \times N_n \times N_m$; $N_l \times N_t \times N_n$)} \; +\BlankLine +} +\caption{Main components in solution process} +\end{algorithm} + + +* Simulation inputs :noexport: +trim +---- +qh = 0 +qe becomes unknown +qalpha != 0 +qhdot = f(gamma2) + eta_h(q0, qe) = 0 # rigid bodies () +qalphadot = f(gamma2, q2) + eta_alpha(q0, qe) # rigid bodies +q1dot = f(gamma2) + eta_q(q0, qe) + + + Connection with High Fidelity structural model +** Initial Model +- Clamped wing -> good for steady loads and aircraft stability. +- Full A/C model -> Needed for dynamic loads. +- Mass model: Both continuous mass model and lumped masses are suitable for analysis. + Engines and other components definitely as lumped masses. +** Input requirements + +#+ATTR_ORG: :width 600 +[[./FEM3d.png]] +*** Load paths +- interpolation elements to connect to other FE nodes. +- aerodynamic forces applied along these paths + +*** Condensed stiffness and mass matrices along load paths +- Should be suitable for eigenvalue analysis + +*** Aerodynamic model via GAFs. +- Preliminary DLM model. Automatically built from wing-box? +- Steady loads: Corrections may be needed. + +** Output requirements +*** Sectional loads along load-paths +Steady and dynamic aeroelastic loads due to trimmed flight, gusts etc. +*** Aeroelastic stability of configuration +Flutter and divergence points +*** Potentially derivatives of the above via AD. +Critical for large design optimisation problems. +** Data workflow +Well in place for Nastran Models except for the derivatives provided by Nastran using Sol 200. + + + +* System based solutions :noexport: +TODO: make automatic label as the first +| Type | Target | Gravity | BC1 | ModalAero | SteadyAero | UnsteadyAero | Point loads | q0 approx | Rigid-body | Nonlinearities | residualised | +|--------------+--------+------------+------------+-----------+--------------+--------------+-------------+-----------+----------------------+------------------------+--------------| +| 1 static | Level | False: "g" | Clamped | None | None | None | None | via q2 | 1-quaternion+strains | All -> "" | None -> "" | +| 2 Dynamic | TRIM1 | True: "G" | Free | Rogers | qalpha | gust | follower | via q1 | All-quaternions | Linear sys -> "l" | True -> "r" | +| 3 staticAD | | | | | | | | | | | | +| 4 dynamicAD | | | | | | | | | | | | +| 3 staticPL | TRIM2 | | Prescribed | Loewner | qx (control) | controls | dead | | | Linear sys+disp -> "L" | | +| 3 dynamicPL | TRIM2 | | Prescribed | Loewner | qx (control) | controls | dead | | | Linear sys+disp -> "L" | | +| 3 staticPLAD | TRIM2 | | Prescribed | Loewner | qx (control) | controls | dead | | | Linear sys+disp -> "L" | | +| 3 | TRIM2 | | Prescribed | Loewner | qx (control) | controls | dead | | | Linear sys+disp -> "L" | | + +| 3 Stability | TRIM2 | | Prescribed | Loewner | qx (control) | controls | dead | | | Linear sys+disp -> "L" | | +| 4 Multibody | | | | | | | | | | | | +| 5 Control | | | | | | | | | | | | + +| Sol name | | label | Imp | +|----------+-------------------------------------------------+-----------------------+-----| +| 10G1 | Structural static under Gravity | [1,0,'G'] | Y | +| 10g11 | Structural static with follower point forces | [1,0,'g',0,0,0,0,1] | Y | +| 10g121 | Structural static with dead point forces | [1,0,'g',0,0,0,0,2] | Y | +| 10g1331 | Structural static with follower+dead forces | [1,0,'g',0,0,0,0,3] | N | +| 10g15 | Manoeuvre under qalpha | [1,0,'g',0,1,1] | Y | +| 10G15 | Manoeuvre under qalpha and Gravity | [1,0,'G',0,1,1] | N | +| 10g75 | Manoeuvre under qalpha and controls | [1,0,'g',0,1,2] | N | +| 10G75 | Manoeuvre under qalpha+controls+Gravity | [1,0,'G',0,1,2] | N | +| 20g1 | Clamped Structural dynamics, free vibrations | [2,0,'g'] | Y | +| 20G2 | Free Structural dynamic with gravity forces | [2,0,'G',1] | Y | +| 20g11 | Structural dynamic follower point forces | [2,0,'g',0,0,0,0,1] | Y | +| 20g121 | Structural dynamic dead point forces | [2,0,'g',0,0,0,0,2] | Y | +| 20g22 | Free Structural dynamic follower point forces | [2,0,'g',1,0,0,0,1] | Y | +| 20g242 | Free Structural dynamic dead point forces | [2,0,'g',1,0,0,0,2] | Y | +| 11G6 | Static trimmed State (elevator-qalpha, | [1,1,'G',1,1] | Y | +| | no gravity updating) | | | +| 12G2 | Static trimmed State (elevator-qalpha, | [1,2,'G',1] | N | +| | gravity updating) | | | +| 21G150 | Dynamic trimmed State | [2,1,'G',1,1,2] | N | +| 20g21 | Gust response | [2,0,'g',0,1,0,1] | Y | +| 20g273 | Gust response, q0 obtained via integrator q1 | [2,0,'g',0,1,0,1,0,1] | Y | +| 20g105 | Gust response with steady qalpha | [2,0,'g',0,1,1,1] | N | +| 20g42 | Gust response Free-flight | [2,0,'g',1,1,0,1] | N | +| 20G42 | Gust response Free-flight and gravity (X error) | [2,0,'G',1,1,0,1] | N | +| 20G1050 | Gust response Free-flight, gravity, controls | [2,0,'G',1,1,2,1] | N | +| | | | | + +#+begin_src python :session py1 :results output + import feniax.intrinsic.functions as functions + label = functions.label_generator([2,0,'g',0,1,0,1,0,1]) + print(label) +#+end_src + +#+RESULTS: +: 20g546 +:  + + diff --git a/examples/BUG/AERO/ADd1c7F3Scao-100p5.npy b/examples/BUG/AERO/ADd1c7F3Scao-100p5.npy new file mode 100644 index 0000000..833c8ef Binary files /dev/null and b/examples/BUG/AERO/ADd1c7F3Scao-100p5.npy differ diff --git a/examples/BUG/AERO/ADd1c7F3Scao-50p5.npy b/examples/BUG/AERO/ADd1c7F3Scao-50p5.npy new file mode 100644 index 0000000..f4dbea8 Binary files /dev/null and b/examples/BUG/AERO/ADd1c7F3Scao-50p5.npy differ diff --git a/examples/BUG/AERO/ADd1c7F3Seao-100p5.npy b/examples/BUG/AERO/ADd1c7F3Seao-100p5.npy new file mode 100644 index 0000000..76611c2 Binary files /dev/null and b/examples/BUG/AERO/ADd1c7F3Seao-100p5.npy differ diff --git a/examples/BUG/AERO/ADd1c7F3Seao-50p5.npy b/examples/BUG/AERO/ADd1c7F3Seao-50p5.npy new file mode 100644 index 0000000..93ca826 Binary files /dev/null and b/examples/BUG/AERO/ADd1c7F3Seao-50p5.npy differ diff --git a/examples/BUG/AERO/Collocation_d1c7.npy b/examples/BUG/AERO/Collocation_d1c7.npy new file mode 100644 index 0000000..8da3c9b Binary files /dev/null and b/examples/BUG/AERO/Collocation_d1c7.npy differ diff --git a/examples/BUG/AERO/DDd1c7F3Scao-100p5.npy b/examples/BUG/AERO/DDd1c7F3Scao-100p5.npy new file mode 100644 index 0000000..95aefa9 Binary files /dev/null and b/examples/BUG/AERO/DDd1c7F3Scao-100p5.npy differ diff --git a/examples/BUG/AERO/DDd1c7F3Scao-50p5.npy b/examples/BUG/AERO/DDd1c7F3Scao-50p5.npy new file mode 100644 index 0000000..b0916f7 Binary files /dev/null and b/examples/BUG/AERO/DDd1c7F3Scao-50p5.npy differ diff --git a/examples/BUG/AERO/DDd1c7F3Seao-100p5.npy b/examples/BUG/AERO/DDd1c7F3Seao-100p5.npy new file mode 100644 index 0000000..7708298 Binary files /dev/null and b/examples/BUG/AERO/DDd1c7F3Seao-100p5.npy differ diff --git a/examples/BUG/AERO/DDd1c7F3Seao-50p5.npy b/examples/BUG/AERO/DDd1c7F3Seao-50p5.npy new file mode 100644 index 0000000..e6e4f9b Binary files /dev/null and b/examples/BUG/AERO/DDd1c7F3Seao-50p5.npy differ diff --git a/examples/BUG/AERO/Dihedral_d1c7.npy b/examples/BUG/AERO/Dihedral_d1c7.npy new file mode 100644 index 0000000..6069c22 Binary files /dev/null and b/examples/BUG/AERO/Dihedral_d1c7.npy differ diff --git a/examples/BUG/AERO/PolesDd1c7F1Scao-50p10.npy b/examples/BUG/AERO/PolesDd1c7F1Scao-50p10.npy new file mode 100644 index 0000000..f17ca63 Binary files /dev/null and b/examples/BUG/AERO/PolesDd1c7F1Scao-50p10.npy differ diff --git a/examples/BUG/AERO/PolesDd1c7F1Scao-50p5.npy b/examples/BUG/AERO/PolesDd1c7F1Scao-50p5.npy new file mode 100644 index 0000000..0f9bd87 Binary files /dev/null and b/examples/BUG/AERO/PolesDd1c7F1Scao-50p5.npy differ diff --git a/examples/BUG/AERO/PolesDd1c7F1Seao-50p10.npy b/examples/BUG/AERO/PolesDd1c7F1Seao-50p10.npy new file mode 100644 index 0000000..f17ca63 Binary files /dev/null and b/examples/BUG/AERO/PolesDd1c7F1Seao-50p10.npy differ diff --git a/examples/BUG/AERO/PolesDd1c7F1Seao-50p5.npy b/examples/BUG/AERO/PolesDd1c7F1Seao-50p5.npy new file mode 100644 index 0000000..ae8d1ef Binary files /dev/null and b/examples/BUG/AERO/PolesDd1c7F1Seao-50p5.npy differ diff --git a/examples/BUG/AERO/PolesDd1c7F2Scao-50p10.npy b/examples/BUG/AERO/PolesDd1c7F2Scao-50p10.npy new file mode 100644 index 0000000..f17ca63 Binary files /dev/null and b/examples/BUG/AERO/PolesDd1c7F2Scao-50p10.npy differ diff --git a/examples/BUG/AERO/PolesDd1c7F2Scao-50p5.npy b/examples/BUG/AERO/PolesDd1c7F2Scao-50p5.npy new file mode 100644 index 0000000..fcccfe1 Binary files /dev/null and b/examples/BUG/AERO/PolesDd1c7F2Scao-50p5.npy differ diff --git a/examples/BUG/AERO/PolesDd1c7F2Seao-50p5.npy b/examples/BUG/AERO/PolesDd1c7F2Seao-50p5.npy new file mode 100644 index 0000000..ae8d1ef Binary files /dev/null and b/examples/BUG/AERO/PolesDd1c7F2Seao-50p5.npy differ diff --git a/examples/BUG/AERO/PolesDd1c7F3Scao-100p5.npy b/examples/BUG/AERO/PolesDd1c7F3Scao-100p5.npy new file mode 100644 index 0000000..fcccfe1 Binary files /dev/null and b/examples/BUG/AERO/PolesDd1c7F3Scao-100p5.npy differ diff --git a/examples/BUG/AERO/PolesDd1c7F3Scao-50p5.npy b/examples/BUG/AERO/PolesDd1c7F3Scao-50p5.npy new file mode 100644 index 0000000..fcccfe1 Binary files /dev/null and b/examples/BUG/AERO/PolesDd1c7F3Scao-50p5.npy differ diff --git a/examples/BUG/AERO/PolesDd1c7F3Seao-100p5.npy b/examples/BUG/AERO/PolesDd1c7F3Seao-100p5.npy new file mode 100644 index 0000000..03493c1 Binary files /dev/null and b/examples/BUG/AERO/PolesDd1c7F3Seao-100p5.npy differ diff --git a/examples/BUG/AERO/PolesDd1c7F3Seao-50p5.npy b/examples/BUG/AERO/PolesDd1c7F3Seao-50p5.npy new file mode 100644 index 0000000..ae8d1ef Binary files /dev/null and b/examples/BUG/AERO/PolesDd1c7F3Seao-50p5.npy differ diff --git a/examples/BUG/AERO/QahDd1c7F2Scao-50.npy b/examples/BUG/AERO/QahDd1c7F2Scao-50.npy new file mode 100644 index 0000000..b9ddbe0 Binary files /dev/null and b/examples/BUG/AERO/QahDd1c7F2Scao-50.npy differ diff --git a/examples/BUG/AERO/QahDd1c7F3Scao-100.npy b/examples/BUG/AERO/QahDd1c7F3Scao-100.npy new file mode 100644 index 0000000..347293a Binary files /dev/null and b/examples/BUG/AERO/QahDd1c7F3Scao-100.npy differ diff --git a/examples/BUG/AERO/QahDd1c7F3Scao-50.npy b/examples/BUG/AERO/QahDd1c7F3Scao-50.npy new file mode 100644 index 0000000..ebed8cb Binary files /dev/null and b/examples/BUG/AERO/QahDd1c7F3Scao-50.npy differ diff --git a/examples/BUG/AERO/QaxDd1c7F2Scao-50.npy b/examples/BUG/AERO/QaxDd1c7F2Scao-50.npy new file mode 100644 index 0000000..877ee6a Binary files /dev/null and b/examples/BUG/AERO/QaxDd1c7F2Scao-50.npy differ diff --git a/examples/BUG/AERO/QaxDd1c7F3Scao-100.npy b/examples/BUG/AERO/QaxDd1c7F3Scao-100.npy new file mode 100644 index 0000000..6775257 Binary files /dev/null and b/examples/BUG/AERO/QaxDd1c7F3Scao-100.npy differ diff --git a/examples/BUG/AERO/QaxDd1c7F3Scao-50.npy b/examples/BUG/AERO/QaxDd1c7F3Scao-50.npy new file mode 100644 index 0000000..6775257 Binary files /dev/null and b/examples/BUG/AERO/QaxDd1c7F3Scao-50.npy differ diff --git a/examples/BUG/AERO/QhxDd1c7F2Scao-50.npy b/examples/BUG/AERO/QhxDd1c7F2Scao-50.npy new file mode 100644 index 0000000..237e80c Binary files /dev/null and b/examples/BUG/AERO/QhxDd1c7F2Scao-50.npy differ diff --git a/examples/BUG/AERO/QhxDd1c7F3Scao-100.npy b/examples/BUG/AERO/QhxDd1c7F3Scao-100.npy new file mode 100644 index 0000000..4d3aade Binary files /dev/null and b/examples/BUG/AERO/QhxDd1c7F3Scao-100.npy differ diff --git a/examples/BUG/AERO/QhxDd1c7F3Scao-50.npy b/examples/BUG/AERO/QhxDd1c7F3Scao-50.npy new file mode 100644 index 0000000..dfcc7f7 Binary files /dev/null and b/examples/BUG/AERO/QhxDd1c7F3Scao-50.npy differ diff --git a/examples/BUG/FEM/Ka.npy b/examples/BUG/FEM/Ka_ca.npy similarity index 100% rename from examples/BUG/FEM/Ka.npy rename to examples/BUG/FEM/Ka_ca.npy diff --git a/examples/BUG/FEM/Ka_ea.npy b/examples/BUG/FEM/Ka_ea.npy new file mode 100644 index 0000000..de4b6eb Binary files /dev/null and b/examples/BUG/FEM/Ka_ea.npy differ diff --git a/examples/BUG/FEM/Ma.npy b/examples/BUG/FEM/Ma_ca.npy similarity index 100% rename from examples/BUG/FEM/Ma.npy rename to examples/BUG/FEM/Ma_ca.npy diff --git a/examples/BUG/FEM/Ma_ea.npy b/examples/BUG/FEM/Ma_ea.npy new file mode 100644 index 0000000..08eb4ca Binary files /dev/null and b/examples/BUG/FEM/Ma_ea.npy differ diff --git a/examples/BUG/FEM/eigenvals_50.npy b/examples/BUG/FEM/eigenvals_cao100.npy similarity index 100% rename from examples/BUG/FEM/eigenvals_50.npy rename to examples/BUG/FEM/eigenvals_cao100.npy diff --git a/examples/BUG/FEM/eigenvals_cao300.npy b/examples/BUG/FEM/eigenvals_cao300.npy new file mode 100644 index 0000000..cfe9333 Binary files /dev/null and b/examples/BUG/FEM/eigenvals_cao300.npy differ diff --git a/examples/BUG/FEM/eigenvals_cao50.npy b/examples/BUG/FEM/eigenvals_cao50.npy new file mode 100644 index 0000000..3c7116d Binary files /dev/null and b/examples/BUG/FEM/eigenvals_cao50.npy differ diff --git a/examples/BUG/FEM/eigenvals_eao100.npy b/examples/BUG/FEM/eigenvals_eao100.npy new file mode 100644 index 0000000..9fe4db6 Binary files /dev/null and b/examples/BUG/FEM/eigenvals_eao100.npy differ diff --git a/examples/BUG/FEM/eigenvals_eao50.npy b/examples/BUG/FEM/eigenvals_eao50.npy new file mode 100644 index 0000000..5a22ed2 Binary files /dev/null and b/examples/BUG/FEM/eigenvals_eao50.npy differ diff --git a/examples/BUG/FEM/eigenvecs_50.npy b/examples/BUG/FEM/eigenvecs_50.npy deleted file mode 100644 index 9e3f9a1..0000000 Binary files a/examples/BUG/FEM/eigenvecs_50.npy and /dev/null differ diff --git a/examples/BUG/FEM/eigenvecs_cao100.npy b/examples/BUG/FEM/eigenvecs_cao100.npy new file mode 100644 index 0000000..7523e3e Binary files /dev/null and b/examples/BUG/FEM/eigenvecs_cao100.npy differ diff --git a/examples/BUG/FEM/eigenvecs_cao300.npy b/examples/BUG/FEM/eigenvecs_cao300.npy new file mode 100644 index 0000000..4b37199 Binary files /dev/null and b/examples/BUG/FEM/eigenvecs_cao300.npy differ diff --git a/examples/BUG/FEM/eigenvecs_cao50.npy b/examples/BUG/FEM/eigenvecs_cao50.npy new file mode 100644 index 0000000..21fe805 Binary files /dev/null and b/examples/BUG/FEM/eigenvecs_cao50.npy differ diff --git a/examples/BUG/FEM/eigenvecs_eao100.npy b/examples/BUG/FEM/eigenvecs_eao100.npy new file mode 100644 index 0000000..8338403 Binary files /dev/null and b/examples/BUG/FEM/eigenvecs_eao100.npy differ diff --git a/examples/BUG/FEM/eigenvecs_eao50.npy b/examples/BUG/FEM/eigenvecs_eao50.npy new file mode 100644 index 0000000..8bfca61 Binary files /dev/null and b/examples/BUG/FEM/eigenvecs_eao50.npy differ diff --git a/examples/BUG/FEM/structuralGrid b/examples/BUG/FEM/structuralGrid_ca similarity index 100% rename from examples/BUG/FEM/structuralGrid rename to examples/BUG/FEM/structuralGrid_ca diff --git a/examples/BUG/FEM/structuralGrid_ea b/examples/BUG/FEM/structuralGrid_ea new file mode 100644 index 0000000..47a53a7 --- /dev/null +++ b/examples/BUG/FEM/structuralGrid_ea @@ -0,0 +1,100 @@ +# VARIABLES = x y z fe_order components +16.5335 0.0 0.0 5 FusBack +22.8143 0.0 0.0 6 FusBack +27.4511 0.0 0.0 7 FusBack +30.419 0.0 0.2989 8 FusBack +34.1749 0.0 0.6771 9 FusBack +15.6286 0.0 0.0 4 FusFront +11.4071 0.0 0.0 3 FusFront +7.6048 0.0 0.0 2 FusFront +3.8062 0.0 0.0 1 FusFront +0.0 0.0 -0.6028 0 FusFront +16.5335 1.00156 2.11455 11 RWing +16.6244 2.45678 2.10285 12 RWing +16.8062 3.36411 2.07944 13 RWing +17.079 4.7251 2.04434 14 RWing +17.2676 5.6656 2.02122 15 RWing +17.4629 6.63927 1.99839 16 RWing +17.6582 7.61294 1.97556 17 RWing +17.8535 8.58661 1.95273 18 RWing +18.0488 9.56028 1.92989 19 RWing +18.2441 10.5339 1.90706 20 RWing +18.4394 11.5076 1.88423 21 RWing +18.6347 12.4813 1.8614 22 RWing +18.83 13.455 1.83857 23 RWing +19.0253 14.4286 1.81574 24 RWing +19.2359 15.3927 1.79321 25 RWing +19.4617 16.3472 1.771 26 RWing +19.6876 17.3018 1.74878 27 RWing +19.9134 18.2563 1.72656 28 RWing +20.1393 19.2108 1.70434 29 RWing +20.3651 20.1653 1.68212 30 RWing +20.591 21.1198 1.65991 31 RWing +20.8168 22.0744 1.63769 32 RWing +21.0427 23.0289 1.61547 33 RWing +21.2685 23.9834 1.59325 34 RWing +21.4944 24.9379 1.57104 35 RWing +21.7202 25.8924 1.54882 36 RWing +16.5335 -1.00156 2.11455 60 LWing +16.6244 -2.45678 2.10285 61 LWing +16.8062 -3.36411 2.07944 62 LWing +17.079 -4.7251 2.04434 63 LWing +17.2676 -5.6656 2.02122 64 LWing +17.4629 -6.63927 1.99839 65 LWing +17.6582 -7.61294 1.97556 66 LWing +17.8535 -8.58661 1.95273 67 LWing +18.0488 -9.56028 1.92989 68 LWing +18.2441 -10.5339 1.90706 69 LWing +18.4394 -11.5076 1.88423 70 LWing +18.6347 -12.4813 1.8614 71 LWing +18.83 -13.455 1.83857 72 LWing +19.0253 -14.4286 1.81574 73 LWing +19.2359 -15.3927 1.79321 74 LWing +19.4617 -16.3472 1.771 75 LWing +19.6876 -17.3018 1.74878 76 LWing +19.9134 -18.2563 1.72656 77 LWing +20.1393 -19.2108 1.70434 78 LWing +20.3651 -20.1653 1.68212 79 LWing +20.591 -21.1198 1.65991 80 LWing +20.8168 -22.0744 1.63769 81 LWing +21.0427 -23.0289 1.61547 82 LWing +21.2685 -23.9834 1.59325 83 LWing +21.4944 -24.9379 1.57104 84 LWing +21.7202 -25.8924 1.54882 85 LWing +38.0238 0.0 1.0646 10 FusTail +34.6321 0.0 1.7946 37 VTP +35.0892 0.0 2.3205 38 VTP +35.5464 0.0 2.8464 39 VTP +36.0036 0.0 3.3723 40 VTP +36.4607 0.0 3.8982 41 VTP +36.9003 0.0 4.4039 42 VTP +37.3398 0.0 4.9095 43 VTP +37.7793 0.0 5.4151 44 VTP +38.3537 0.0 5.9207 46 HTP +38.7465 0.0 6.5277 45 VTPTail +38.3537 0.1261 5.9207 47 RHTP +38.3537 0.2522 5.9207 48 RHTP +38.5245 0.7876 5.8927 49 RHTP +38.6953 1.323 5.8646 50 RHTP +38.8661 1.8583 5.8366 51 RHTP +39.037 2.3937 5.8085 52 RHTP +39.2078 2.9291 5.7805 53 RHTP +39.3635 3.4172 5.7549 54 RHTP +39.5193 3.9053 5.7293 55 RHTP +39.675 4.3934 5.7037 56 RHTP +39.8308 4.8816 5.6781 57 RHTP +39.9865 5.3697 5.6525 58 RHTP +40.1423 5.8578 5.627 59 RHTP +38.3537 -0.1261 5.9207 86 LHTP +38.3537 -0.2522 5.9207 87 LHTP +38.5245 -0.7876 5.8927 88 LHTP +38.6953 -1.323 5.8646 89 LHTP +38.8661 -1.8583 5.8366 90 LHTP +39.037 -2.3937 5.8085 91 LHTP +39.2078 -2.9291 5.7805 92 LHTP +39.3635 -3.4172 5.7549 93 LHTP +39.5193 -3.9053 5.7293 94 LHTP +39.675 -4.3934 5.7037 95 LHTP +39.8308 -4.8816 5.6781 96 LHTP +39.9865 -5.3697 5.6525 97 LHTP +40.1423 -5.8578 5.627 98 LHTP diff --git a/examples/BUG/NASTRAN/Asets/asets_free_reduced.bdf b/examples/BUG/NASTRAN/Asets/asets_free_reduced.bdf index feb722d..f127296 100644 --- a/examples/BUG/NASTRAN/Asets/asets_free_reduced.bdf +++ b/examples/BUG/NASTRAN/Asets/asets_free_reduced.bdf @@ -4,12 +4,12 @@ $pyNastran: encoding=utf-8 $pyNastran: nnodes=0 $pyNastran: nelements=0 $SETS -ASET1 123456 1000 THRU 1004 1006 THRU 1010 2001 - 2003 2005 2008 2010 2012 2014 2016 2018 - 2020 2022 2024 2026 2028 2030 2032 2034 - 2036 2038 2040 2042 2044 2046 2048 2050 - 2052 3001 THRU 3008 3010 4000 4001 THRU - 401310002001100020031000200510002008100020101000201210002014 - 1000201610002018100020201000202210002024100020261000202810002030 - 1000203210002034100020361000203810002040100020421000204410002046 - 10002048100020501000205210004001 THRU10004013 1005 +ASET1 123456 1000 THRU 1010 2001 2003 2005 2008 + 2010 2012 2014 2016 2018 2020 2022 2024 + 2026 2028 2030 2032 2034 2036 2038 2040 + 2042 2044 2046 2048 2050 2052 3001 THRU + 3008 3010 4000 4001 THRU 40131000200110002003 + 1000200510002008100020101000201210002014100020161000201810002020 + 1000202210002024100020261000202810002030100020321000203410002036 + 1000203810002040100020421000204410002046100020481000205010002052 + 10004001 THRU10004013 diff --git a/examples/BUG/NASTRAN/BUG_103cam.bdf b/examples/BUG/NASTRAN/BUG_103cam.bdf deleted file mode 100644 index 7e6c192..0000000 --- a/examples/BUG/NASTRAN/BUG_103cam.bdf +++ /dev/null @@ -1,76 +0,0 @@ -$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ -$ EXECUTIVE CONTROL $ -$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ -$ TIME 100 $(Max execution time) -assign output4='./output_fem/Phia.op4',formatted,UNIT=11 -assign output4='./output_fem/Ma.op4',formatted,UNIT=12 -assign output4='./output_fem/Ka.op4',formatted,UNIT=13 -SOL 103 -DIAG 20 -COMPILE SEMODES -ALTER 'CALL.*SUPER3.*CASECC.*CASEM.*PHA1.*LAMA2' -OUTPUT4 PHA,,,,//0/11///20 -OUTPUT4 MAA,,,,//0/12///20 -OUTPUT4 KAA,,,,//0/13///20 -CEND - -$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ -$ CASE CONTROL $ -$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ -TITLE=BUG model # -ECHO=NONE -SPC = 1 -SPCF = ALL -DISPLACEMENT=ALL -METHOD = 100 - -$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ -$ BULK $ -$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ - -BEGIN BULK -PARAM,AUTOMSET,YES -$ PARAM,BAILOUT,-1 -$ PARAM,GRDPNT,0 -$ PARAM,K6ROT,1.0 -$ PARAM,SNORM,20.0 -PARAM,POST,-1 -$ PARAM,MAXRATIO,1.0E07 -$ PARAM,EXTOUT,DMIGPCH -EIGRL,100,,,100 - -$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ -$ MODEL $ -$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ - -INCLUDE './Parts/BUG_Fuselage_VTP.bdf' -INCLUDE './Parts/BUG_WING_LWBOX.bdf' -INCLUDE './Parts/BUG_WING_RWBOX.bdf' -INCLUDE './Parts/fuselage_LWBOX_rbe.bdf' -INCLUDE './Parts/fuselage_RWBOX_rbe.bdf' -INCLUDE './Parts/MTOW_FUEL_LWBOX.bdf' -INCLUDE './Parts/MTOW_FUEL_RWBOXmod.bdf' -INCLUDE './Parts/BUG_LHTP.bdf' -INCLUDE './Parts/BUG_RHTP.bdf' - -$$$$$$$$$$$$$$$$$$$$$$$$$ -$ SPCS $ -$$$$$$$$$$$$$$$$$$$$$$$$$ - -$ INCLUDE './Config/spcs.bdf' - -$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ -$ CLAMPING NODE $ -$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ - -SPC1 1 123456 1005 - - -$$$$$$$$$$$$$$$$$$$$$$$$$ -$ ASETs $ -$$$$$$$$$$$$$$$$$$$$$$$$$ - -INCLUDE './Config/asets_clamped_reduced.bdf' - -$$$$$$$$$$$$$$$$$$$$$$$$$$$ -ENDDATA diff --git a/examples/BUG/NASTRAN/BUG_103cfo.bdf b/examples/BUG/NASTRAN/BUG_103cfo.bdf deleted file mode 100644 index d2b9366..0000000 --- a/examples/BUG/NASTRAN/BUG_103cfo.bdf +++ /dev/null @@ -1,61 +0,0 @@ -$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ -$ EXECUTIVE CONTROL $ -$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ - -SOL 103 -CEND -$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ -$ CASE CONTROL $ -$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ - -TITLE=BUG model # -ECHO=NONE -SPC = 1 -SPCF = ALL -DISPLACEMENT=ALL -METHOD = 100 - -$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ -$ BULK $ -$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ - -BEGIN BULK -PARAM,AUTOMSET,YES -$ PARAM,BAILOUT,-1 -$ PARAM,GRDPNT,0 -$ PARAM,K6ROT,1.0 -$ PARAM,SNORM,20.0 -PARAM,POST,-1 -$ PARAM,MAXRATIO,1.0E07 -$ PARAM,EXTOUT,DMIGPCH -EIGRL,100,,,100 - -$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ -$ MODEL $ -$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ - -INCLUDE './Parts/BUG_Fuselage_VTP.bdf' -INCLUDE './Parts/BUG_WING_LWBOX.bdf' -INCLUDE './Parts/BUG_WING_RWBOX.bdf' -INCLUDE './Parts/fuselage_LWBOX_rbe.bdf' -INCLUDE './Parts/fuselage_RWBOX_rbe.bdf' -INCLUDE './Parts/MTOW_FUEL_LWBOX.bdf' -INCLUDE './Parts/MTOW_FUEL_RWBOXmod.bdf' -INCLUDE './Parts/BUG_LHTP.bdf' -INCLUDE './Parts/BUG_RHTP.bdf' - -$$$$$$$$$$$$$$$$$$$$$$$$$ -$ SPCS $ -$$$$$$$$$$$$$$$$$$$$$$$$$ - -$ INCLUDE './Config/spcs.bdf' - -$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ -$ CLAMPING NODE $ -$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ - -SPC1 1 123456 1005 - - -$$$$$$$$$$$$$$$$$$$$$$$$$$$ -ENDDATA diff --git a/examples/BUG/NASTRAN/BUG_103eao.bdf b/examples/BUG/NASTRAN/BUG_103eao.bdf deleted file mode 100644 index 765ca5b..0000000 --- a/examples/BUG/NASTRAN/BUG_103eao.bdf +++ /dev/null @@ -1,79 +0,0 @@ -$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ -$ EXECUTIVE CONTROL $ -$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ - -SOL 103 -CEND -$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ -$ CASE CONTROL $ -$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ -SET 1 = 1006,1007,1008,1009,1004,1003,1002,1001, -1000,2001,2003,2005,2008,2010,2012,2014, -2016,2018,2020,2022,2024,2026,2028,2030, -2032,2034,2036,2038,2040,2042,2044,2046, -2048,2050,2052,10002001,10002003,10002005,10002008,10002010, -10002012,10002014,10002016,10002018,10002020,10002022,10002024,10002026, -10002028,10002030,10002032,10002034,10002036,10002038,10002040,10002042, -10002044,10002046,10002048,10002050,10002052,1010,3001,3002, -3003,3004,3005,3006,3007,3008,4000,3010, -4001,4002,4003,4004,4005,4006,4007,4008, -4009,4010,4011,4012,4013,10004001,10004002,10004003, -10004004,10004005,10004006,10004007,10004008,10004009,10004010,10004011, -10004012,10004013,1005 -TITLE=BUG model # -ECHO=NONE -$ SPC = 1 -$ SPCF = ALL -DISPLACEMENT=1 -METHOD = 100 - -$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ -$ BULK $ -$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ - -BEGIN BULK -PARAM,AUTOMSET,YES -$ PARAM,BAILOUT,-1 -$ PARAM,GRDPNT,0 -$ PARAM,K6ROT,1.0 -$ PARAM,SNORM,20.0 -PARAM,POST,-1 -$ PARAM,MAXRATIO,1.0E07 -$ PARAM,EXTOUT,DMIGPCH -EIGRL,100,,,100 - -$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ -$ MODEL $ -$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ - -INCLUDE './Parts/BUG_Fuselage_VTP.bdf' -INCLUDE './Parts/BUG_WING_LWBOX.bdf' -INCLUDE './Parts/BUG_WING_RWBOX.bdf' -INCLUDE './Parts/fuselage_LWBOX_rbe.bdf' -INCLUDE './Parts/fuselage_RWBOX_rbe.bdf' -INCLUDE './Parts/MTOW_FUEL_LWBOX.bdf' -INCLUDE './Parts/MTOW_FUEL_RWBOXmod.bdf' -INCLUDE './Parts/BUG_LHTP.bdf' -INCLUDE './Parts/BUG_RHTP.bdf' - -$$$$$$$$$$$$$$$$$$$$$$$$$ -$ SPCS $ -$$$$$$$$$$$$$$$$$$$$$$$$$ - -$ INCLUDE './Config/spcs.bdf' - -$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ -$ CLAMPING NODE $ -$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ - -$ SPC1 1 123456 1005 - - -$$$$$$$$$$$$$$$$$$$$$$$$$ -$ ASETs $ -$$$$$$$$$$$$$$$$$$$$$$$$$ - -INCLUDE './Config/asets_free_reduced.bdf' - -$$$$$$$$$$$$$$$$$$$$$$$$$$$ -ENDDATA diff --git a/examples/BUG/NASTRAN/BUG_103efo.bdf b/examples/BUG/NASTRAN/BUG_103efo.bdf deleted file mode 100644 index 8a396c1..0000000 --- a/examples/BUG/NASTRAN/BUG_103efo.bdf +++ /dev/null @@ -1,55 +0,0 @@ -$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ -$ EXECUTIVE CONTROL $ -$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ - -SOL 103 -CEND -$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ -$ CASE CONTROL $ -$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ - -TITLE=BUG model # -ECHO=NONE -$ SPC = 1 -$ SPCF = ALL -DISPLACEMENT=ALL -METHOD = 100 - -$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ -$ BULK $ -$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ - -BEGIN BULK -PARAM,AUTOMSET,YES -$ PARAM,BAILOUT,-1 -$ PARAM,GRDPNT,0 -$ PARAM,K6ROT,1.0 -$ PARAM,SNORM,20.0 -PARAM,POST,-1 -$ PARAM,MAXRATIO,1.0E07 -$ PARAM,EXTOUT,DMIGPCH -EIGRL,100,,,100 - -$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ -$ MODEL $ -$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ - -INCLUDE './Parts/BUG_Fuselage_VTP.bdf' -INCLUDE './Parts/BUG_WING_LWBOX.bdf' -INCLUDE './Parts/BUG_WING_RWBOX.bdf' -INCLUDE './Parts/fuselage_LWBOX_rbe.bdf' -INCLUDE './Parts/fuselage_RWBOX_rbe.bdf' -INCLUDE './Parts/MTOW_FUEL_LWBOX.bdf' -INCLUDE './Parts/MTOW_FUEL_RWBOXmod.bdf' -INCLUDE './Parts/BUG_LHTP.bdf' -INCLUDE './Parts/BUG_RHTP.bdf' - -$$$$$$$$$$$$$$$$$$$$$$$$$ -$ SPCS $ -$$$$$$$$$$$$$$$$$$$$$$$$$ - -$ INCLUDE './Config/spcs.bdf' - - -$$$$$$$$$$$$$$$$$$$$$$$$$$$ -ENDDATA diff --git a/examples/BUG/NASTRAN/BUGaero1.bdf b/examples/BUG/NASTRAN/BUGaero1.bdf index 0219e80..da59dfc 100644 --- a/examples/BUG/NASTRAN/BUGaero1.bdf +++ b/examples/BUG/NASTRAN/BUGaero1.bdf @@ -7,24 +7,8 @@ CEND $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ $ CASE CONTROL $ $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ -SET 1 = 1006,1007,1008,1009,1004,1003,1002,1001, -1000,2001,2003,2005,2008,2010,2012,2014, -2016,2018,2020,2022,2024,2026,2028,2030, -2032,2034,2036,2038,2040,2042,2044,2046, -2048,2050,2052,10002001,10002003,10002005,10002008,10002010, -10002012,10002014,10002016,10002018,10002020,10002022,10002024,10002026, -10002028,10002030,10002032,10002034,10002036,10002038,10002040,10002042, -10002044,10002046,10002048,10002050,10002052,1010,3001,3002, -3003,3004,3005,3006,3007,3008,4000,3010, -4001,4002,4003,4004,4005,4006,4007,4008, -4009,4010,4011,4012,4013,10004001,10004002,10004003, -10004004,10004005,10004006,10004007,10004008,10004009,10004010,10004011, -10004012,10004013,1005 TITLE=BUG model # ECHO=NONE -SPC = 1 -SPCF = ALL -DISPLACEMENT=1 METHOD = 100 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ @@ -33,47 +17,15 @@ $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ BEGIN BULK PARAM,AUTOMSET,YES -$ PARAM,BAILOUT,-1 -$ PARAM,GRDPNT,0 -$ PARAM,K6ROT,1.0 -$ PARAM,SNORM,20.0 PARAM,POST,-1 -$ PARAM,MAXRATIO,1.0E07 -$ PARAM,EXTOUT,DMIGPCH EIGRL,100,,,100 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ $ MODEL $ $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ -INCLUDE './Parts/BUG_Fuselage_VTP.bdf' -INCLUDE './Parts/BUG_WING_LWBOX.bdf' -INCLUDE './Parts/BUG_WING_RWBOX.bdf' -INCLUDE './Parts/fuselage_LWBOX_rbe.bdf' -INCLUDE './Parts/fuselage_RWBOX_rbe.bdf' -INCLUDE './Parts/MTOW_FUEL_LWBOX.bdf' -INCLUDE './Parts/MTOW_FUEL_RWBOXmod.bdf' -INCLUDE './Parts/BUG_LHTP.bdf' -INCLUDE './Parts/BUG_RHTP.bdf' +INCLUDE './base_model.bdf' INCLUDE './OldAero/aero_rightonly.bdf' -$$$$$$$$$$$$$$$$$$$$$$$$$ -$ SPCS $ -$$$$$$$$$$$$$$$$$$$$$$$$$ - -$ INCLUDE './Config/spcs.bdf' - -$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ -$ CLAMPING NODE $ -$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ - -SPC1 1 123456 1005 - - -$$$$$$$$$$$$$$$$$$$$$$$$$ -$ ASETs $ -$$$$$$$$$$$$$$$$$$$$$$$$$ - -INCLUDE './Config/asets_clamped_reduced.bdf' $$$$$$$$$$$$$$$$$$$$$$$$$$$ ENDDATA diff --git a/examples/BUG/Streamlit/pages/7_Shards.py b/examples/BUG/Streamlit/pages/7_Shards.py index 6856387..aa5d616 100644 --- a/examples/BUG/Streamlit/pages/7_Shards.py +++ b/examples/BUG/Streamlit/pages/7_Shards.py @@ -4,7 +4,7 @@ import streamlit as st import importlib -# importlib.reload(sti) +importlib.reload(sti) st.set_page_config( page_title="Shards systems", @@ -41,7 +41,7 @@ config = configuration.Config.from_file(f"{si}/config.yaml") Csol[si.name] = sol Cconfig[si.name] = config - st.session_state['Csol'] = Csol + st.session_state['Csol'] = Csol st.session_state['Cconfig'] = Cconfig sti.systems_shard(st.session_state['Csolshard'], diff --git a/examples/BUG/modelgeneration.org b/examples/BUG/modelgeneration.org index ba35996..b0741b5 100644 --- a/examples/BUG/modelgeneration.org +++ b/examples/BUG/modelgeneration.org @@ -5,14 +5,14 @@ (pyvenv-workon "feniax") #+end_src -* Folder structure +* Folder and file structure #+begin_src markdown :tangle "./README.md" :results none # Folder structure - Parts: BDF files with the structural model parts - DLMs: generated DLM models #+end_src -#+begin_src python +#+begin_src python :results none import pathlib pathlib.Path('./FEM').mkdir(parents=True, exist_ok=True) pathlib.Path('./AERO').mkdir(parents=True, exist_ok=True) @@ -23,14 +23,24 @@ pathlib.Path('./NASTRAN/simulations_out').mkdir(parents=True, exist_ok=True) pathlib.Path('./NASTRAN/Gusts').mkdir(parents=True, exist_ok=True) pathlib.Path('./figs').mkdir(parents=True, exist_ok=True) - #+end_src +#+end_src -#+RESULTS: -: None +#+begin_src shell :tangle ./workflow.sh + python P1_modalsolution.py + source P2_runmodal.sh # Run Nastran modal solution + python P3_dlm.py + source P4_rundihedral.sh + python P5_getdihedral.py + python P6_outputgafs.py + source P7_rungafs.sh + source P8_rungafSteady.sh + python P9_getgafSteady.py + python P10_rogerRFA.py +#+end_src * Modal Solution :PROPERTIES: -:header-args: :tangle ./modal_solution.py :session *pybug* :comments yes +:header-args: :tangle ./P1_modalsolution.py :session *pybug* :comments yes :END: ** INPUT PARAMETERS @@ -42,7 +52,7 @@ import feniax.plotools.nastranvtk.bdfdef as bdfdef import numpy as np - num_modes = 50 + num_modes = 100 sol = "eao" # {c,e}{a,f}{o,p} WRITE_GRID= True WRITE_ASETS= True @@ -278,10 +288,9 @@ ENDDATA #+end_src ** Run Nastran -Running Nastran using the tailored functions in run_nastra.sh which moves output files and checks for fatal errors. +Running Nastran using the tailored functions in run_nastra.sh which moves output files and checks for fatal errors. -#+begin_src bash :session shell1 :noweb yes :tangle run_modal.sh :results none - cd NASTRAN +#+begin_src bash :session shell1 :noweb yes :tangle P2_runmodal.sh :results none :dir ./NASTRAN source run_nastran.sh run_nastran BUG103_<>.bdf move_outputs BUG103_<>.bdf @@ -335,6 +344,7 @@ Running Nastran using the tailored functions in run_nastra.sh which moves output #+end_src *** Plot VTK modes + Plot the modal shapes in Paraview #+begin_src python op2_file = f"./NASTRAN/simulations_out/BUG103_{sol}.op2" @@ -360,18 +370,21 @@ Read the pch file and save FE matrices to FEM folder if SAVE_FE: np.save(f"./FEM/Ka_{soli}.npy", stiffnessMatrix) np.save(f"./FEM/Ma_{soli}.npy", massMatrix) - assert len(asets_indexes) == len(stiffnessMatrix), "the FE matrices size does not match the indexes used to build the aset modes from the full set" + try: + assert len(asets_indexes) == len(stiffnessMatrix), "the FE matrices size does not match the indexes used to build the aset modes from the full set" + except NameError: + print("Careful, no aset-matrix sizes checked") #+end_src #+RESULTS: * DLM generation :PROPERTIES: - :header-args: :session *pybug* :tangle ./dlm.py :comments yes + :header-args: :session *pybug* :comments yes :END: ** INPUT PARAMETERS #+NAME: parameters_dlm0 -#+begin_src python +#+begin_src python :tangle ./P3_dlm.py import json import feniax.unastran.aero as aero from pyNastran.bdf.bdf import BDF @@ -424,7 +437,7 @@ Read the pch file and save FE matrices to FEM folder ** Build Build the aero model based on discretisation and the right-hand side aero built initially in BUGaero1.bdf #+NAME: DLMbuild -#+begin_src python :results none +#+begin_src python :results none :tangle ./P3_dlm.py # Read old model with right side of CAEROS bdfaero = BDF()#debug=False) @@ -471,7 +484,7 @@ Build the aero model based on discretisation and the right-hand side aero built #+end_src -** Paraview postprocessing +** DLM grid: Collocation points and Paraview plotting Old method: build panel coordinates out of corner points #+NAME: DLMparaview @@ -480,9 +493,9 @@ Old method: build panel coordinates out of corner points gridmesh = panels.build_gridmesh(grid, label_dlm, save_dir=f"./paraview/aero{label_dlm}") # write paraview mesh #+end_src -Use pynastran via get_collocation. Note: should push fix to pyNastran -#+NAME: DLM -#+begin_src python :results none +Use pynastran via get_collocation. Note: should push fix to pyNastran. +#+NAME: DLMGrid +#+begin_src python :results none :tangle ./P3_dlm.py dlmgrid = aero.GenDLMGrid(dlm.model) dlmgrid.plot_pyvista(f"./paraview/dlm{label_dlm}") @@ -505,7 +518,6 @@ Alternative: export to cquads and use old codes #bdfdef.vtkRef("./NASTRAN/Paraview/BUG_103cao.bdf") # write full FE paraview #+end_src - ** Dihedral extraction Basically extracting the value of for the normal of each panel that needs to be multiplied by *** bdf models @@ -1343,7 +1355,7 @@ Basically extracting the value of for the normal of each panel that needs to be *** Run Nastran Running Nastran using the tailored functions in run_nastra.sh which moves output files and checks for fatal errors. -#+begin_src bash :session shell1 :results none :tangle run_gafs.sh +#+begin_src bash :session shell1 :results none :tangle P4_rundihedral.sh cd NASTRAN source run_nastran.sh run_nastran BUGdihedral.bdf @@ -1352,8 +1364,9 @@ Running Nastran using the tailored functions in run_nastra.sh which moves output *** Read op4 -#+begin_src python :results none :noweb yes - +#+begin_src python :results none :noweb yes :tangle P5_getdihedral.py + import numpy as np + import feniax.unastran.op4handler as op4handler dihedral = op4handler.read_data(f'./NASTRAN/data_out/Dihedral.op4', 'WJ') SAVE_DIHEDRAL = True @@ -1364,7 +1377,7 @@ Running Nastran using the tailored functions in run_nastra.sh which moves output * GAFs extraction :PROPERTIES: -:header-args: :session *pybug* :tangle ./gafs.py :comments yes :noweb yes +:header-args: :session *pybug* :tangle ./P6_outputgafs.py :comments yes :noweb yes :END: ** INPUT PARAMETERS #+NAME: parameters_gafs0 @@ -1372,11 +1385,13 @@ Running Nastran using the tailored functions in run_nastra.sh which moves output import numpy as np import feniax.unastran.aero as nasaero import feniax.unastran.op4handler as op4handler + from feniax.utils import standard_atmosphere import pickle import itertools sol = "eao" - num_modes = 50 - mach = 0.8 + num_modes = 100 + mach = 0.7 + altitude = 10000 # meters Mach = str(mach).replace('.','_') machs = [mach] reduced_freqs = np.hstack([1e-6, np.linspace(1e-5,1e-1, 25), @@ -1385,12 +1400,12 @@ Running Nastran using the tailored functions in run_nastra.sh which moves output reduced_freqs = np.hstack([1e-5, np.linspace(1e-4, 1, 100) ]) #reduced_freqs = np.geomspace(1e-5, 1, 100, endpoint=True) - num_modes = 50 flutter_id = 9010 mach_fact = machs kv_fact = [200., 220.] - u_inf = 200. - rho_inf = 1.2 + T, rho_inf, P, a = standard_atmosphere(altitude) + u_inf = mach * a + #rho_inf = 1.2 density_fact = [rho_inf] chord_ref = 3. span_ref = 24. * 2 # always full span @@ -1401,7 +1416,7 @@ Running Nastran using the tailored functions in run_nastra.sh which moves output flutter_sett = dict() aero_sett = dict() label_dlm = "<>" - label_flow = f"F1" + label_flow = f"F3" label_gaf = f"D{label_dlm}{label_flow}S{sol}-{num_modes}" input_dict = dict(reduced_freqs=list(reduced_freqs), mach=mach, u_inf=u_inf, rho_inf=rho_inf) with open(f"./NASTRAN/GAFs/input_{label_flow}.json", "w") as fp: @@ -1525,24 +1540,169 @@ TODO: add ASETs and check whether it affects the results *** Run Nastran Running Nastran using the tailored functions in run_nastra.sh which moves output files and checks for fatal errors. -#+begin_src bash :session shell1 :results none :tangle run_gafs.sh - cd NASTRAN +#+begin_src bash :session shell1 :results none :tangle P7_rungafs.sh :dir NASTRAN source run_nastran.sh run_nastran BUGgafs_<>.bdf move_outputs BUGgafs_<>.bdf #+end_src +** Steady +*** bdf models +:PROPERTIES: +:header-args: :noweb yes :comments no +:END: + +#+NAME: GAFsSteady_setup +#+begin_src org :tangle no + $ + $--------------------------------------------------------------------------- + $ AERODYNAMIC DOFS + $--------------------------------------------------------------------------- + $ + AESTAT 1 ANGLEA + AESTAT 2 SIDES + AESTAT 3 PITCH + AESTAT 4 ROLL + AESTAT 5 YAW + AESTAT 6 URDD1 + AESTAT 7 URDD2 + AESTAT 8 URDD3 + AESTAT 9 URDD4 + AESTAT 10 URDD5 + AESTAT 11 URDD6 + + $AEROS 4.163 44.8 146.6 + AEROS,,,<>,<>,<> + $ + $ + $ TRIM 960 0.81 15762.81 URDD1 0. URDD2 0. 1. + $ URDD3 0. URDD4 0. URDD5 0. URDD6 0. + $ ROLL 0. YAW 0. SIDES 0. PITCH 0. + $ Flprn_r 0. WTAil_r 0. Elev_r 0. ANGLEA 0.261799 + TRIM,960,<>,<>,URDD1,0.,URDD2,0.,1., + ,URDD3,0.,URDD4,0.,URDD5,0.,URDD6,0., + ,ROLL,0.,YAW,0.,SIDES,0.,PITCH,0. + $ANGLEA 0.261799 + +#+end_src + +#+NAME: GAFsSteady_case +#+begin_src org :tangle no + NASTRAN QUARTICDLM=1 + SOL 144 + INCLUDE './DMAPs/Qhx.bdf' + CEND + + $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + $ CASE CONTROL $ + $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + + TITLE=BUG model # + TRIM = 960 + $LOAD = 2000 + $ + DISP = ALL + FORCE = ALL + AEROF = ALL + MONITOR = ALL + TRIMF = ALL + OLOAD(CID) = ALL + ECHO=NONE +#+end_src + +#+NAME: GAFsSteady_bulk +#+begin_src org :tangle no + BEGIN BULK + PARAM,BAILOUT,0 + PARAM,GRDPNT,0 + PARAM,K6ROT,1.0 + PARAM,SNORM,20.0 + PARAM,POST,0 + $PARAM,MAXRATIO,1.0E07 $Default anyway + $PARAM AUTOSPC YES + MDLPRM MLTSPLIN 1 $Aero grids can be defined in multiple splines (dafault 0) + PARAM WTMASS 1.0 + + INCLUDE './base_model.bdf' + INCLUDE './DLMs/<>.bdf' +#+end_src + +**** cao +#+begin_src org :tangle "./NASTRAN/BUGgafsSteady_cao.bdf" + assign OUTPUT4='./data_out/Qax<>.op4',formatted,UNIT=11 + assign OUTPUT4='./data_out/Qah<>.op4',formatted,UNIT=12 + assign OUTPUT4='./data_out/Qhx<>.op4',formatted,UNIT=13 + assign INPUTT4='./data_out/Phi<>_<>.op4',formatted,UNIT=90 + + <> + SPC = 1 + <> + <> + + $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + $ CLAMPING NODE $ + $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + + SPC1 1 12346 1005 + SUPORT 1005 5 +#+end_src + +*** Run Nastran +Running Nastran using the tailored functions in run_nastra.sh which moves output files and checks for fatal errors. + +#+begin_src bash :session shell1 :results none :noweb yes :tangle P8_rungafSteady.sh + cd NASTRAN + source run_nastran.sh + run_nastran BUGgafsSteady_<>.bdf + move_outputs BUGgafsSteady_<>.bdf +#+end_src + +*** Read op4 + +#+begin_src python :results none :noweb yes :tangle P9_getgafSteady.py + import numpy as np + import feniax.unastran.op4handler as op4handler + + Qax_name = "Qax<>" + Qah_name = "Qah<>" + Qhx_name = "Qhx<>" + Qax = op4handler.read_data(f'./NASTRAN/data_out/{Qax_name}.op4', + 'Q_AX') + Qah = op4handler.read_data(f'./NASTRAN/data_out/{Qah_name}.op4', + 'Q_AH') + Qhx = op4handler.read_data(f'./NASTRAN/data_out/{Qhx_name}.op4', + 'Q_HX') + SAVE_Qx = True + if SAVE_Qx: + np.save(f"./AERO/{Qax_name}.npy", Qax) + np.save(f"./AERO/{Qah_name}.npy", Qah) + np.save(f"./AERO/{Qhx_name}.npy", Qhx) + +#+end_src -*** Roger RFA +* Roger RFA :PROPERTIES: :header-args: :session *pybug* :noweb yes :END: -#+begin_src python :results none :tangle rogerRFA.py +** INPUT PARAMETERS +#+NAME: parameters_rfa0 +#+begin_src python + num_poles = 5 + Dhj_file = f"D{label_gaf}p{num_poles}" + Ahh_file = f"A{label_gaf}p{num_poles}" + Poles_file = f"Poles{label_gaf}p{num_poles}" +#+end_src + +#+RESULTS: parameters_rfa0 + +#+begin_src python :results none :tangle P10_rogerRFA.py import numpy as np import feniax.unastran.aero as nasaero import feniax.unastran.op4handler as op4handler import pickle - + import importlib + import feniax.aeromodal.roger as roger + importlib.reload(roger) #op4m = op4.OP4() #Qop4 = op4m.read_op4(file_name) @@ -1552,17 +1712,11 @@ Running Nastran using the tailored functions in run_nastra.sh which moves output 'Q_HJ') #+end_src -**** RFA for Qhh and Qhj +** RFA for Qhh and Qhj - Pole optimisation - #+begin_src python :results none :noweb yes :tangle rogerRFA.py - import importlib - import feniax.aeromodal.roger as roger - importlib.reload(roger) - num_poles = 5 - Dhj_file = f"D{label_gaf}p{num_poles}" - Ahh_file = f"A{label_gaf}p{num_poles}" - Poles_file = f"Poles{label_gaf}p{num_poles}" + #+begin_src python :results none :tangle P10_rogerRFA.py + <> optpoles = roger.OptimisePoles(reduced_freqs, Qhh, num_poles_=num_poles, poles_step_=0.3, @@ -1590,7 +1744,7 @@ Running Nastran using the tailored functions in run_nastra.sh which moves output Qroger3 = qhhr3.eval_array(reduced_freqs) #+end_src - #+begin_src python :results none :noweb yes :tangle rogerRFA.py + #+begin_src python :results none :tangle P10_rogerRFA.py poles = qhhr3.poles #jnp.load("./AERO/PolesDd1c7F1Scao-50.npy") rogerhj = roger.ComputeRoger(Qhj, reduced_freqs, poles, 2) np.save(f"./AERO/{Dhj_file}.npy", rogerhj.roger_matrices) @@ -1599,7 +1753,7 @@ Running Nastran using the tailored functions in run_nastra.sh which moves output #+end_src - Double number of poles - #+begin_src python :results none :noweb yes :tangle rogerRFA.py + #+begin_src python :results none num_poles *= 2 Dhj_file = f"D{label_gaf}p{num_poles}" Ahh_file = f"A{label_gaf}p{num_poles}" @@ -1619,7 +1773,7 @@ Running Nastran using the tailored functions in run_nastra.sh which moves output Qroger4 = qhhr4.eval_array(reduced_freqs) #+end_src - #+begin_src python :results none :noweb yes :tangle rogerRFA.py + #+begin_src python :results none poles = qhhr4.poles #jnp.load("./AERO/PolesDd1c7F1Scao-50.npy") rogerhj = roger.ComputeRoger(Qhj, reduced_freqs, poles, 2) np.save(f"./AERO/{Dhj_file}.npy", rogerhj.roger_matrices) @@ -1628,7 +1782,7 @@ Running Nastran using the tailored functions in run_nastra.sh which moves output #+end_src - n times more poles, no optimisation - #+begin_src python :results none :noweb yes :tangle rogerRFA.py + #+begin_src python :results none n = 1.4 num_poles = 7 # *= n num_poles = int(num_poles) @@ -1648,22 +1802,22 @@ Running Nastran using the tailored functions in run_nastra.sh which moves output Qrogerhj = rogerhjeval.eval_array(reduced_freqs) #+end_src - #+begin_src python :results none :noweb yes :tangle rogerRFA.py + #+begin_src python :results none roger.plot_gafs(0, 0, Qhh, [Qrogerhh]) #roger.plot_gafs(20, 2, Qhh, [Qroger, Qroger2, Qroger3]) #+end_src -**** Plot Qhh study -#+begin_src python :results none :noweb yes :tangle rogerRFA.py +** Plot Qhh study +#+begin_src python :results none roger.plot_gafs(0, 1, Qhh, [Qroger1, Qroger2, Qroger3, Qroger4]) #roger.plot_gafs(20, 2, Qhh, [Qroger, Qroger2, Qroger3]) #+end_src -#+begin_src python :results none :noweb yes :tangle rogerRFA.py +#+begin_src python :results none iterate_vect = list(range(10)) plot_prod = list(itertools.product(iterate_vect,iterate_vect)) for li in plot_prod: @@ -1671,7 +1825,7 @@ Running Nastran using the tailored functions in run_nastra.sh which moves output roger.plot_gafs(li[0], li[1], Qhh, [Qroger1]) #+end_src -#+begin_src python :results none :noweb yes :tangle rogerRFA.py +#+begin_src python :results none iterate_vect = list(range(10)) plot_prod = list(itertools.product(iterate_vect,iterate_vect)) for li in plot_prod: @@ -1679,32 +1833,15 @@ Running Nastran using the tailored functions in run_nastra.sh which moves output roger.plot_gafs(li[0], li[1], Qhh, [Qroger1, Qroger2, Qroger3, Qroger4]) #+end_src - -**** Plot Qhh -#+begin_src python :results none :noweb yes :tangle rogerRFA.py - - roger.plot_gafs(0, 1, Qhh, [Qroger1, Qroger2, Qroger3, Qroger4]) - - #roger.plot_gafs(20, 2, Qhh, [Qroger, Qroger2, Qroger3]) -#+end_src - -#+begin_src python :results none :noweb yes :tangle rogerRFA.py - iterate_vect = list(range(10)) - plot_prod = list(itertools.product(iterate_vect,iterate_vect)) - for li in plot_prod: - if np.linalg.norm(Qhh[:, li[0], li[1]]) > 1e-4: - roger.plot_gafs(li[0], li[1], Qhh, [Qrogerhh]) -#+end_src - -**** Plot Qhj -#+begin_src python :results none :noweb yes :tangle rogerRFA.py +** Plot Qhj +#+begin_src python :results none roger.plot_gafs(0, 1, Qhj, [Qrogerhj]) #roger.plot_gafs(20, 2, Qhh, [Qroger, Qroger2, Qroger3]) #+end_src -#+begin_src python :results none :noweb yes :tangle rogerRFA.py +#+begin_src python :results none iterate_vect = list(range(10)) plot_prod = list(range(20)) panel = 20 @@ -1713,159 +1850,46 @@ Running Nastran using the tailored functions in run_nastra.sh which moves output roger.plot_gafs(li, panel, Qhj, [Qrogerhj]) #+end_src -** Steady -*** bdf models +* FENIAX :PROPERTIES: -:header-args: :noweb yes :comments no +:header-args: :session *pybug* :comments yes :noweb yes :results none :END: -#+NAME: GAFsSteady_setup -#+begin_src org :tangle no - $ - $--------------------------------------------------------------------------- - $ AERODYNAMIC DOFS - $--------------------------------------------------------------------------- - $ - AESTAT 1 ANGLEA - AESTAT 2 SIDES - AESTAT 3 PITCH - AESTAT 4 ROLL - AESTAT 5 YAW - AESTAT 6 URDD1 - AESTAT 7 URDD2 - AESTAT 8 URDD3 - AESTAT 9 URDD4 - AESTAT 10 URDD5 - AESTAT 11 URDD6 - - $AEROS 4.163 44.8 146.6 - AEROS,,,<>,<>,<> - $ - $ - $ TRIM 960 0.81 15762.81 URDD1 0. URDD2 0. 1. - $ URDD3 0. URDD4 0. URDD5 0. URDD6 0. - $ ROLL 0. YAW 0. SIDES 0. PITCH 0. - $ Flprn_r 0. WTAil_r 0. Elev_r 0. ANGLEA 0.261799 - TRIM,960,<>,<>,URDD1,0.,URDD2,0.,1., - ,URDD3,0.,URDD4,0.,URDD5,0.,URDD6,0., - ,ROLL,0.,YAW,0.,SIDES,0.,PITCH,0. - $ANGLEA 0.261799 - +#+NAME: ImportsFeniax +#+begin_src python + import pathlib + import time + import jax.numpy as jnp + import numpy as np + import feniax.preprocessor.configuration as configuration # import Config, dump_to_yaml + from feniax.preprocessor.inputs import Inputs + import feniax.feniax_main + import feniax.plotools.reconstruction as reconstruction #+end_src -#+NAME: GAFsSteady_case -#+begin_src org :tangle no - NASTRAN QUARTICDLM=1 - SOL 144 - INCLUDE './DMAPs/Qhx.bdf' - CEND - - $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ - $ CASE CONTROL $ - $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ +#+NAME: ImportsFeniaxShard +#+begin_src python + import pathlib + import time + #import jax.numpy as jnp + import numpy as np + from feniax.preprocessor.inputs import Inputs + import feniax.feniax_shardmain +#+end_src - TITLE=BUG model # - TRIM = 960 - $LOAD = 2000 - $ - DISP = ALL - FORCE = ALL - AEROF = ALL - MONITOR = ALL - TRIMF = ALL - OLOAD(CID) = ALL - ECHO=NONE +#+NAME: MainBugParameters +#+begin_src python + label_dlm = "<>" + sol = "<>" + label_gaf = "<>" + num_modes = <> + c_ref = <> + u_inf = <> + rho_inf = <> #+end_src -#+NAME: GAFsSteady_bulk -#+begin_src org :tangle no - BEGIN BULK - PARAM,BAILOUT,0 - PARAM,GRDPNT,0 - PARAM,K6ROT,1.0 - PARAM,SNORM,20.0 - PARAM,POST,0 - $PARAM,MAXRATIO,1.0E07 $Default anyway - $PARAM AUTOSPC YES - MDLPRM MLTSPLIN 1 $Aero grids can be defined in multiple splines (dafault 0) - PARAM WTMASS 1.0 - - INCLUDE './base_model.bdf' - INCLUDE './DLMs/<>.bdf' -#+end_src - -**** cao -#+begin_src org :tangle "./NASTRAN/BUGgafsSteady_cao.bdf" - assign OUTPUT4='./data_out/Qax<>.op4',formatted,UNIT=11 - assign OUTPUT4='./data_out/Qah<>.op4',formatted,UNIT=12 - assign OUTPUT4='./data_out/Qhx<>.op4',formatted,UNIT=13 - assign INPUTT4='./data_out/Phi<>_<>.op4',formatted,UNIT=90 - - <> - SPC = 1 - <> - <> - - $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ - $ CLAMPING NODE $ - $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ - - SPC1 1 12346 1005 - SUPORT 1005 5 -#+end_src - -*** Run Nastran -Running Nastran using the tailored functions in run_nastra.sh which moves output files and checks for fatal errors. - -#+begin_src bash :session shell1 :results none :noweb yes :tangle run_gafs.sh - cd NASTRAN - source run_nastran.sh - run_nastran BUGgafsSteady_<>.bdf - move_outputs BUGgafsSteady_<>.bdf -#+end_src - -*** Read op4 - -#+begin_src python :results none :noweb yes :tangle rogerRFA.py - - Qax_name = "Qax<>" - Qah_name = "Qah<>" - Qhx_name = "Qhx<>" - Qax = op4handler.read_data(f'./NASTRAN/data_out/{Qax_name}.op4', - 'Q_AX') - Qah = op4handler.read_data(f'./NASTRAN/data_out/{Qah_name}.op4', - 'Q_AH') - Qhx = op4handler.read_data(f'./NASTRAN/data_out/{Qhx_name}.op4', - 'Q_HX') - SAVE_Qx = True - if SAVE_Qx: - np.save(f"./AERO/{Qax_name}.npy", Qax) - np.save(f"./AERO/{Qah_name}.npy", Qah) - np.save(f"./AERO/{Qhx_name}.npy", Qhx) - -#+end_src - -* FENIAX -:PROPERTIES: -:header-args: :session *pybug* :comments yes :noweb yes -:END: - -** Manoeuvre -*** Run -#+begin_src python :tangle ./settings_manoeuvre.py - import pathlib - import jax.numpy as jnp - import feniax.preprocessor.configuration as configuration # import Config, dump_to_yaml - from feniax.preprocessor.inputs import Inputs - import feniax.feniax_main - import feniax.plotools.reconstruction as reconstruction - - label_gaf = "<>" - num_poles = 5 - Dhj_file = f"D{label_gaf}p{num_poles}" - Ahh_file = f"A{label_gaf}p{num_poles}" - Poles_file = f"Poles{label_gaf}p{num_poles}" - +#+NAME: MainBug +#+begin_src python inp = Inputs() inp.engine = "intrinsicmodal" inp.fem.eig_type = "inputs" @@ -1884,141 +1908,183 @@ Running Nastran using the tailored functions in run_nastra.sh which moves output RHTP=None, LHTP=None, ) - - inp.fem.folder = pathlib.Path('./FEM/') - inp.fem.eig_names = ["eigenvals_50.npy", "eigenvecs_50.npy"] - inp.fem.num_modes = 50 + inp.fem.grid = f"./FEM/structuralGrid_{sol[:-1]}" + #inp.fem.folder = pathlib.Path('./FEM/') + inp.fem.Ka_name = f"./FEM/Ka_{sol[:-1]}.npy" + inp.fem.Ma_name = f"./FEM/Ma_{sol[:-1]}.npy" + inp.fem.eig_names = [f"./FEM/eigenvals_{sol}{num_modes}.npy", + f"./FEM/eigenvecs_{sol}{num_modes}.npy"] inp.driver.typeof = "intrinsic" - - #inp.driver.sol_path = pathlib.Path( - # f"./results_{datetime.datetime.now().strftime('%Y-%m-%d_%H:%M:%S')}") + inp.fem.num_modes = num_modes + inp.driver.typeof = "intrinsic" +#+end_src +** Discrete loads +Forces and moments at node 35 and 61 + +#+NAME: DiscreteLoads +#+begin_src python :tangle settings_DiscreteLoads.py + <> + sol = "cao" + num_modes = 300 + <> inp.driver.sol_path = pathlib.Path( - "./results1manoeuvre") + "./results/DiscreteLoads1") + inp.simulation.typeof = "single" - inp.systems.sett.s1.solution = "static" - inp.systems.sett.s1.solver_library = "diffrax" - inp.systems.sett.s1.solver_function = "newton" - inp.systems.sett.s1.solver_settings = dict(rtol=1e-6, + inp.system.name = "s1" + inp.system.solution = "static" + inp.system.solver_library = "diffrax" + inp.system.solver_function = "newton" + inp.system.solver_settings = dict(rtol=1e-6, + atol=1e-6, + max_steps=100, + norm="linalg_norm", + kappa=0.01) + inp.system.xloads.follower_forces = True + inp.system.xloads.x = [0, 1, 2, 3, 4, 5] + inp.system.t = [0.5, 1, 1.5, 2, 2.5, 3, 3.25, 3.5, 3.75, 4, 4.25, 4.5, 4.75, 5.] + lz1 = 5e4 * 0.5 + lz2 = 9e4 * 0.5 + lz3 = 2e5 * 0.5 + lz4 = 4e5 * 0.5 + lz5 = 5e5 * 0.5 + lx1 = lz1 * 5 + lx2 = lz2 * 5 + lx3 = lz3 * 5 + lx4 = lz4 * 5 + lx5 = lz5 * 5 + ly1 = lz1 * 7 + ly2 = lz2 * 7 + ly3 = lz3 * 7 + ly4 = lz4 * 7 + ly5 = lz5 * 7 + + # rwing: 14-35 + # lwing: 40-61 + inp.system.xloads.follower_points = [[35, 2], [61, 2], [35, 4], [61, 4]] + inp.system.xloads.follower_interpolation = [[0., lz1, lz2, lz3, lz4, lz5], + [0., lz1, lz2, lz3, lz4, lz5], + [0., lx1, lx2, lx3, lx4, lx5], + [0., lx1, lx2, lx3, lx4, lx5]] + t1 = time.time() + sol = feniax.feniax_main.main(input_dict=inp) + t2 = time.time() + print(f"Time DiscreteLoads: {t2 - t1}") +#+end_src + +** Manoeuvre + +!!Warning: label_gaf may have a solution with a free model (eao for instance) instead of a clamped, in which case it is not correct. + +#+NAME: ManoeuvreMain +#+begin_src python + <> + <> + <> + inp.simulation.typeof = "single" + inp.system.name = "s1" + inp.system.solution = "static" + inp.system.solver_library = "diffrax" + inp.system.solver_function = "newton" + inp.system.solver_settings = dict(rtol=1e-6, atol=1e-6, max_steps=100, norm="linalg_norm", kappa=0.01) - inp.systems.sett.s1.xloads.modalaero_forces = True - inp.systems.sett.s1.xloads.x = [0.,1.] - inp.systems.sett.s1.t = [0.25, 0.5, 0.75, 1] - inp.systems.sett.s1.aero.c_ref = <> - inp.systems.sett.s1.aero.u_inf = 234 # a0(7000) =312 - inp.systems.sett.s1.aero.rho_inf = 0.6 - inp.systems.sett.s1.aero.poles = f"./AERO/{Poles_file}.npy" - inp.systems.sett.s1.aero.A = f"./AERO/{Ahh_file}.npy" - inp.systems.sett.s1.aero.Q0_rigid = f"./AERO/Qhx{label_gaf}.npy" - inp.systems.sett.s1.aero.qalpha = jnp.array(([0., 0., 0, 0, 0, 0], - [0., 4 * jnp.pi / 180, 0, 0, 0, 0])) + inp.system.xloads.modalaero_forces = True + inp.system.xloads.x = [0.,1.] + inp.system.t = [1/6, 1/3, 1/2, 2/3, 5/6, 1]#[0.25, 0.5, 0.75, 1] + inp.system.aero.c_ref = c_ref + inp.system.aero.u_inf = u_inf # a0(7000) =312 + inp.system.aero.rho_inf = rho_inf + inp.system.aero.poles = f"./AERO/{Poles_file}.npy" + inp.system.aero.A = f"./AERO/{Ahh_file}.npy" + inp.system.aero.Q0_rigid = f"./AERO/Qhx{label_gaf}.npy" + inp.system.aero.qalpha = [[0., 0., 0, 0, 0, 0], + [0., 6 * np.pi / 180, 0, 0, 0, 0]] # interpolation: x=0 => qalpha=0 + # x=1 => qalpha = 4 +#+end_src +*** Run +#+begin_src python :tangle ./settings_manoeuvre1.py + <> + <> + inp.driver.sol_path = pathlib.Path( + "./results/manoeuvre2") config = configuration.Config(inp) + t1 = time.time() solstatic1 = feniax.feniax_main.main(input_obj=config) + t2 = time.time() + print(f"Time Manoeuvre: {t2 - t1}") + #+end_src -#+RESULTS: - *** Plot #+NAME: 3Dstatic -#+begin_src python :tangle ./settings_manoeuvre.py +#+begin_src python :tangle ./settings_manoeuvre1.py rintrinsic, uintrinsic = reconstruction.rbf_based( nastran_bdf="./NASTRAN/BUG103.bdf", X=config.fem.X, - time=range(len(inp.systems.sett.s1.t)), - ra=solstatic1.staticsystem_s1.ra, - Rab=solstatic1.staticsystem_s1.Cab, + time=range(len(inp.system.t)), + ra=solstatic1.staticsystem_sys1.ra, + Rab=solstatic1.staticsystem_sys1.Cab, R0ab=solstatic1.modes.C0ab, - vtkpath="./paraview/solstatic1/bug", + vtkpath=inp.driver.sol_path / "paraview/solstatic1/bug", plot_timeinterval=1, - plot_ref=True, + plot_ref=False, tolerance=1e-3, size_cards=8, rbe3s_full=False, ra_movie=None) #+end_src - ** Gust -*** Run -#+begin_src python :tangle settings_gust1.py - import pathlib - import jax.numpy as jnp - import feniax.preprocessor.configuration as configuration # import Config, dump_to_yaml - from feniax.preprocessor.inputs import Inputs - import feniax.feniax_main - import feniax.plotools.reconstruction as reconstruction - - label_gaf = "<>" - num_poles = 5 - Dhj_file = f"D{label_gaf}p{num_poles}" - Ahh_file = f"A{label_gaf}p{num_poles}" - Poles_file = f"Poles{label_gaf}p{num_poles}" - - inp = Inputs() - inp.engine = "intrinsicmodal" - inp.fem.eig_type = "inputs" - - inp.fem.connectivity = dict(# FusWing=['RWing', - # 'LWing'], - FusBack=['FusTail', - 'VTP'], - FusFront=None, - RWing=None, - LWing=None, - FusTail=None, - VTP=['HTP', 'VTPTail'], - HTP=['RHTP', 'LHTP'], - VTPTail=None, - RHTP=None, - LHTP=None, - ) - - inp.fem.grid = f"./FEM/structuralGrid_{sol[:-1]}" - #inp.fem.folder = pathlib.Path('./FEM/') - inp.fem.Ka_name = f"./FEM/Ka_{sol[:-1]}.npy" - inp.fem.Ma_name = f"./FEM/Ma_{sol[:-1]}.npy" - inp.fem.eig_names = [f"./FEM/eigenvals_{sol}{num_modes}.npy", - f"./FEM/eigenvecs_{sol}{num_modes}.npy"] - inp.fem.num_modes = 50 - inp.driver.typeof = "intrinsic" - - #inp.driver.sol_path = pathlib.Path( - # f"./results_{datetime.datetime.now().strftime('%Y-%m-%d_%H:%M:%S')}") - inp.driver.sol_path = pathlib.Path( - "./results1gust") +#+NAME: GustMain +#+begin_src python + <> + <> + <> inp.simulation.typeof = "single" inp.system.name = "s1" - inp.system.bc1 = 'free' + if sol[0] == "e": # free model, otherwise clamped + inp.system.bc1 = 'free' + inp.system.q0treatment = 1 inp.system.solution = "dynamic" - inp.system.t1 = 1 - inp.system.tn = 1001 + inp.system.t1 = 2 + inp.system.tn = 2001 inp.system.solver_library = "runge_kutta" inp.system.solver_function = "ode" inp.system.solver_settings = dict(solver_name="rk4") inp.system.xloads.modalaero_forces = True - inp.system.q0treatment = 1 - inp.system.aero.c_ref = <> - inp.system.aero.u_inf = 200. - inp.system.aero.rho_inf = 0.4 + inp.system.aero.c_ref = c_ref + inp.system.aero.u_inf = u_inf + inp.system.aero.rho_inf = rho_inf inp.system.aero.poles = f"./AERO/{Poles_file}.npy" inp.system.aero.A = f"./AERO/{Ahh_file}.npy" inp.system.aero.D = f"./AERO/{Dhj_file}.npy" inp.system.aero.gust_profile = "mc" - inp.system.aero.gust.intensity = 10*2 - inp.system.aero.gust.length = 100. + inp.system.aero.gust.intensity = 20 + inp.system.aero.gust.length = 150. inp.system.aero.gust.step = 0.1 inp.system.aero.gust.shift = 0. - inp.system.aero.gust.panels_dihedral = jnp.load("./AERO/Dihedral_d1c7.npy") - inp.system.aero.gust.collocation_points = "./AERO/Collocation_d1c7.npy" + inp.system.aero.gust.panels_dihedral = f"./AERO/Dihedral_{label_dlm}.npy" + inp.system.aero.gust.collocation_points = f"./AERO/Collocation_{label_dlm}.npy" +#+end_src + +*** Run +#+begin_src python :tangle settings_gust1.py + <> + <> + inp.driver.sol_path = pathlib.Path( + f"./results/gust2_{sol}") config = configuration.Config(inp) + t1 = time.time() sol1 = feniax.feniax_main.main(input_obj=config) -#+end_src + t2 = time.time() + print(f"Time gust: {t2 - t1}") -#+RESULTS: +#+end_src *** Plot #+NAME: Tipdisplacement @@ -2035,7 +2101,7 @@ Running Nastran using the tailored functions in run_nastra.sh which moves output # dash=dash[i % 3]) ), - dict(xaxis_range=[0,1],yaxis_range=[0,5])) + dict(xaxis_range=[0,2],yaxis_range=[0,5])) figname = f"./figs/try1.png" fig.write_image(figname, scale=1) fig.show() @@ -2091,9 +2157,9 @@ Running Nastran using the tailored functions in run_nastra.sh which moves output 'RBE3.vtu', 'CBAR.vtu'] - folder = pathlib.Path('./paraview/soldyn1') + folder = pathlib.Path("./results1gust/paraview/") # pathlib.Path('./paraview/soldyn1') folder_out = folder / "merged" - folder_out.mkdir(exist_ok=True) + folder_out.mkdir(exist_ok=True, parents=True) for i, fi in enumerate(folder.glob("bug_*")): if fi.is_dir() and (fi / paraview_files[0]).is_file(): #print(fi) @@ -2114,17 +2180,677 @@ Running Nastran using the tailored functions in run_nastra.sh which moves output #+end_src -#+begin_src bash :results none +merge all vtu files +#+begin_src bash :results none :session shell4 ~/Downloads/ParaView-5.10.1-MPI-Linux-Python3.9-x86_64/bin/pvpython paraview_3Ddynamic.py #+end_src -#+begin_src bash :session shell3 - cd ./paraview/soldyn1/merged - ffmpeg -framerate 20 -pattern_type glob -i '*.png' -c:v mpeg4 -qscale:v 1 BugGust.mp4 +make gust video +#+begin_src bash :session shell4 + cd ./results1gust/paraview/merged + # ffmpeg -y -framerate 20 -pattern_type glob -i '*.png' -c:v mpeg4 -qscale:v 1 BugGust.mp4 + ffmpeg -y -framerate 20 -pattern_type glob -i '*.png' -c:v mjpeg -qscale:v 1 BugGust.mp4 + cd - +#+end_src + +make gif +#+begin_src bash :session shell4 + cd ./results1gust/paraview/merged + ffmpeg -i BugGust.mp4 -filter_complex "[0:v] split [a][b];[a] palettegen [p];[b][p] paletteuse" BugGust.gif cd - #+end_src + +#+begin_src bash :session shell4 + cd ./results1gust/paraview/merged + # ffmpeg -y -framerate 20 -pattern_type glob -i '*.png' -c:v mpeg4 -qscale:v 1 BugGust.mp4 + ffmpeg -y -framerate 20 -pattern_type glob -i '*.png' -c:v mjpeg -qscale:v 1 BugGust.mp4 + cd - +#+end_src + + #+RESULTS: +: bash: cd: ./results1gust/paraview/merged: No such file or directory + +#+begin_src bash :session shell4 + cd ./paraview/soldyn1/merged + # ffmpeg -y -framerate 20 -pattern_type glob -i '*.png' BugGustWrecked.gif + #ffmpeg -y -framerate 20 -pattern_type glob -i '*.png' -filter_complex "[0:v] split [a][b];[a] palettegen [p];[b][p] paletteuse" BugGustWrecked2.gif + ffmpeg -i BugGust_wrecked.mp4 -filter_complex "[0:v] split [a][b];[a] palettegen [p];[b][p] paletteuse" BugGustWrecked2.gif + cd - +#+end_src + +** Sharding + +*** Discrete loads +Forces and moments at node 35 and 61 + +#+NAME: DiscreteLoads +#+begin_src python :tangle settings_ShardDiscreteLoads.py :session *pyshard3* + <> + sol = "cao" + num_modes = 300 + <> + inp.driver.sol_path = pathlib.Path( + "./results/ShardDiscreteLoads1") + + inp.simulation.typeof = "single" + inp.system.name = "s1" + inp.system.solution = "static" + inp.system.solver_library = "diffrax" + inp.system.solver_function = "newton" + inp.system.solver_settings = dict(rtol=1e-6, + atol=1e-6, + max_steps=100, + norm="linalg_norm", + kappa=0.01) + inp.system.xloads.follower_forces = True + inp.system.xloads.x = [0, 1, 2, 3, 4, 5] + inp.system.t = [0.5, 1, 1.5, 2, 2.5, 3, 3.25, 3.5, 3.75, 4, 4.25, 4.5, 4.75, 5.] + lz1 = 5e4 * 0.5 + lz2 = 9e4 * 0.5 + lz3 = 2e5 * 0.5 + lz4 = 4e5 * 0.5 + lz5 = 5e5 * 0.5 + lx1 = lz1 * 5 + lx2 = lz2 * 5 + lx3 = lz3 * 5 + lx4 = lz4 * 5 + lx5 = lz5 * 5 + ly1 = lz1 * 7 + ly2 = lz2 * 7 + ly3 = lz3 * 7 + ly4 = lz4 * 7 + ly5 = lz5 * 7 + + # rwing: 14-35 + # lwing: 40-61 + # [[[node_i, component_j]..(total_forces per run)],...(parallel forces)[[node_i, component_j]..]] + inputforces = dict(follower_points=[[[35, 0], [61, 0], [35, 1], [61, 1]], + [[35, 1], [61, 1], [35, 0], [61, 0]], + [[35, 2], [61, 2], [35, 0], [61, 0]], + [[35, 3], [61, 3], [35, 0], [61, 0]], + [[35, 4], [61, 4], [35, 0], [61, 0]], + [[35, 5], [61, 5], [35, 0], [61, 0]], + [[35, 0], [61, 0], [35, 4], [61, 4]], + [[35, 2], [61, 2], [35, 4], [61, 4]], + ], + # [[[0,...interpolation points]..(total_forces per run)],...(parallel forces)[[0,...]..]] + follower_interpolation= [[[0., lx1, lx2, lx3, lx4, lx5], + [0., lx1, lx2, lx3, lx4, lx5], + [0., 0., 0., 0., 0., 0.], + [0., 0., 0., 0., 0., 0.]], + [[0., ly1, ly2, ly3, ly4, ly5], + [0., ly1, ly2, ly3, ly4, ly5], + [0., 0., 0., 0., 0., 0.], + [0., 0., 0., 0., 0., 0.]], + [[0., lz1, lz2, lz3, lz4, lz5], + [0., lz1, lz2, lz3, lz4, lz5], + [0., 0., 0., 0., 0., 0.], + [0., 0., 0., 0., 0., 0.]], + [[0., lz1, lz2, lz3, lz4, lz5], + [0., lz1, lz2, lz3, lz4, lz5], + [0., 0., 0., 0., 0., 0.], + [0., 0., 0., 0., 0., 0.]], + [[0., lx1, lx2, lx3, lx4, lx5], + [0., lx1, lx2, lx3, lx4, lx5], + [0., 0., 0., 0., 0., 0.], + [0., 0., 0., 0., 0., 0.]], + [[0., lx1, lx2, lx3, lx4, lx5], + [0., lx1, lx2, lx3, lx4, lx5], + [0., 0., 0., 0., 0., 0.], + [0., 0., 0., 0., 0., 0.]], + [[0., lz1, lz2, lz3, lz4, lz5], + [0., lz1, lz2, lz3, lz4, lz5], + [0., lx1, lx2, lx3, lx4, lx5], + [0., lx1, lx2, lx3, lx4, lx5] + ], + [[0., lz1, lz2, lz3, lz4, lz5], + [0., lz1, lz2, lz3, lz4, lz5], + [0., lx1, lx2, lx3, lx4, lx5], + [0., lx1, lx2, lx3, lx4, lx5]] + ] + ) + inp.system.shard = dict(input_type="pointforces", + inputs=inputforces) + t1 = time.time() + sol = feniax.feniax_shardmain.main(input_dict=inp, device_count=device_count) + t2 = time.time() + print(f"Time DiscreteLoads: {t2 - t1}") +#+end_src + +**** Plot results + +#+begin_src python :session *pyshard3* + import pyvista as pv + import feniax.intrinsic.xloads as xloads + import importlib + importlib.reload(xloads) + import pathlib + import feniax.preprocessor.configuration as configuration + def glyph_forcefolllower(t, x, force, ra, Rab, scale=1, save_path=None, folder="", name="force", X=None, R0ab=None, plot=False): + ra = np.array(ra) + Rab = np.array(Rab) + if save_path is not None: + save_path = pathlib.Path(save_path) + for i, ti in enumerate(t): + _ffollower = np.array(xloads.force_pointfollower(t=ti, x=np.array(inp.system.xloads.x), force_follower=force, Rab=Rab[i])) + ffollower = _ffollower[:3].T / scale + mesh = pv.PolyData(ra[i].T) + mesh["vectors"] = ffollower + # Create glyphs to represent vectors + glyphs = mesh.glyph() #orient="vectors", scale=False, factor=1e-3) + if save_path is not None: + path = save_path / f"{folder}_{i}.vtk" + path.mkdir(parents=True, exist_ok=True) + glyphs.save(path / f"{name}.vtk", binary=False) + if X is not None and R0ab is not None: + _ffollower = np.array(xloads.force_pointfollower(t=0, x=np.array(inp.system.xloads.x), + force_follower=force, Rab=np.array(R0ab))) + ffollower = _ffollower[:3].T / scale + mesh = pv.PolyData(np.array(X)) + mesh["vectors"] = ffollower + # Create glyphs to represent vectors + glyphs = mesh.glyph() #orient="vectors", scale=False, factor=1e-3) + if save_path is not None: + path = save_path / f"{folder}_ref.vtk/" + path.mkdir(parents=True, exist_ok=True) + glyphs.save(path / f"{name}.vtk", binary=False) + if plot: + # Plot the vector field + plotter = pv.Plotter() + plotter.add_mesh(glyphs, color='red') + plotter.add_mesh(mesh) + plotter.show() + + config = configuration.Config(inp) + X=config.fem.X + for li in range(8): + f = sol.shards_sys1.points[li] + R0ab=sol.modes.C0ab + glyph_forcefolllower(t=inp.system.t, + x = inp.system.xloads.x, + force=f, + ra=sol.staticsystem_sys1.ra[li], + Rab=sol.staticsystem_sys1.Cab[li], + scale=1e4*4, + save_path=inp.driver.sol_path / f"paraviewM{num_modes}/L{li}/", + folder="bug", + X=X, + R0ab=R0ab) +#+end_src + +#+begin_src python :session *pyshard3* + import feniax.plotools.reconstruction as reconstruction + import feniax.preprocessor.configuration as configuration + config = configuration.Config(inp) + for i in range(8): + rintrinsic, uintrinsic = reconstruction.rbf_based( + nastran_bdf="./NASTRAN/BUG103.bdf", + X=config.fem.X, + time=range(len(inp.system.t)), + ra=sol.staticsystem_sys1.ra[i], + Rab=sol.staticsystem_sys1.Cab[i], + R0ab=sol.modes.C0ab, + vtkpath=inp.driver.sol_path / f"paraviewM{num_modes}/L{i}/bug", + plot_timeinterval=1, + plot_ref=False, + tolerance=1e-3, + size_cards=8, + rbe3s_full=False, + ra_movie=None) +#+end_src + +#+begin_src python :tangle ./results/ShardDiscreteLoads1/paraview/merge.py :session *pyshard3* + + from paraview.simple import * + import pathlib + + def merge_paraview(file_list, file_out): + # Create a reader for each file + #readers = [XMLUnstructuredGridReader(FileName=file) for file in file_list] + readers = [OpenDataFile(file) for file in file_list] + # Append the readers to merge the datasets + appended = AppendDatasets(Input=readers) + + # Save the merged dataset + #writer = XMLUnstructuredGridWriter(Input=appended, FileName=file_out) + #writer.UpdatePipeline() + SaveData(file_out, proxy=appended) + + paraview_files = ['CQUAD4.vtu', + 'CONM2.vtu', + #'CBUSH.vtu', + 'RBE2.vtu', + 'CTRIA3.vtu', + 'RBE3.vtu', + 'CBAR.vtu', + 'force.vtk' + ] + + for i in range(8): + folder = pathlib.Path(f"./L{i}") + folder_out = folder / "merged" + folder_out.mkdir(exist_ok=True, parents=True) + for i, fi in enumerate(folder.glob("bug_*")): + if fi.is_dir() and (fi / paraview_files[0]).is_file(): + print(fi) + name_len = len("bug_") + index = fi.name[name_len:].split('.')[0] + file_list = [str(fi / pvf) for pvf in paraview_files] + #print(file_list) + file_out = str(folder_out / f"bug_{index}.vtu") + # readers = [XMLUnstructuredGridReader(FileName=file) for file in file_list] + + # # Append the readers to merge the datasets + # appended = AppendDatasets(Input=readers) + # print(file_out) + # # Save the merged dataset + # writer = XMLUnstructuredGridWriter(Input=appended, FileName=file_out) + # writer.UpdatePipeline() + merge_paraview(file_list, file_out) + +#+end_src + +merge all vtu files +#+begin_src bash :results none :session shell1 + # cd results/ShardDiscreteLoads1/paraview/ + ~/Downloads/ParaView-5.10.1-MPI-Linux-Python3.9-x86_64/bin/pvpython merge.py +#+end_src +one needs to go to paraview now and save images + +make gust video +#+begin_src bash :session shell1 + # ffmpeg -y -framerate 20 -pattern_type glob -i '*.png' -c:v mpeg4 -qscale:v 1 BugGust.mp4 + ffmpeg -y -framerate 1 -pattern_type glob -i '*.png' -c:v mjpeg -qscale:v 1 BugGust.mp4 + cd - +#+end_src + +*** Discrete loads Montecarlo +**** High loading +#+begin_src python :tangle settings_DiscreteLoadsMC.py :session *pyshard4* + <> + sol = "cao" + paths = 8*200 + num_modes = 100 + device_count = 1 + <> + inp.driver.sol_path = pathlib.Path( + f"./results/DiscreteMC1high{num_modes}") + + inp.simulation.typeof = "single" + inp.system.name = "s1" + inp.system.solution = "static" + inp.system.solver_library = "diffrax" + inp.system.solver_function = "newton" + inp.system.solver_settings = dict(rtol=1e-6, + atol=1e-6, + max_steps=50, + norm="linalg_norm", + kappa=0.01) + inp.system.xloads.follower_forces = True + inp.system.xloads.x = [0, 1, 2, 3, 4, 5] + inp.system.t = [1, 2, 3, 4, 5] + # rwing: 14-35 + # lwing: 40-61 + points = [] + interpolation = [] + _interpolation = [0., 3.e3, 7e3, 9e3, 1e4, 1.5e4] # 1.5e4, 2e4..4e4] #[0., 0., 0., 0.] + _interpolation_torsion = [i*2 for i in _interpolation] #[0., 4e3, 1e4, 2e4, 4e4, 5e4] + for ri,li in zip(range(14, 36),range(40,62)): + points.append([ri, 2]) + points.append([ri, 4]) + points.append([li, 2]) + points.append([li, 4]) + for i, _ in enumerate(range(len(points))): + + if i % 2 == 0: + interpolation.append(_interpolation) + else: + interpolation.append(_interpolation_torsion) + + interpolation = np.array(interpolation) # num_pointforces x num_interpolation + sigma0 = 0.15 # percentage of mu for sigma + mu = _interpolation[-1] + sigma = (sigma0) * _interpolation[-1] + rand = np.random.normal(mu, sigma, paths) + mu_torsion = _interpolation_torsion[-1] + sigma_torsion = (sigma0) * _interpolation_torsion[-1] + rand_torsion = np.random.normal(mu_torsion, sigma_torsion, paths) + follower_interpolation = [] + for i, ri in enumerate(rand): + interpolationrand = np.copy(interpolation) + interpolationrand[::2, -1] = ri + interpolationrand[1::2, -1] = rand_torsion[i] + follower_interpolation.append(interpolationrand) + #follower_interpolation = [interpolation * ri for ri in rand] + follower_points = [points]*paths + inputforces = dict(follower_points=follower_points, + follower_interpolation=follower_interpolation + ) + inp.system.shard = dict(input_type="pointforces", + inputs=inputforces) + + t1 = time.time() + sol1 = feniax.feniax_shardmain.main(input_dict=inp, device_count=device_count) + t2 = time.time() + print(f"Time DiscreteLoads MC: {t2 - t1}") + + # np.mean(sol.staticsystem_sys1.ra[:,-1,2,35]) + # np.std(sol.staticsystem_sys1.ra[:,-1,2,35]) +#+end_src + +***** Plot results + +#+begin_src python :session *pyshard4* + import feniax.plotools.reconstruction as reconstruction + import feniax.preprocessor.configuration as configuration + config = configuration.Config(inp) + for li in [0, paths-1]: + rintrinsic, uintrinsic = reconstruction.rbf_based( + nastran_bdf="./NASTRAN/BUG103.bdf", + X=config.fem.X, + time=range(len(inp.system.t)), + ra=sol1.staticsystem_sys1.ra[li], + Rab=sol1.staticsystem_sys1.Cab[li], + R0ab=sol1.modes.C0ab, + vtkpath=inp.driver.sol_path / f"paraviewM{num_modes}/L{li}/bug", + plot_timeinterval=1, + plot_ref=True, + tolerance=1e-3, + size_cards=8, + rbe3s_full=False, + ra_movie=None + ) +#+end_src + + +#+NAME: discrete loads +#+begin_src python :session *pyshard4* + import pyvista as pv + import feniax.intrinsic.xloads as xloads + import importlib + importlib.reload(xloads) + import pathlib + import feniax.preprocessor.configuration as configuration + def glyph_forcefolllower(t, x, force, ra, Rab, scale=1, save_path=None, folder="", name="force", X=None, R0ab=None, plot=False): + ra = np.array(ra) + Rab = np.array(Rab) + if save_path is not None: + save_path = pathlib.Path(save_path) + for i, ti in enumerate(t): + _ffollower = np.array(xloads.force_pointfollower(t=ti, x=np.array(inp.system.xloads.x), force_follower=force, Rab=Rab[i])) + ffollower = _ffollower[:3].T / scale + mesh = pv.PolyData(ra[i].T) + mesh["vectors"] = ffollower + # Create glyphs to represent vectors + glyphs = mesh.glyph() #orient="vectors", scale=False, factor=1e-3) + if save_path is not None: + path = save_path / f"{folder}_{i}.vtk" + path.mkdir(parents=True, exist_ok=True) + glyphs.save(path / f"{name}.vtk", binary=False) + if X is not None and R0ab is not None: + _ffollower = np.array(xloads.force_pointfollower(t=0, x=np.array(inp.system.xloads.x), + force_follower=force, Rab=np.array(R0ab))) + ffollower = _ffollower[:3].T / scale + mesh = pv.PolyData(np.array(X)) + mesh["vectors"] = ffollower + # Create glyphs to represent vectors + glyphs = mesh.glyph() #orient="vectors", scale=False, factor=1e-3) + if save_path is not None: + path = save_path / f"{folder}_ref.vtk/" + path.mkdir(parents=True, exist_ok=True) + glyphs.save(path / f"{name}.vtk", binary=False) + if plot: + # Plot the vector field + plotter = pv.Plotter() + plotter.add_mesh(glyphs, color='red') + plotter.add_mesh(mesh) + plotter.show() + + config = configuration.Config(inp) + X=config.fem.X + for li in [0, paths-1]: + f = sol1.shards_sys1.points[li] + R0ab=sol1.modes.C0ab + glyph_forcefolllower(t=inp.system.t, + x = inp.system.xloads.x, + force=f, + ra=sol1.staticsystem_sys1.ra[li], + Rab=sol1.staticsystem_sys1.Cab[li], + scale=1e4*0.5, + save_path=inp.driver.sol_path / f"paraviewM{num_modes}/L{li}/", + folder="bug", + X=X, + R0ab=R0ab) +#+end_src + +**** Small loading + +#+begin_src python :tangle settings_DiscreteLoadsMC.py :session *pyshard4* + <> + sol = "cao" + <> + inp.driver.sol_path = pathlib.Path( + "./results/DiscreteMC1small") + + inp.simulation.typeof = "single" + inp.system.name = "s1" + inp.system.solution = "static" + inp.system.solver_library = "diffrax" + inp.system.solver_function = "newton" + inp.system.solver_settings = dict(rtol=1e-6, + atol=1e-6, + max_steps=50, + norm="linalg_norm", + kappa=0.01) + inp.system.xloads.follower_forces = True + inp.system.xloads.x = [0, 1, 2, 3, 4, 5] + inp.system.t = [1, 2, 3, 4, 5] + # rwing: 14-35 + # lwing: 40-61 + points = [] + interpolation = [] + _interpolationsmall = [i*1e-2 for i in _interpolation] + _interpolationsmall_torsion = [i*1e-2 for i in _interpolation_torsion] + for ri,li in zip(range(14, 36),range(40,62)): + points.append([ri, 2]) + points.append([ri, 4]) + points.append([li, 2]) + points.append([li, 4]) + for i, _ in enumerate(range(len(points))): + + if i % 2 == 0: + interpolation.append(_interpolationsmall) + else: + interpolation.append(_interpolationsmall_torsion) + + interpolation = np.array(interpolation) # num_pointforces x num_interpolation + sigma0 = 0.15 # percentage of mu for sigma + mu = _interpolationsmall[-1] + sigma = (sigma0) * _interpolationsmall[-1] + rand = np.random.normal(mu, sigma, paths) + mu_torsion = _interpolationsmall_torsion[-1] + sigma_torsion = (sigma0) * _interpolationsmall_torsion[-1] + rand_torsion = np.random.normal(mu_torsion, sigma_torsion, paths) + follower_interpolation = [] + for i, ri in enumerate(rand): + interpolationrand = np.copy(interpolation) + interpolationrand[::2, -1] = ri + interpolationrand[1::2, -1] = rand_torsion[i] + follower_interpolation.append(interpolationrand) + #follower_interpolation = [interpolation * ri for ri in rand] + follower_points = [points]*paths + inputforces = dict(follower_points=follower_points, + follower_interpolation=follower_interpolation + ) + inp.system.shard = dict(input_type="pointforces", + inputs=inputforces) + t1 = time.time() + sol2 = feniax.feniax_shardmain.main(input_dict=inp, device_count=device_count) + t2 = time.time() + print(f"Time DiscreteLoads MC: {t2 - t1}") + # np.mean(sol.staticsystem_sys1.ra[:,-1,2,35]) + # np.std(sol.staticsystem_sys1.ra[:,-1,2,35]) +#+end_src + +**** Very Small loading + #+begin_src python :tangle settings_DiscreteLoadsMC.py :session *pyshard4* + <> + sol = "cao" + <> + inp.driver.sol_path = pathlib.Path( + "./results/DiscreteMC1verysmall") + + inp.simulation.typeof = "single" + inp.system.name = "s1" + inp.system.solution = "static" + inp.system.solver_library = "diffrax" + inp.system.solver_function = "newton" + inp.system.solver_settings = dict(rtol=1e-6, + atol=1e-6, + max_steps=50, + norm="linalg_norm", + kappa=0.01) + inp.system.xloads.follower_forces = True + inp.system.xloads.x = [0, 1, 2, 3, 4, 5] + inp.system.t = [1, 2, 3, 4, 5] + # rwing: 14-35 + # lwing: 40-61 + points = [] + interpolation = [] + _interpolationverysmall = [i*1e-3 for i in _interpolation] + _interpolationverysmall_torsion = [i*1e-3 for i in _interpolation_torsion] + + for ri,li in zip(range(14, 36),range(40,62)): + points.append([ri, 2]) + points.append([ri, 4]) + points.append([li, 2]) + points.append([li, 4]) + for i, _ in enumerate(range(len(points))): + + if i % 2 == 0: + interpolation.append(_interpolationverysmall) + else: + interpolation.append(_interpolationverysmall_torsion) + + interpolation = np.array(interpolation) # num_pointforces x num_interpolationverysmall + sigma0 = 0.15 # percentage of mu for sigma + mu = _interpolationverysmall[-1] + sigma = (sigma0) * _interpolationverysmall[-1] + rand = np.random.normal(mu, sigma, paths) + mu_torsion = _interpolationverysmall_torsion[-1] + sigma_torsion = (sigma0) * _interpolationverysmall_torsion[-1] + rand_torsion = np.random.normal(mu_torsion, sigma_torsion, paths) + follower_interpolation = [] + for i, ri in enumerate(rand): + interpolationrand = np.copy(interpolation) + interpolationrand[::2, -1] = ri + interpolationrand[1::2, -1] = rand_torsion[i] + follower_interpolation.append(interpolationrand) + #follower_interpolation = [interpolation * ri for ri in rand] + follower_points = [points]*paths + inputforces = dict(follower_points=follower_points, + follower_interpolation=follower_interpolation + ) + inp.system.shard = dict(input_type="pointforces", + inputs=inputforces) + t1 = time.time() + sol3 = feniax.feniax_shardmain.main(input_dict=inp, device_count=device_count) + t2 = time.time() + print(f"Time DiscreteLoads MC: {t2 - t1}") + + # np.mean(sol.staticsystem_sys1.ra[:,-1,2,35]) + # np.std(sol.staticsystem_sys1.ra[:,-1,2,35]) +#+end_src +**** Statistics +#+begin_src python :tangle settings_DiscreteLoadsMC.py :session *pyshard4* + u_mean = np.mean(sol1.staticsystem_sys1.ra[:,-1,2,35] - config.fem.X[35,2]) + u_std = np.std(sol1.staticsystem_sys1.ra[:,-1,2,35]) + + print(f"Mean displacement node 35: {u_mean}") + print(f"std displacement node 35: {u_std}") + print(f"Ratio displacement node 35: {u_mean/u_std}") + print("***************") + + u_mean2 = np.mean(sol2.staticsystem_sys1.ra[:,-1,2,35] - config.fem.X[35,2]) + u_std2 = np.std(sol2.staticsystem_sys1.ra[:,-1,2,35]) + + print(f"Mean displacement node 35: {u_mean2}") + print(f"std displacement node 35: {u_std2}") + print(f"ratio displacement node 35: {u_mean2/u_std2}") + print("***************") + + u_mean3 = np.mean(sol3.staticsystem_sys1.ra[:,-1,2,35] - config.fem.X[35,2]) + u_std3 = np.std(sol3.staticsystem_sys1.ra[:,-1,2,35]) + + print(f"Mean displacement node 35: {u_mean3}") + print(f"std displacement node 35: {u_std3}") + print(f"ratio displacement node 35: {u_mean3/u_std3}") + print("***************") + +#+end_src + +*** Manoeuvre +#+begin_src python :tangle ./settings_manoeuvre1shard.py :session *pyshard1* + <> + <> + device_count = 8 + #rho_rand = np.random.normal(0.6, 0.6*0.15, 500) + inp.driver.sol_path = pathlib.Path( + "./results/manoeuvre1Shard") + inputflow = dict(u_inf=np.linspace(167.7, 251.6, 16), + rho_inf=np.linspace(0.33, 0.5, 16)) + inp.system.shard = dict(input_type="steadyalpha", + inputs=inputflow) + t1 = time.time() + solstatic1shard = feniax.feniax_shardmain.main(input_dict=inp, device_count=device_count) + t2 = time.time() + print(f"Time Manoeuvre: {t2 - t1}") + +#+end_src + +*** Gust +for local testing: +#+begin_src python :session *pyshard2* + <> + <> + device_count = 4 + inp.driver.sol_path = pathlib.Path( + f"./results/gust2_{sol}Shard") + inp.system.aero.gust.fixed_discretisation = [150, u_inf] + # Shard inputs + inputflow = dict(length=np.linspace(25,265,8), + intensity= np.linspace(0.1, 20, 8), + rho_inf = np.linspace(0.3, 0.48, 8) + ) + inp.system.shard = dict(input_type="gust1", + inputs=inputflow) + t1 = time.time() + solgust21shard = feniax.feniax_shardmain.main(input_dict=inp, device_count=device_count) + t2 = time.time() + print(f"Time SHARD Gust: {t2 - t1}") + +#+end_src + +#+begin_src python :tangle settings_gust1shard.py :session *pyshard2* + <> + <> + + inp.driver.sol_path = pathlib.Path( + f"./results/gust2_{sol}Shard") + inp.system.aero.gust.fixed_discretisation = [150, u_inf] + # Shard inputs + inputflow = dict(length=np.linspace(25,265,13), + intensity=np.linspace(0.1, 3, 11), + rho_inf = np.linspace(0.34,0.48,8) + ) + inp.system.shard = dict(input_type="gust1", + inputs=inputflow) + + solgust21shard = feniax.feniax_shardmain.main(input_dict=inp, device_count=device_count) +#+end_src + * Gust solution (146) :PROPERTIES: diff --git a/examples/BUG/settings_DiscreteLoads.py b/examples/BUG/settings_DiscreteLoads.py new file mode 100644 index 0000000..7992b10 --- /dev/null +++ b/examples/BUG/settings_DiscreteLoads.py @@ -0,0 +1,123 @@ +# [[file:modelgeneration.org::DiscreteLoads][DiscreteLoads]] +import pathlib +import time +#import jax.numpy as jnp +import numpy as np +from feniax.preprocessor.inputs import Inputs +import feniax.feniax_shardmain +sol = "cao" +num_modes = 300 +inp = Inputs() +inp.engine = "intrinsicmodal" +inp.fem.eig_type = "inputs" + +inp.fem.connectivity = dict(# FusWing=['RWing', + # 'LWing'], + FusBack=['FusTail', + 'VTP'], + FusFront=None, + RWing=None, + LWing=None, + FusTail=None, + VTP=['HTP', 'VTPTail'], + HTP=['RHTP', 'LHTP'], + VTPTail=None, + RHTP=None, + LHTP=None, + ) +inp.fem.grid = f"./FEM/structuralGrid_{sol[:-1]}" +#inp.fem.folder = pathlib.Path('./FEM/') +inp.fem.Ka_name = f"./FEM/Ka_{sol[:-1]}.npy" +inp.fem.Ma_name = f"./FEM/Ma_{sol[:-1]}.npy" +inp.fem.eig_names = [f"./FEM/eigenvals_{sol}{num_modes}.npy", + f"./FEM/eigenvecs_{sol}{num_modes}.npy"] +inp.driver.typeof = "intrinsic" +inp.fem.num_modes = num_modes +inp.driver.typeof = "intrinsic" +inp.driver.sol_path = pathlib.Path( + "./results/DiscreteLoads1") + +inp.simulation.typeof = "single" +inp.system.name = "s1" +inp.system.solution = "static" +inp.system.solver_library = "diffrax" +inp.system.solver_function = "newton" +inp.system.solver_settings = dict(rtol=1e-6, + atol=1e-6, + max_steps=100, + norm="linalg_norm", + kappa=0.01) +inp.system.xloads.follower_forces = True +inp.system.xloads.x = [0, 1, 2, 3, 4, 5] +inp.system.t = [0.5, 1, 1.5, 2, 2.5, 3, 3.25, 3.5, 3.75, 4, 4.25, 4.5, 4.75, 5.] +lz1 = 5e4 * 0.5 +lz2 = 9e4 * 0.5 +lz3 = 2e5 * 0.5 +lz4 = 4e5 * 0.5 +lz5 = 5e5 * 0.5 +lx1 = lz1 * 5 +lx2 = lz2 * 5 +lx3 = lz3 * 5 +lx4 = lz4 * 5 +lx5 = lz5 * 5 +ly1 = lz1 * 7 +ly2 = lz2 * 7 +ly3 = lz3 * 7 +ly4 = lz4 * 7 +ly5 = lz5 * 7 + +# rwing: 14-35 +# lwing: 40-61 + # [[[node_i, component_j]..(total_forces per run)],...(parallel forces)[[node_i, component_j]..]] +inputforces = dict(follower_points=[[[35, 0], [61, 0], [35, 1], [61, 1]], + [[35, 1], [61, 1], [35, 0], [61, 0]], + [[35, 2], [61, 2], [35, 0], [61, 0]], + [[35, 3], [61, 3], [35, 0], [61, 0]], + [[35, 4], [61, 4], [35, 0], [61, 0]], + [[35, 5], [61, 5], [35, 0], [61, 0]], + [[35, 0], [61, 0], [35, 4], [61, 4]], + [[35, 2], [61, 2], [35, 4], [61, 4]], + ], + # [[[0,...interpolation points]..(total_forces per run)],...(parallel forces)[[0,...]..]] + follower_interpolation= [[[0., lx1, lx2, lx3, lx4, lx5], + [0., lx1, lx2, lx3, lx4, lx5], + [0., 0., 0., 0., 0., 0.], + [0., 0., 0., 0., 0., 0.]], + [[0., ly1, ly2, ly3, ly4, ly5], + [0., ly1, ly2, ly3, ly4, ly5], + [0., 0., 0., 0., 0., 0.], + [0., 0., 0., 0., 0., 0.]], + [[0., lz1, lz2, lz3, lz4, lz5], + [0., lz1, lz2, lz3, lz4, lz5], + [0., 0., 0., 0., 0., 0.], + [0., 0., 0., 0., 0., 0.]], + [[0., lz1, lz2, lz3, lz4, lz5], + [0., lz1, lz2, lz3, lz4, lz5], + [0., 0., 0., 0., 0., 0.], + [0., 0., 0., 0., 0., 0.]], + [[0., lx1, lx2, lx3, lx4, lx5], + [0., lx1, lx2, lx3, lx4, lx5], + [0., 0., 0., 0., 0., 0.], + [0., 0., 0., 0., 0., 0.]], + [[0., lx1, lx2, lx3, lx4, lx5], + [0., lx1, lx2, lx3, lx4, lx5], + [0., 0., 0., 0., 0., 0.], + [0., 0., 0., 0., 0., 0.]], + [[0., lz1, lz2, lz3, lz4, lz5], + [0., lz1, lz2, lz3, lz4, lz5], + [0., lx1, lx2, lx3, lx4, lx5], + [0., lx1, lx2, lx3, lx4, lx5] + ], + [[0., lz1, lz2, lz3, lz4, lz5], + [0., lz1, lz2, lz3, lz4, lz5], + [0., lx1, lx2, lx3, lx4, lx5], + [0., lx1, lx2, lx3, lx4, lx5]] + ] + ) +inp.system.shard = dict(input_type="pointforces", + inputs=inputforces) +t1 = time.time() +sol = feniax.feniax_shardmain.main(input_dict=inp, device_count=8) +t2 = time.time() +print(f"Time DiscreteLoads: {t2 - t1}") +# DiscreteLoads ends here diff --git a/examples/BUG/settings_DiscreteLoadsMC.py b/examples/BUG/settings_DiscreteLoadsMC.py new file mode 100644 index 0000000..5785d6a --- /dev/null +++ b/examples/BUG/settings_DiscreteLoadsMC.py @@ -0,0 +1,324 @@ +# [[file:modelgeneration.org::*High loading][High loading:1]] +import pathlib +import time +#import jax.numpy as jnp +import numpy as np +from feniax.preprocessor.inputs import Inputs +import feniax.feniax_shardmain +sol = "cao" +num_modes = 100 +inp = Inputs() +inp.engine = "intrinsicmodal" +inp.fem.eig_type = "inputs" + +inp.fem.connectivity = dict(# FusWing=['RWing', + # 'LWing'], + FusBack=['FusTail', + 'VTP'], + FusFront=None, + RWing=None, + LWing=None, + FusTail=None, + VTP=['HTP', 'VTPTail'], + HTP=['RHTP', 'LHTP'], + VTPTail=None, + RHTP=None, + LHTP=None, + ) +inp.fem.grid = f"./FEM/structuralGrid_{sol[:-1]}" +#inp.fem.folder = pathlib.Path('./FEM/') +inp.fem.Ka_name = f"./FEM/Ka_{sol[:-1]}.npy" +inp.fem.Ma_name = f"./FEM/Ma_{sol[:-1]}.npy" +inp.fem.eig_names = [f"./FEM/eigenvals_{sol}{num_modes}.npy", + f"./FEM/eigenvecs_{sol}{num_modes}.npy"] +inp.driver.typeof = "intrinsic" +inp.fem.num_modes = num_modes +inp.driver.typeof = "intrinsic" +inp.driver.sol_path = pathlib.Path( + f"./results/DiscreteMC1high{num_modes}") + +inp.simulation.typeof = "single" +inp.system.name = "s1" +inp.system.solution = "static" +inp.system.solver_library = "diffrax" +inp.system.solver_function = "newton" +inp.system.solver_settings = dict(rtol=1e-6, + atol=1e-6, + max_steps=50, + norm="linalg_norm", + kappa=0.01) +inp.system.xloads.follower_forces = True +inp.system.xloads.x = [0, 1, 2, 3, 4, 5] +inp.system.t = [1, 2, 3, 4, 5] +# rwing: 14-35 +# lwing: 40-61 +points = [] +interpolation = [] +_interpolation = [0., 3.e3, 7e3, 9e3, 1e4, 1.5e4] # 1.5e4, 2e4..4e4] #[0., 0., 0., 0.] +_interpolation_torsion = [i*2 for i in _interpolation] #[0., 4e3, 1e4, 2e4, 4e4, 5e4] +for ri,li in zip(range(14, 36),range(40,62)): + points.append([ri, 2]) + points.append([ri, 4]) + points.append([li, 2]) + points.append([li, 4]) +for i, _ in enumerate(range(len(points))): + + if i % 2 == 0: + interpolation.append(_interpolation) + else: + interpolation.append(_interpolation_torsion) + +interpolation = np.array(interpolation) # num_pointforces x num_interpolation +paths = 8*200 +sigma0 = 0.15 # percentage of mu for sigma +mu = _interpolation[-1] +sigma = (sigma0) * _interpolation[-1] +rand = np.random.normal(mu, sigma, paths) +mu_torsion = _interpolation_torsion[-1] +sigma_torsion = (sigma0) * _interpolation_torsion[-1] +rand_torsion = np.random.normal(mu_torsion, sigma_torsion, paths) +follower_interpolation = [] +for i, ri in enumerate(rand): + interpolationrand = np.copy(interpolation) + interpolationrand[::2, -1] = ri + interpolationrand[1::2, -1] = rand_torsion[i] + follower_interpolation.append(interpolationrand) +#follower_interpolation = [interpolation * ri for ri in rand] +follower_points = [points]*paths +inputforces = dict(follower_points=follower_points, + follower_interpolation=follower_interpolation + ) +inp.system.shard = dict(input_type="pointforces", + inputs=inputforces) + +t1 = time.time() +sol1 = feniax.feniax_shardmain.main(input_dict=inp, device_count=8) +t2 = time.time() +print(f"Time DiscreteLoads MC: {t2 - t1}") + +# np.mean(sol.staticsystem_sys1.ra[:,-1,2,35]) +# np.std(sol.staticsystem_sys1.ra[:,-1,2,35]) +# High loading:1 ends here + +# [[file:modelgeneration.org::*Small loading][Small loading:1]] +import pathlib +import time +#import jax.numpy as jnp +import numpy as np +from feniax.preprocessor.inputs import Inputs +import feniax.feniax_shardmain +sol = "cao" +num_modes = 50 +inp = Inputs() +inp.engine = "intrinsicmodal" +inp.fem.eig_type = "inputs" + +inp.fem.connectivity = dict(# FusWing=['RWing', + # 'LWing'], + FusBack=['FusTail', + 'VTP'], + FusFront=None, + RWing=None, + LWing=None, + FusTail=None, + VTP=['HTP', 'VTPTail'], + HTP=['RHTP', 'LHTP'], + VTPTail=None, + RHTP=None, + LHTP=None, + ) +inp.fem.grid = f"./FEM/structuralGrid_{sol[:-1]}" +#inp.fem.folder = pathlib.Path('./FEM/') +inp.fem.Ka_name = f"./FEM/Ka_{sol[:-1]}.npy" +inp.fem.Ma_name = f"./FEM/Ma_{sol[:-1]}.npy" +inp.fem.eig_names = [f"./FEM/eigenvals_{sol}{num_modes}.npy", + f"./FEM/eigenvecs_{sol}{num_modes}.npy"] +inp.driver.typeof = "intrinsic" +inp.fem.num_modes = num_modes +inp.driver.typeof = "intrinsic" +inp.driver.sol_path = pathlib.Path( + "./results/DiscreteMC1small") + +inp.simulation.typeof = "single" +inp.system.name = "s1" +inp.system.solution = "static" +inp.system.solver_library = "diffrax" +inp.system.solver_function = "newton" +inp.system.solver_settings = dict(rtol=1e-6, + atol=1e-6, + max_steps=50, + norm="linalg_norm", + kappa=0.01) +inp.system.xloads.follower_forces = True +inp.system.xloads.x = [0, 1, 2, 3, 4, 5] +inp.system.t = [1, 2, 3, 4, 5] +# rwing: 14-35 +# lwing: 40-61 +points = [] +interpolation = [] +_interpolationsmall = [i*1e-2 for i in _interpolation] +_interpolationsmall_torsion = [i*1e-2 for i in _interpolation_torsion] +for ri,li in zip(range(14, 36),range(40,62)): + points.append([ri, 2]) + points.append([ri, 4]) + points.append([li, 2]) + points.append([li, 4]) +for i, _ in enumerate(range(len(points))): + + if i % 2 == 0: + interpolation.append(_interpolationsmall) + else: + interpolation.append(_interpolationsmall_torsion) + +interpolation = np.array(interpolation) # num_pointforces x num_interpolation +paths = 8*10 +sigma0 = 0.15 # percentage of mu for sigma +mu = _interpolationsmall[-1] +sigma = (sigma0) * _interpolationsmall[-1] +rand = np.random.normal(mu, sigma, paths) +mu_torsion = _interpolationsmall_torsion[-1] +sigma_torsion = (sigma0) * _interpolationsmall_torsion[-1] +rand_torsion = np.random.normal(mu_torsion, sigma_torsion, paths) +follower_interpolation = [] +for i, ri in enumerate(rand): + interpolationrand = np.copy(interpolation) + interpolationrand[::2, -1] = ri + interpolationrand[1::2, -1] = rand_torsion[i] + follower_interpolation.append(interpolationrand) +#follower_interpolation = [interpolation * ri for ri in rand] +follower_points = [points]*paths +inputforces = dict(follower_points=follower_points, + follower_interpolation=follower_interpolation + ) +inp.system.shard = dict(input_type="pointforces", + inputs=inputforces) + +sol2 = feniax.feniax_shardmain.main(input_dict=inp, device_count=8) +# np.mean(sol.staticsystem_sys1.ra[:,-1,2,35]) +# np.std(sol.staticsystem_sys1.ra[:,-1,2,35]) +# Small loading:1 ends here + +# [[file:modelgeneration.org::*Very Small loading][Very Small loading:1]] +import pathlib +import time +#import jax.numpy as jnp +import numpy as np +from feniax.preprocessor.inputs import Inputs +import feniax.feniax_shardmain +sol = "cao" +num_modes = 50 +inp = Inputs() +inp.engine = "intrinsicmodal" +inp.fem.eig_type = "inputs" + +inp.fem.connectivity = dict(# FusWing=['RWing', + # 'LWing'], + FusBack=['FusTail', + 'VTP'], + FusFront=None, + RWing=None, + LWing=None, + FusTail=None, + VTP=['HTP', 'VTPTail'], + HTP=['RHTP', 'LHTP'], + VTPTail=None, + RHTP=None, + LHTP=None, + ) +inp.fem.grid = f"./FEM/structuralGrid_{sol[:-1]}" +#inp.fem.folder = pathlib.Path('./FEM/') +inp.fem.Ka_name = f"./FEM/Ka_{sol[:-1]}.npy" +inp.fem.Ma_name = f"./FEM/Ma_{sol[:-1]}.npy" +inp.fem.eig_names = [f"./FEM/eigenvals_{sol}{num_modes}.npy", + f"./FEM/eigenvecs_{sol}{num_modes}.npy"] +inp.driver.typeof = "intrinsic" +inp.fem.num_modes = num_modes +inp.driver.typeof = "intrinsic" +inp.driver.sol_path = pathlib.Path( + "./results/DiscreteMC1verysmall") + +inp.simulation.typeof = "single" +inp.system.name = "s1" +inp.system.solution = "static" +inp.system.solver_library = "diffrax" +inp.system.solver_function = "newton" +inp.system.solver_settings = dict(rtol=1e-6, + atol=1e-6, + max_steps=50, + norm="linalg_norm", + kappa=0.01) +inp.system.xloads.follower_forces = True +inp.system.xloads.x = [0, 1, 2, 3, 4, 5] +inp.system.t = [1, 2, 3, 4, 5] +# rwing: 14-35 +# lwing: 40-61 +points = [] +interpolation = [] +_interpolationverysmall = [i*1e-3 for i in _interpolation] +_interpolationverysmall_torsion = [i*1e-3 for i in _interpolation_torsion] + +for ri,li in zip(range(14, 36),range(40,62)): + points.append([ri, 2]) + points.append([ri, 4]) + points.append([li, 2]) + points.append([li, 4]) +for i, _ in enumerate(range(len(points))): + + if i % 2 == 0: + interpolation.append(_interpolationverysmall) + else: + interpolation.append(_interpolationverysmall_torsion) + +interpolation = np.array(interpolation) # num_pointforces x num_interpolationverysmall +paths = 8*10 +sigma0 = 0.15 # percentage of mu for sigma +mu = _interpolationverysmall[-1] +sigma = (sigma0) * _interpolationverysmall[-1] +rand = np.random.normal(mu, sigma, paths) +mu_torsion = _interpolationverysmall_torsion[-1] +sigma_torsion = (sigma0) * _interpolationverysmall_torsion[-1] +rand_torsion = np.random.normal(mu_torsion, sigma_torsion, paths) +follower_interpolation = [] +for i, ri in enumerate(rand): + interpolationrand = np.copy(interpolation) + interpolationrand[::2, -1] = ri + interpolationrand[1::2, -1] = rand_torsion[i] + follower_interpolation.append(interpolationrand) +#follower_interpolation = [interpolation * ri for ri in rand] +follower_points = [points]*paths +inputforces = dict(follower_points=follower_points, + follower_interpolation=follower_interpolation + ) +inp.system.shard = dict(input_type="pointforces", + inputs=inputforces) + +sol3 = feniax.feniax_shardmain.main(input_dict=inp, device_count=8) +# np.mean(sol.staticsystem_sys1.ra[:,-1,2,35]) +# np.std(sol.staticsystem_sys1.ra[:,-1,2,35]) +# Very Small loading:1 ends here + +# [[file:modelgeneration.org::*Statistics][Statistics:1]] +u_mean = np.mean(sol1.staticsystem_sys1.ra[:,-1,2,35] - config.fem.X[35,2]) +u_std = np.std(sol1.staticsystem_sys1.ra[:,-1,2,35]) + +print(f"Mean displacement node 35: {u_mean}") +print(f"std displacement node 35: {u_std}") +print(f"Ratio displacement node 35: {u_mean/u_std}") +print("***************") + +u_mean2 = np.mean(sol2.staticsystem_sys1.ra[:,-1,2,35] - config.fem.X[35,2]) +u_std2 = np.std(sol2.staticsystem_sys1.ra[:,-1,2,35]) + +print(f"Mean displacement node 35: {u_mean2}") +print(f"std displacement node 35: {u_std2}") +print(f"ratio displacement node 35: {u_mean2/u_std2}") +print("***************") + +u_mean3 = np.mean(sol3.staticsystem_sys1.ra[:,-1,2,35] - config.fem.X[35,2]) +u_std3 = np.std(sol3.staticsystem_sys1.ra[:,-1,2,35]) + +print(f"Mean displacement node 35: {u_mean3}") +print(f"std displacement node 35: {u_std3}") +print(f"ratio displacement node 35: {u_mean3/u_std3}") +print("***************") +# Statistics:1 ends here diff --git a/examples/BUG/settings_gust1.py b/examples/BUG/settings_gust1.py index 0acfa96..690e2fc 100644 --- a/examples/BUG/settings_gust1.py +++ b/examples/BUG/settings_gust1.py @@ -1,16 +1,23 @@ # [[file:modelgeneration.org::*Run][Run:1]] import pathlib +import time import jax.numpy as jnp +import numpy as np import feniax.preprocessor.configuration as configuration # import Config, dump_to_yaml from feniax.preprocessor.inputs import Inputs import feniax.feniax_main - -label_gaf = "Dd1c7F1Scao-50" +import feniax.plotools.reconstruction as reconstruction +label_dlm = "d1c7" +sol = "eao" +label_gaf = "Dd1c7F3Seao-100" +num_modes = 100 +c_ref = 3.0 +u_inf = 209.62786434059765 +rho_inf = 0.41275511341689247 num_poles = 5 Dhj_file = f"D{label_gaf}p{num_poles}" Ahh_file = f"A{label_gaf}p{num_poles}" Poles_file = f"Poles{label_gaf}p{num_poles}" - inp = Inputs() inp.engine = "intrinsicmodal" inp.fem.eig_type = "inputs" @@ -29,42 +36,43 @@ RHTP=None, LHTP=None, ) - -inp.fem.grid = "structuralGridclamped" -inp.fem.folder = pathlib.Path('./FEM/') -inp.fem.eig_names = ["eigenvals_50.npy", "eigenvecs_50.npy"] -inp.fem.num_modes = 50 +inp.fem.grid = f"./FEM/structuralGrid_{sol[:-1]}" +#inp.fem.folder = pathlib.Path('./FEM/') +inp.fem.Ka_name = f"./FEM/Ka_{sol[:-1]}.npy" +inp.fem.Ma_name = f"./FEM/Ma_{sol[:-1]}.npy" +inp.fem.eig_names = [f"./FEM/eigenvals_{sol}{num_modes}.npy", + f"./FEM/eigenvecs_{sol}{num_modes}.npy"] +inp.driver.typeof = "intrinsic" +inp.fem.num_modes = num_modes inp.driver.typeof = "intrinsic" - -#inp.driver.sol_path = pathlib.Path( -# f"./results_{datetime.datetime.now().strftime('%Y-%m-%d_%H:%M:%S')}") -inp.driver.sol_path = pathlib.Path( - "./results1gust") inp.simulation.typeof = "single" inp.system.name = "s1" +if sol[0] == "e": # free model, otherwise clamped + inp.system.bc1 = 'free' + inp.system.q0treatment = 1 inp.system.solution = "dynamic" -inp.system.t1 = 7.5 -inp.system.tn = 4001 +inp.system.t1 = 2 +inp.system.tn = 2501 inp.system.solver_library = "runge_kutta" inp.system.solver_function = "ode" inp.system.solver_settings = dict(solver_name="rk4") inp.system.xloads.modalaero_forces = True -inp.system.q0treatment = 2 -inp.system.aero.c_ref = 3.0 -inp.system.aero.u_inf = 180. -inp.system.aero.rho_inf = 1. +inp.system.aero.c_ref = c_ref +inp.system.aero.u_inf = u_inf +inp.system.aero.rho_inf = rho_inf inp.system.aero.poles = f"./AERO/{Poles_file}.npy" inp.system.aero.A = f"./AERO/{Ahh_file}.npy" inp.system.aero.D = f"./AERO/{Dhj_file}.npy" inp.system.aero.gust_profile = "mc" -inp.system.aero.gust.intensity = 0.01 # 14.0732311562*0.001 -inp.system.aero.gust.length = 30. -inp.system.aero.gust.step = 0.05 +inp.system.aero.gust.intensity = 20 +inp.system.aero.gust.length = 150. +inp.system.aero.gust.step = 0.1 inp.system.aero.gust.shift = 0. -inp.system.aero.gust.panels_dihedral = jnp.load("./AERO/Dihedral_d1c7.npy") -inp.system.aero.gust.collocation_points = "./AERO/Collocation_d1c7.npy" +inp.system.aero.gust.panels_dihedral = f"./AERO/Dihedral_{label_dlm}.npy" +inp.system.aero.gust.collocation_points = f"./AERO/Collocation_{label_dlm}.npy" +inp.driver.sol_path = pathlib.Path( + f"./results/gust2_{sol}") config = configuration.Config(inp) sol1 = feniax.feniax_main.main(input_obj=config) - # Run:1 ends here diff --git a/examples/BUG/settings_gust1shard.py b/examples/BUG/settings_gust1shard.py index f60212d..9f17e96 100644 --- a/examples/BUG/settings_gust1shard.py +++ b/examples/BUG/settings_gust1shard.py @@ -1,15 +1,21 @@ -# [[file:modelgeneration.org::*Run][Run:1]] +# [[file:modelgeneration.org::*Gust][Gust:2]] import pathlib +import time #import jax.numpy as jnp +import numpy as np from feniax.preprocessor.inputs import Inputs import feniax.feniax_shardmain - -label_gaf = "Dd1c7F1Scao-50" +label_dlm = "d1c7" +sol = "eao" +label_gaf = "Dd1c7F3Seao-100" +num_modes = 100 +c_ref = 3.0 +u_inf = 209.62786434059765 +rho_inf = 0.41275511341689247 num_poles = 5 Dhj_file = f"D{label_gaf}p{num_poles}" Ahh_file = f"A{label_gaf}p{num_poles}" Poles_file = f"Poles{label_gaf}p{num_poles}" - inp = Inputs() inp.engine = "intrinsicmodal" inp.fem.eig_type = "inputs" @@ -28,87 +34,52 @@ RHTP=None, LHTP=None, ) - -inp.fem.grid = "structuralGridclamped" -inp.fem.folder = pathlib.Path('./FEM/') -inp.fem.eig_names = ["eigenvals_50.npy", "eigenvecs_50.npy"] -inp.fem.num_modes = 50 +inp.fem.grid = f"./FEM/structuralGrid_{sol[:-1]}" +#inp.fem.folder = pathlib.Path('./FEM/') +inp.fem.Ka_name = f"./FEM/Ka_{sol[:-1]}.npy" +inp.fem.Ma_name = f"./FEM/Ma_{sol[:-1]}.npy" +inp.fem.eig_names = [f"./FEM/eigenvals_{sol}{num_modes}.npy", + f"./FEM/eigenvecs_{sol}{num_modes}.npy"] +inp.driver.typeof = "intrinsic" +inp.fem.num_modes = num_modes inp.driver.typeof = "intrinsic" - -#inp.driver.sol_path = pathlib.Path( -# f"./results_{datetime.datetime.now().strftime('%Y-%m-%d_%H:%M:%S')}") -inp.driver.sol_path = pathlib.Path( - "./results_shardgust1") inp.simulation.typeof = "single" inp.system.name = "s1" +if sol[0] == "e": # free model, otherwise clamped + inp.system.bc1 = 'free' + inp.system.q0treatment = 1 inp.system.solution = "dynamic" -inp.system.t1 = 7.5 -inp.system.tn = 4001 +inp.system.t1 = 2 +inp.system.tn = 2501 inp.system.solver_library = "runge_kutta" inp.system.solver_function = "ode" inp.system.solver_settings = dict(solver_name="rk4") inp.system.xloads.modalaero_forces = True -inp.system.q0treatment = 2 -inp.system.aero.c_ref = 3.0 -inp.system.aero.u_inf = 180. -inp.system.aero.rho_inf = 1. +inp.system.aero.c_ref = c_ref +inp.system.aero.u_inf = u_inf +inp.system.aero.rho_inf = rho_inf inp.system.aero.poles = f"./AERO/{Poles_file}.npy" inp.system.aero.A = f"./AERO/{Ahh_file}.npy" inp.system.aero.D = f"./AERO/{Dhj_file}.npy" inp.system.aero.gust_profile = "mc" -#inp.system.aero.gust.intensity = 14.0732311562*0.001 -inp.system.aero.gust.fixed_discretisation = [50, 180] -#inp.system.aero.gust.length = 67. -inp.system.aero.gust.step = 1 +inp.system.aero.gust.intensity = 20 +inp.system.aero.gust.length = 150. +inp.system.aero.gust.step = 0.1 inp.system.aero.gust.shift = 0. -inp.system.aero.gust.panels_dihedral = "./AERO/Dihedral_d1c7.npy" -inp.system.aero.gust.collocation_points = "./AERO/Collocation_d1c7.npy" +inp.system.aero.gust.panels_dihedral = f"./AERO/Dihedral_{label_dlm}.npy" +inp.system.aero.gust.collocation_points = f"./AERO/Collocation_{label_dlm}.npy" - -inputflow = dict(length=[30,67], - intensity= [0.01,0.08, 0.2, 0.7], - #u_inf = [inp.system.aero.u_inf], - # rho_inf = [inp.system.aero.rho_inf] - ) +inp.driver.sol_path = pathlib.Path( + f"./results/gust2_{sol}Shard") +inp.system.aero.gust.fixed_discretisation = [150, u_inf] +# Shard inputs +inputflow = dict(length=np.linspace(25,265,13), + intensity=np.linspace(0.1, 3, 11), + rho_inf = np.linspace(0.34,0.48,8) + ) inp.system.shard = dict(input_type="gust1", - inputs=inputflow) - -sol = feniax.feniax_shardmain.main(input_dict=inp, device_count=8) - -# import jax -# import jax.numpy as jnp -# from functools import partial -# @partial(jax.jit, static_argnames=["a","d"]) -# def r(a,d): -# amin = min(a) -# amax = max(a) -# out = jnp.arange(amin,amax,d) -# return out - - -# import jax -# import jax.numpy as jnp -# from jax import lax - -# def compute_range_params(start, stop, step): -# # Compute range parameters outside the jitted function -# num_elements = (stop - start + step - 1) // step -# return num_elements - -# @jax.jit -# def dynamic_range_fn(start, step, num_elements): -# # Convert num_elements to a static integer -# num_elements = lax.convert_element_type(num_elements, jnp.int32) - -# # Use lax.iota to create a range -# return start + lax.iota(jnp.int32, num_elements) * step - -# # Example usage with precomputed num_elements -# start = 0 -# stop = 10 -# step = 2 -# num_elements = compute_range_params(start, stop, step) -# result = dynamic_range_fn(start, step, num_elements) -# print(result) - + inputs=inputflow) +num_gpus = 8 +solgust21shard = feniax.feniax_shardmain.main(input_dict=inp, device_count=num_gpus) +# Gust:2 ends here diff --git a/examples/BUG/settings_manoeuvre.py b/examples/BUG/settings_manoeuvre1.py similarity index 51% rename from examples/BUG/settings_manoeuvre.py rename to examples/BUG/settings_manoeuvre1.py index d9dadec..e2a62ba 100644 --- a/examples/BUG/settings_manoeuvre.py +++ b/examples/BUG/settings_manoeuvre1.py @@ -1,17 +1,23 @@ # [[file:modelgeneration.org::*Run][Run:1]] import pathlib +import time import jax.numpy as jnp +import numpy as np import feniax.preprocessor.configuration as configuration # import Config, dump_to_yaml from feniax.preprocessor.inputs import Inputs import feniax.feniax_main import feniax.plotools.reconstruction as reconstruction - -label_gaf = "Dd1c7F1Scao-50" +label_dlm = "d1c7" +sol = "eao" +label_gaf = "Dd1c7F3Seao-100" +num_modes = 100 +c_ref = 3.0 +u_inf = 209.62786434059765 +rho_inf = 0.41275511341689247 num_poles = 5 Dhj_file = f"D{label_gaf}p{num_poles}" Ahh_file = f"A{label_gaf}p{num_poles}" Poles_file = f"Poles{label_gaf}p{num_poles}" - inp = Inputs() inp.engine = "intrinsicmodal" inp.fem.eig_type = "inputs" @@ -30,54 +36,56 @@ RHTP=None, LHTP=None, ) - -inp.fem.folder = pathlib.Path('./FEM/') -inp.fem.eig_names = ["eigenvals_50.npy", "eigenvecs_50.npy"] -inp.fem.num_modes = 50 +inp.fem.grid = f"./FEM/structuralGrid_{sol[:-1]}" +#inp.fem.folder = pathlib.Path('./FEM/') +inp.fem.Ka_name = f"./FEM/Ka_{sol[:-1]}.npy" +inp.fem.Ma_name = f"./FEM/Ma_{sol[:-1]}.npy" +inp.fem.eig_names = [f"./FEM/eigenvals_{sol}{num_modes}.npy", + f"./FEM/eigenvecs_{sol}{num_modes}.npy"] +inp.driver.typeof = "intrinsic" +inp.fem.num_modes = num_modes inp.driver.typeof = "intrinsic" - -#inp.driver.sol_path = pathlib.Path( -# f"./results_{datetime.datetime.now().strftime('%Y-%m-%d_%H:%M:%S')}") -inp.driver.sol_path = pathlib.Path( - "./results1manoeuvre") inp.simulation.typeof = "single" -inp.systems.sett.s1.solution = "static" -inp.systems.sett.s1.solver_library = "diffrax" -inp.systems.sett.s1.solver_function = "newton" -inp.systems.sett.s1.solver_settings = dict(rtol=1e-6, +inp.system.name = "s1" +inp.system.solution = "static" +inp.system.solver_library = "diffrax" +inp.system.solver_function = "newton" +inp.system.solver_settings = dict(rtol=1e-6, atol=1e-6, max_steps=100, norm="linalg_norm", kappa=0.01) -inp.systems.sett.s1.xloads.modalaero_forces = True -inp.systems.sett.s1.xloads.x = [0.,1.] -inp.systems.sett.s1.t = [0.25, 0.5, 0.75, 1] -inp.systems.sett.s1.aero.c_ref = 3.0 -inp.systems.sett.s1.aero.u_inf = 170. -inp.systems.sett.s1.aero.rho_inf = 0.778 -inp.systems.sett.s1.aero.poles = f"./AERO/{Poles_file}.npy" -inp.systems.sett.s1.aero.A = f"./AERO/{Ahh_file}.npy" -inp.systems.sett.s1.aero.Q0_rigid = f"./AERO/Qhx{label_gaf}.npy" -inp.systems.sett.s1.aero.qalpha = jnp.array(([0., 0., 0, 0, 0, 0], - [0., 6 * jnp.pi / 180, 0, 0, 0, 0])) - - +inp.system.xloads.modalaero_forces = True +inp.system.xloads.x = [0.,1.] +inp.system.t = [1/3, 2/3, 1]#[0.25, 0.5, 0.75, 1] +inp.system.aero.c_ref = c_ref +inp.system.aero.u_inf = u_inf # a0(7000) =312 +inp.system.aero.rho_inf = rho_inf +inp.system.aero.poles = f"./AERO/{Poles_file}.npy" +inp.system.aero.A = f"./AERO/{Ahh_file}.npy" +inp.system.aero.Q0_rigid = f"./AERO/Qhx{label_gaf}.npy" +inp.system.aero.qalpha = [[0., 0., 0, 0, 0, 0], + [0., 6 * np.pi / 180, 0, 0, 0, 0]] # interpolation: x=0 => qalpha=0 + # x=1 => qalpha = 4 +inp.driver.sol_path = pathlib.Path( + "./results/manoeuvre2") config = configuration.Config(inp) solstatic1 = feniax.feniax_main.main(input_obj=config) # Run:1 ends here -plot3D= False -if plot3D: - reconstruction.rbf_based( + +# [[file:modelgeneration.org::3Dstatic][3Dstatic]] +rintrinsic, uintrinsic = reconstruction.rbf_based( nastran_bdf="./NASTRAN/BUG103.bdf", X=config.fem.X, - time=range(len(inp.systems.sett.s1.t)), - ra=solstatic1.staticsystem_s1.ra, - Rab=solstatic1.staticsystem_s1.Cab, + time=range(len(inp.system.t)), + ra=solstatic1.staticsystem_sys1.ra, + Rab=solstatic1.staticsystem_sys1.Cab, R0ab=solstatic1.modes.C0ab, - vtkpath="./paraview/solstatic1", + vtkpath=inp.driver.sol_path / "paraview/solstatic1/bug", plot_timeinterval=1, - plot_ref=True, + plot_ref=False, tolerance=1e-3, size_cards=8, rbe3s_full=False, ra_movie=None) +# 3Dstatic ends here diff --git a/examples/BUG/settings_manoeuvre1shard.py b/examples/BUG/settings_manoeuvre1shard.py new file mode 100644 index 0000000..fd47c12 --- /dev/null +++ b/examples/BUG/settings_manoeuvre1shard.py @@ -0,0 +1,78 @@ +# [[file:modelgeneration.org::*Manoeuvre][Manoeuvre:1]] +import pathlib +import time +#import jax.numpy as jnp +import numpy as np +from feniax.preprocessor.inputs import Inputs +import feniax.feniax_shardmain +label_dlm = "d1c7" +sol = "eao" +label_gaf = "Dd1c7F3Seao-100" +num_modes = 100 +c_ref = 3.0 +u_inf = 209.62786434059765 +rho_inf = 0.41275511341689247 +num_poles = 5 +Dhj_file = f"D{label_gaf}p{num_poles}" +Ahh_file = f"A{label_gaf}p{num_poles}" +Poles_file = f"Poles{label_gaf}p{num_poles}" +inp = Inputs() +inp.engine = "intrinsicmodal" +inp.fem.eig_type = "inputs" + +inp.fem.connectivity = dict(# FusWing=['RWing', + # 'LWing'], + FusBack=['FusTail', + 'VTP'], + FusFront=None, + RWing=None, + LWing=None, + FusTail=None, + VTP=['HTP', 'VTPTail'], + HTP=['RHTP', 'LHTP'], + VTPTail=None, + RHTP=None, + LHTP=None, + ) +inp.fem.grid = f"./FEM/structuralGrid_{sol[:-1]}" +#inp.fem.folder = pathlib.Path('./FEM/') +inp.fem.Ka_name = f"./FEM/Ka_{sol[:-1]}.npy" +inp.fem.Ma_name = f"./FEM/Ma_{sol[:-1]}.npy" +inp.fem.eig_names = [f"./FEM/eigenvals_{sol}{num_modes}.npy", + f"./FEM/eigenvecs_{sol}{num_modes}.npy"] +inp.driver.typeof = "intrinsic" +inp.fem.num_modes = num_modes +inp.driver.typeof = "intrinsic" +inp.simulation.typeof = "single" +inp.system.name = "s1" +inp.system.solution = "static" +inp.system.solver_library = "diffrax" +inp.system.solver_function = "newton" +inp.system.solver_settings = dict(rtol=1e-6, + atol=1e-6, + max_steps=100, + norm="linalg_norm", + kappa=0.01) +inp.system.xloads.modalaero_forces = True +inp.system.xloads.x = [0.,1.] +inp.system.t = [1/3, 2/3, 1]#[0.25, 0.5, 0.75, 1] +inp.system.aero.c_ref = c_ref +inp.system.aero.u_inf = u_inf # a0(7000) =312 +inp.system.aero.rho_inf = rho_inf +inp.system.aero.poles = f"./AERO/{Poles_file}.npy" +inp.system.aero.A = f"./AERO/{Ahh_file}.npy" +inp.system.aero.Q0_rigid = f"./AERO/Qhx{label_gaf}.npy" +inp.system.aero.qalpha = [[0., 0., 0, 0, 0, 0], + [0., 6 * np.pi / 180, 0, 0, 0, 0]] # interpolation: x=0 => qalpha=0 + # x=1 => qalpha = 4 + +#rho_rand = np.random.normal(0.6, 0.6*0.15, 500) +inp.driver.sol_path = pathlib.Path( + "./results/manoeuvre1Shard") +inputflow = dict(u_inf=np.linspace(190, 240, 8), + rho_inf=np.linspace(0.41, 0.81, 8)) +inp.system.shard = dict(input_type="steadyalpha", + inputs=inputflow) +num_gpus = 8 +solstatic1shard = feniax.feniax_shardmain.main(input_dict=inp, device_count=num_gpus) +# Manoeuvre:1 ends here diff --git a/examples/SailPlane/settingsShard.py b/examples/SailPlane/settingsShard.py index c2ab2cb..3703ef6 100644 --- a/examples/SailPlane/settingsShard.py +++ b/examples/SailPlane/settingsShard.py @@ -55,12 +55,12 @@ [[20, 5], [43, 5]], [[15, 5], [38, 5]], [[10, 5], [33, 5]], - ], + ], # [[[node_i, component_j]..(total_forces per run)],...(parallel forces)[[node_i, component_j]..]] follower_interpolation= [[[0.,2e5,2.5e5,3.e5,4.e5,4.8e5,5.3e5], [0.,2e5,2.5e5,3.e5,4.e5,4.8e5,5.3e5] ] - ]*8 - ) + ]*8 # [[[0,...interpolation points]..(total_forces per run)],...(parallel forces)[[0,...]..]] + ) inp.system.shard = dict(input_type="pointforces", inputs=inputforces) #config = configuration.Config(inp) diff --git a/feniax/intrinsic/args.py b/feniax/intrinsic/args.py index ef4b15e..c0328a8 100644 --- a/feniax/intrinsic/args.py +++ b/feniax/intrinsic/args.py @@ -690,6 +690,7 @@ def arg_20g546( *args, **kwargs, ): + eta_0 = kwargs["eta_0"] gamma1 = sol.data.couplings.gamma1 gamma2 = sol.data.couplings.gamma2 @@ -711,20 +712,20 @@ def arg_20g546( omega, phi1l, states, + aero.poles, num_modes, num_poles, + gust.x, + c_ref, aero.A0hat, aero.A1hat, aero.A2hatinv, aero.A3hat, u_inf, - c_ref, - aero.poles, - gust.x, F1g, Flg, ) - + @catter2library def arg_20G546( diff --git a/feniax/intrinsic/argshard.py b/feniax/intrinsic/argshard.py index cbb2560..62f548b 100644 --- a/feniax/intrinsic/argshard.py +++ b/feniax/intrinsic/argshard.py @@ -90,3 +90,46 @@ def arg_20g21( ) ) + +def arg_20g546( + sol: solution.IntrinsicSolution, + system: intrinsicmodal.Dsystem, + fem: intrinsicmodal.Dfem, + *args, + **kwargs, +): + + eta_0 = kwargs["eta_0"] + phi1l = sol.data.modes.phi1l + phi2l = sol.data.modes.phi2l + psi2l = sol.data.modes.psi2l + X_xdelta = sol.data.modes.X_xdelta + omega = sol.data.modes.omega + X_xdelta = sol.data.modes.X_xdelta + C0ab = sol.data.modes.C0ab + gamma1 = sol.data.couplings.gamma1 + gamma2 = sol.data.couplings.gamma2 + states = system.states + c_ref = system.aero.c_ref + num_modes = fem.num_modes + num_poles = system.aero.num_poles + poles = system.aero.poles + A = system.aero.A + D = system.aero.D + xgust = system.aero.gust.time + collocation_x = system.aero.gust.collocation_points[:,0] + return (phi1l, phi2l, psi2l, X_xdelta, C0ab, A, D, c_ref, + ( + eta_0, + gamma1, + gamma2, + omega, + phi1l, + states, + poles, + num_modes, + num_poles, + xgust, + ) + + ) diff --git a/feniax/intrinsic/dq_dynamic.py b/feniax/intrinsic/dq_dynamic.py index 847013c..e9a4143 100644 --- a/feniax/intrinsic/dq_dynamic.py +++ b/feniax/intrinsic/dq_dynamic.py @@ -352,23 +352,43 @@ def dq_20g273(t, q, *args): def dq_20g546(t, q, *args): """Gust response free flight, q0 obtained via integrator q1.""" + # ( + # eta_0, + # gamma1, + # gamma2, + # omega, + # phi1l, + # states, + # num_modes, + # num_poles, + # A0hat, + # A1hat, + # A2hatinv, + # A3hat, + # u_inf, + # c_ref, + # poles, + # xgust, + # F1gust, + # Flgust, + # ) = args[0] ( eta_0, gamma1, gamma2, omega, - phi1l, + phi1l, states, + poles, num_modes, num_poles, + xgust, + c_ref, A0hat, A1hat, A2hatinv, A3hat, u_inf, - c_ref, - poles, - xgust, F1gust, Flgust, ) = args[0] diff --git a/feniax/intrinsic/dynamicShard.py b/feniax/intrinsic/dynamicShard.py index 5594bd0..714eb9b 100644 --- a/feniax/intrinsic/dynamicShard.py +++ b/feniax/intrinsic/dynamicShard.py @@ -151,3 +151,110 @@ def _main_20g21_3(inp): main_vmap = jax.vmap(_main_20g21_3) results = main_vmap(inputs) return results + +@partial(jax.jit, static_argnames=["config"]) +def main_20g546_3( + inputs, # + q0, + config, + args, + **kwargs, +): + + X = config.fem.X + phi1l, phi2l, psi2l, X_xdelta, C0ab, A, D, c_ref, _dqargs = args + states = _dqargs[5] + q1_index = states["q1"] + q2_index = states["q2"] + # q1_index = states["q1"] + # q2_index = states["q2"] + + collocation_points = config.system.aero.gust.collocation_points + gust_shift = config.system.aero.gust.shift + dihedral = config.system.aero.gust.panels_dihedral + fshape_span = igust._get_spanshape(config.system.aero.gust.shape) + # gust_totaltime = config.system.aero.gust.totaltime + time_gust = config.system.aero.gust.time + # @jax.jit + def _main_20g546_3(inp): + + rho_inf = inp[0] + u_inf = inp[1] + gust_length = inp[2] + gust_intensity = inp[3] + q_inf = 0.5 * rho_inf * u_inf**2 + gust_totaltime = gust_length / u_inf + + A0hat = q_inf * A[0] + A1hat = c_ref * rho_inf * u_inf / 4 * A[1] + A2hat = c_ref**2 * rho_inf / 8 * A[2] + A3hat = q_inf * A[3:] + A2hatinv = jnp.linalg.inv(jnp.eye(len(A2hat)) - A2hat) + D0hat = q_inf * D[0] + D1hat = c_ref * rho_inf * u_inf / 4 * D[1] + D2hat = c_ref**2 * rho_inf / 8 * D[2] + D3hat = q_inf * D[3:] + + gust, gust_dot, gust_ddot = igust._downwashRogerMc( + u_inf, + gust_length, + gust_intensity, + gust_shift, + collocation_points, + dihedral, # normals, + time_gust, + gust_totaltime, + fshape_span, + ) + Q_w, Q_wdot, Q_wddot, Q_wsum, Ql_wdot = igust._getGAFs( + D0hat, # NbxNm + D1hat, + D2hat, + D3hat, + gust, + gust_dot, + gust_ddot, + ) + + args_inp = (c_ref, + A0hat, + A1hat, + A2hatinv, + A3hat, + u_inf, + Q_wsum, + Ql_wdot + ) + + dq_args = _dqargs + args_inp + states_puller, eqsolver = sollibs.factory( + config.system.solver_library, config.system.solver_function + ) + + sol = eqsolver( + dq_dynamic.dq_20g546, + dq_args, + config.system.solver_settings, + q0=q0, + t0=config.system.t0, + t1=config.system.t1, + tn=config.system.tn, + dt=config.system.dt, + t=config.system.t, + ) + # jax.debug.breakpoint() + q = states_puller(sol) + q1 = q[:, q1_index] + q2 = q[:, q2_index] + tn = len(q) + # X2, X3, ra, Cab = isys.recover_staticfields(q2, tn, X, + # phi2l, psi2l, X_xdelta, C0ab, config.fem) + X1, X2, X3, ra, Cab = isys.recover_fieldsRB( + q1, q2, tn, config.system.dt, X, phi1l, phi2l, psi2l, X_xdelta, C0ab, config + ) + + return dict(q=q, X1=X1,X2=X2, X3=X3, ra=ra, Cab=Cab) + + main_vmap = jax.vmap(_main_20g546_3) + results = main_vmap(inputs) + return results diff --git a/feniax/intrinsic/postprocess.py b/feniax/intrinsic/postprocess.py index e266a9d..4e5f88e 100644 --- a/feniax/intrinsic/postprocess.py +++ b/feniax/intrinsic/postprocess.py @@ -45,7 +45,7 @@ def velocity_ra(): ... def strains_ra(): ... - +@jax.jit def integrate_node0(X1, dt, ra_n0, Rab_n0): v_average = (X1[:-1, :3] + X1[1:, :3]) / 2 theta_average = (X1[:-1, 3:6] + X1[1:, 3:6]) / 2 * dt diff --git a/feniax/intrinsic/xloads.py b/feniax/intrinsic/xloads.py index 167c225..f5ae3d9 100644 --- a/feniax/intrinsic/xloads.py +++ b/feniax/intrinsic/xloads.py @@ -135,6 +135,26 @@ def eta_pointdead(t, phi1, x, force_dead, Rab): eta = jnp.tensordot(phi1, f_fd, axes=([1, 2], [0, 1])) return eta +@jax.jit +def force_pointdead(t, x, force_dead): + f = linear_interpolation(t, x, force_dead) + return f + + +@jax.jit +def force_pointfollower(t, x, force_follower, Rab): + f1 = jax.vmap( + lambda R, x: jnp.vstack( + [jnp.hstack([R, jnp.zeros((3, 3))]), jnp.hstack([jnp.zeros((3, 3)), R])] + ) + @ x, + in_axes=(2, 1), + out_axes=1, + ) + f = linear_interpolation(t, x, force_follower) + f_fd = f1(Rab, f) + return f_fd + @jax.jit def eta_pointdead_const(phi1, f, Rab): diff --git a/feniax/plotools/streamlit/intrinsic.py b/feniax/plotools/streamlit/intrinsic.py index 2a0c8e3..813c5cc 100644 --- a/feniax/plotools/streamlit/intrinsic.py +++ b/feniax/plotools/streamlit/intrinsic.py @@ -86,8 +86,8 @@ def show_vtu(vtu_folder, stream_out="streamlit"): def df_geometry(fem): st.header("Geometry and FE Models") st.subheader("Condensed model") - st.table(fem.df_grid) - + with st.expander("See dataframe with geometry nodes"): + st.table(fem.df_grid) def fe_matrices(fem): st.subheader("Ka and Ma input matrices") @@ -134,7 +134,8 @@ def df_modes(sol, config): df = pd.DataFrame( mvalue[mnumber].T * scale, columns=["x", "y", "z", "rx", "ry", "rz"] ) - st.table(df) + with st.expander("See Modes data-frame"): + st.table(df) st.divider() fig = modes_3Dconfiguration( mode=mvalue[mnumber] * scale, config=config, mode_label=mname, settings=None @@ -323,7 +324,7 @@ def sys_X_comparison(X, t, labels, ylabel, mode="lines"): componenti = None col1, col2 = st.columns(2) ntimes, ncomponents, nnodes = X[0].shape - + dashstyle = [None, 'dash', 'dot', 'dashdot']*10 nodei = col1.selectbox("Select a node", options=range(nnodes)) componenti = col2.selectbox("Select a component", options=range(ncomponents)) # breakpoint() @@ -335,7 +336,7 @@ def sys_X_comparison(X, t, labels, ylabel, mode="lines"): fig, dict( name=labels[i], - # line=dict(color="navy"), + line=dict(dash=dashstyle[i]), mode=mode, ), dict(title="Time evolution"), @@ -355,15 +356,34 @@ def sys_displacements_comp(solsys, config): t = [list(range(len(qi))) for qi in Xvs[0].q] mode = "lines+markers" - ra = [x.ra - config[labels[i]].fem.X.T for i, x in enumerate(Xvs)] - sys_X_comparison(ra, t, labels, "Displacements", mode) + try: + ua = [x.ra - config[labels[i]].fem.X.T for i, x in enumerate(Xvs)] + except KeyError: + label0 = list(config.keys())[0] + ua = [x.ra - config[label0].fem.X.T for i, x in enumerate(Xvs)] + sys_X_comparison(ua, t, labels, "Displacements", mode) +def sys_positions_comp(solsys): + labels = list(solsys.keys()) + Xvs = list(solsys.values()) + if hasattr(Xvs[0], "t"): + t = [x.t for x in Xvs] + mode = "lines" + else: + t = [list(range(len(qi))) for qi in Xvs[0].q] + mode = "lines+markers" + ra = [x.ra for i, x in enumerate(Xvs)] + sys_X_comparison(ra, t, labels, "Positions", mode) + def sys_displacements(solsys, config): sys_X(solsys.ra - config.fem.X.T, solsys, label="Displacements") - +def sys_positions(solsys): + sys_X(solsys.ra, solsys, label="Positions") + def sys_velocities_comp(solsys): + Xvs = list(solsys.values()) if hasattr(Xvs[0], "X1"): labels = list(solsys.keys()) @@ -373,7 +393,6 @@ def sys_velocities_comp(solsys): else: st.text("Static solution!! All velocities are 0") - def sys_velocities(solsys): if hasattr(solsys, "X1"): sys_X(solsys.X1, solsys, label="X1") @@ -604,6 +623,7 @@ def systems(sol, config): [ "STATES", "DISPLACEMENTS", + "POSITIONS", "VELOCITIES", "STRAINS", "INTERNALFORCES", @@ -630,6 +650,8 @@ def systems(sol, config): sys_states(solsys) case show.DISPLACEMENTS.name: sys_displacements(solsys, config) + case show.POSITIONS.name: + sys_positions(solsys) case show.VELOCITIES.name: sys_velocities(solsys) case show.STRAINS.name: @@ -653,6 +675,7 @@ def systems_comparison(sol, config): [ "STATES", "DISPLACEMENTS", + "POSITIONS", "VELOCITIES", "STRAINS", "INTERNALFORCES", @@ -684,6 +707,8 @@ def systems_comparison(sol, config): sys_states_comparison(solsys) case show.DISPLACEMENTS.name: sys_displacements_comp(solsys, config) + case show.POSITIONS.name: + sys_positions(solsys) case show.VELOCITIES.name: sys_velocities_comp(solsys) case show.STRAINS.name: @@ -733,6 +758,7 @@ def systems_shard(sol_shard, config_shard, sol=None, config=None): [ "STATES", "DISPLACEMENTS", + "POSITIONS", "VELOCITIES", "STRAINS", "INTERNALFORCES", @@ -792,6 +818,8 @@ def systems_shard(sol_shard, config_shard, sol=None, config=None): sys_states_comparison(selected_solsys) case show.DISPLACEMENTS.name: sys_displacements_comp(selected_solsys, config_shard) + case show.POSITIONS.name: + sys_positions_comp(selected_solsys) case show.VELOCITIES.name: sys_velocities_comp(selected_solsys) case show.STRAINS.name: diff --git a/feniax/plotools/streamlit/pages/Geometry.py b/feniax/plotools/streamlit/pages/Geometry.py index 52d6a3a..1cddc8a 100644 --- a/feniax/plotools/streamlit/pages/Geometry.py +++ b/feniax/plotools/streamlit/pages/Geometry.py @@ -2,6 +2,9 @@ import streamlit as st import feniax.plotools.streamlit.theory as stt +import importlib +importlib.reload(sti) + st.set_page_config( page_title="Initial model geometry", page_icon="🛫", @@ -19,3 +22,4 @@ if left.button("Click to see FE matrices", use_container_width=True): sti.fe_matrices(st.session_state.config.fem) + diff --git a/feniax/plotools/upyvista.py b/feniax/plotools/upyvista.py index bf256d0..9d141f3 100644 --- a/feniax/plotools/upyvista.py +++ b/feniax/plotools/upyvista.py @@ -18,9 +18,7 @@ def render_wireframe(points, lines, pl: pv.Plotter = None): pl = pv.Plotter() mesh = pv.PolyData(points, lines) - pl.add_mesh( - mesh, color="black", line_width=5, style="wireframe", render_lines_as_tubes=True - ) + pl.add_mesh(mesh, color="black", line_width=5, style="wireframe", render_lines_as_tubes=True) pl.add_points( points, @@ -43,3 +41,32 @@ def render_mesh(points, lines): mesh = pv.PolyData(points, lines) return mesh + +# import pyvista as pv + +# points = np.random.rand(100, 3) +# mesh = pv.PolyData(points) +# mesh.plot(point_size=10, style='points', color='tan') + +# polydata = pv.PolyData(points) + +# # Add the vectors as point data +# polydata["vectors"] = vectors + +# # Create glyphs to represent vectors +# glyphs = polydata.glyph(orient="vectors", scale=False, factor=0.3) + +# # Plot the vector field +# plotter = pv.Plotter() +# plotter.add_mesh(glyphs, color='red') +# plotter.add_mesh(polydata, color='blue', point_size=5, render_points_as_spheres=True) +# plotter.show() + +# mesh = pyvista.PolyData(v, cells) +# mesh.save(folder_path / f"collocation_{k}.ply", binary=False) + + +# X=config.fem.X, +# time=range(len(inp.system.t)), +# ra=sol.staticsystem_sys1.ra[i], +# Rab=sol.staticsystem_sys1.Cab[i], diff --git a/feniax/preprocessor/containers/intrinsicmodal.py b/feniax/preprocessor/containers/intrinsicmodal.py index bfbd74a..5fff03a 100644 --- a/feniax/preprocessor/containers/intrinsicmodal.py +++ b/feniax/preprocessor/containers/intrinsicmodal.py @@ -723,6 +723,8 @@ def __post_init__(self): object.__setattr__(self, "C", jnp.load(self.C)) if isinstance(self.D, (str, pathlib.Path)): object.__setattr__(self, "D", jnp.load(self.D)) + if self.qalpha is not None and not isinstance(self.qalpha, jnp.ndarray): + object.__setattr__(self, "qalpha", jnp.array(self.qalpha)) self._initialize_attributes() diff --git a/feniax/preprocessor/containers/intrinsicsol.py b/feniax/preprocessor/containers/intrinsicsol.py index ea8643a..c8cd478 100644 --- a/feniax/preprocessor/containers/intrinsicsol.py +++ b/feniax/preprocessor/containers/intrinsicsol.py @@ -105,7 +105,8 @@ class GustRoger: @dataclass(slots=True) class Shards: points: jnp.ndarray = None - + device_count: int = None + local_device_count: int = None # import dataclasses diff --git a/feniax/preprocessor/solution.py b/feniax/preprocessor/solution.py index 15cc608..85520c3 100644 --- a/feniax/preprocessor/solution.py +++ b/feniax/preprocessor/solution.py @@ -79,6 +79,8 @@ def save_container(path, container): attr_path = path / attr_name if isinstance(attr, jnp.ndarray): jnp.save(attr_path, attr) + elif isinstance(attr, int) or isinstance(attr, float): + jnp.save(attr_path, jnp.array(attr)) elif isinstance(attr, np.ndarray): jnp.save(attr_path, attr) elif isinstance(attr, (list, dict, tuple)): @@ -102,7 +104,12 @@ def load_container(path: pathlib.Path, Container): if Container.__annotations__[attr_name].__name__ == "Array": kwargs[attr_name] = jnp.load(attr_path.with_suffix(".npy")) elif Container.__annotations__[attr_name].__name__ == "ndarray": - kwargs[attr_name] = np.load(attr_path.with_suffix(".npy")) + kwargs[attr_name] = jnp.load(attr_path.with_suffix(".npy")) + elif (Container.__annotations__[attr_name].__name__ == "int"): + kwargs[attr_name] = int(jnp.load(attr_path.with_suffix(".npy"))) + elif (Container.__annotations__[attr_name].__name__ == "float"): + kwargs[attr_name] = float(jnp.load(attr_path.with_suffix(".npy"))) + elif ( (Container.__annotations__[attr_name].__name__ == "dict") or (Container.__annotations__[attr_name].__name__ == "list") diff --git a/feniax/systems/intrinsicShard.py b/feniax/systems/intrinsicShard.py index 61bde6a..873bb6a 100644 --- a/feniax/systems/intrinsicShard.py +++ b/feniax/systems/intrinsicShard.py @@ -83,6 +83,8 @@ def build_solution(self): "Shards", label="_" + self.name, points=self.xpoints, + device_count=jax.device_count(), + local_device_count=jax.local_device_count() ) if self.settings.save: self.sol.save_container("Shards", label="_" + self.name) diff --git a/feniax/systems/intrinsic_system.py b/feniax/systems/intrinsic_system.py index e7c824a..26d5d8c 100644 --- a/feniax/systems/intrinsic_system.py +++ b/feniax/systems/intrinsic_system.py @@ -38,6 +38,27 @@ def recover_fields(q1, q2, tn, X, phi1l, phi2l, psi2l, X_xdelta, C0ab, config): return X1, X2, X3, ra, Cab +@partial(jax.jit, static_argnames=["config", "tn", "dt"]) +def recover_fieldsRB(q1, q2, tn, dt, X, phi1l, phi2l, psi2l, X_xdelta, C0ab, config): + ra_n0 = X[0] + Rab_n0 = jnp.eye(3) + X1 = postprocess.compute_velocities(phi1l, q1) + + X2 = postprocess.compute_internalforces(phi2l, q2) + X3 = postprocess.compute_strains(psi2l, q2) + + Cab0, ra0 = postprocess.integrate_node0( + X1[:, :, 0], dt, ra_n0, Rab_n0 + ) + Cab, ra = postprocess.integrate_strains_t( + ra0, + Cab0, + X3, + X_xdelta, + C0ab, + config, + ) + return X1, X2, X3, ra, Cab @partial(jax.jit, static_argnames=["config", "tn"]) def recover_staticfields(q2, tn, X, phi2l, psi2l, X_xdelta, C0ab, config): diff --git a/feniax/utils.py b/feniax/utils.py index 069307f..bcac2c9 100644 --- a/feniax/utils.py +++ b/feniax/utils.py @@ -7,7 +7,7 @@ from pyNastran.f06 import parse_flutter as flutter import matplotlib.pyplot as plt import pandas as pd - +import numpy as np def flatten_list(lis): l = list(lis) @@ -29,3 +29,23 @@ def flatten_list(lis): i = i + 1 return l + +def standard_atmosphere(h, k=0.0065, R=287.05, g=9.806, gamma=1.4, T_0=288.16, rho_0=1.225): + n = 1 / (1 - k * R / g) + if h < 11000.0: + T = T_0 - k * h + rho = rho_0 * (T / T_0) ** (1 / (n - 1)) + P = rho * R * T + a = np.sqrt(gamma * R * T) + elif 11000.0 <= h <= 25000.0: + h_11k = 11000.0 + T_11k = T_0 - k * h_11k + rho_11k = rho_0 * (T_11k / T_0) ** (1 / (n - 1)) + P_11k = rho_11k * R * T_11k + + psi = np.exp(-(h - h_11k) * g / (R * T_11k)) + T = T_11k + rho = rho_11k * psi + P = P_11k * psi + a = np.sqrt(gamma * R * T_11k) + return T, rho, P, a diff --git a/pyproject.toml b/pyproject.toml index f0bf879..f8738e9 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -37,7 +37,7 @@ dependencies = ["numpy", "pandas", #"PyYAML", "ruamel.yaml", - "jax<=0.4.31", + "jax", "jaxlib", "diffrax", "jaxopt",