Skip to content

Changes

Page history
Update use_custom authored by Adrien Crovato's avatar Adrien Crovato
Hide whitespace changes
Inline Side-by-side
use_custom.md
View page @ b396595a
## Custom scripts ## Custom scripts
Users can customize their use of SDPM by writing their own scripts. Each script follows a similar pattern: Users can customize their use of SDPM by writing their own scripts. Each script follows a similar pattern:
1. Creating the mesh 1. Creating the mesh
2. Defining the problem 2. Defining the problem
3. Solving the flow 3. Solving the flow
4. Post-processing the results 4. Post-processing the results
_Note_ SDPM can be used in a development or in a production environment. For the former, the code needs to be compiled, and then launched using [run.py](https://gitlab.uliege.be/am-dept/sdpm/blob/master/run.py). For the latter, the code needs to be installed to a location which is in the python path (the default install path is the user's python site-packages). _Note_ SDPM can be used in a development or in a production environment. For the former, the code needs to be compiled, and then launched using [run.py](https://gitlab.uliege.be/am-dept/sdpm/blob/master/run.py). For the latter, the code needs to be installed to a location which is in the python path (the default install path is the user's python site-packages).
### Creating the mesh ### Creating the mesh
The mesh is created using The mesh is created using
```python ```python
msh = sdpm.gmsh.MeshLoader(model_path).run(pars) msh = sdpm.gmsh.MeshLoader(model_path).run(pars)
``` ```
where `model_path` is the path to the file containing the model and `pars` is an optional dictionary of parameters to configure Gmsh. The model can either be a geometry (.geo) or a mesh (.msh). In the latter case, any parameters are ignored. where `model_path` is the path to the file containing the model and `pars` is an optional dictionary of parameters to configure Gmsh. The model can either be a geometry (.geo) or a mesh (.msh). In the latter case, any parameters are ignored.
### Defining the problem ### Defining the problem
The problem is defined as follows The problem is defined as follows
```python ```python
pbl = sdpm.Problem(msh) pbl = sdpm.Problem(msh)
pbl.setGeometry(sref, cref, xref, yref, zref, ysym) pbl.setGeometry(sref, cref, xref, yref, zref, ysym)
pbl.setFreestream(aoa, aos, minf) pbl.setFreestream(aoa, aos, minf)
pbl.setUnsteady(n_modes, kred) pbl.setUnsteady(n_modes, kred)
``` ```
where `sref` is the reference surface area, `cref` is the reference chord length, `xref`, `yref` and `zref` are the coordinates of the reference center around which aerodynamic moments are computed, `ysym` indicates whether only a half model is provided and y-symmetry boundary condition should be applied, `aoa` and `aos` are the angles of attack and of side slip, `minf` is the freestream Mach number, `n_modes` is the number of unsteady motions and `kred` is an array of reduced frequencies. Note that the last line can be ignored for steady problems. where `sref` is the reference surface area, `cref` is the reference chord length, `xref`, `yref` and `zref` are the coordinates of the reference center around which aerodynamic moments are computed, `ysym` indicates whether only a half model is provided and y-symmetry boundary condition should be applied, `aoa` and `aos` are the angles of attack and of side slip, `minf` is the freestream Mach number, `n_modes` is the number of unsteady motions and `kred` is an array of reduced frequencies. Note that the last line can be ignored for steady problems.
Lifting bodies are then added using Lifting bodies are then added using
```python ```python
wing = sdpm.Body(msh, body_name, te_name, cref, n_div, n_modes, n_kred) wing = sdpm.Body(msh, body_name, te_name, cref, n_div, n_modes, n_kred)
pbl.addBody(wing) pbl.addBody(wing)
``` ```
where `body_name` and `te_name` are the names of the lifting body to be added and of its trailing edge in the Gmsh model, `n_div` is the number of chordwise panels, and `n_kred` is the number of reduced frequencies. For steady problems, the last two arguments default to 0 and `n_div` can be set to 1. where `body_name` and `te_name` are the names of the lifting body to be added and of its trailing edge in the Gmsh model, `n_div` is the number of chordwise panels, and `n_kred` is the number of reduced frequencies. For steady problems, the last two arguments default to 0 and `n_div` can be set to 1.
Unsteady rigid body motions can be added using Unsteady rigid body motions can be added using
```python ```python
wing.addMotion() wing.addMotion()
wing.setRigidMotion(i_mode, r_amp, d_amp, x_ref, z_ref) wing.setRigidMotion(i_mode, r_amp, d_amp, x_ref, z_ref)
``` ```
where `i_mode` is the number of the motion being configured, `r_amp` is the amplitude of rotation about y, `d_amp` is the amplitude of displacement along z, and `x_ref` and `z_ref` are the x and z coordinates of the points around which the body rotates. where `i_mode` is the number of the motion being configured, `r_amp` is the amplitude of rotation about y, `d_amp` is the amplitude of displacement along z, and `x_ref` and `z_ref` are the x and z coordinates of the points around which the body rotates.
Unsteady flexible motions (mode shapes) can be added using Unsteady flexible motions (mode shapes) can be added using
```python ```python
x_modes, amp_modes = sdpm.utils.interpolate_modes(modes_path, wing.getNodes()) x_modes, amp_modes = sdpm.utils.interpolate_modes(modes_path, wing.getNodes())
for i in range(n_modes): for i in range(n_modes):
wing.addMotion() wing.addMotion()
wing.setFlexibleMotion(i, 1 / amp_modes[i], x_modes[:, 3*i:3*i+3]) wing.setFlexibleMotion(i, 1 / amp_modes[i], x_modes[:, 3*i:3*i+3])
``` ```
where `modes_path` is the path to the file containing the mode shapes. Each row of this file must contain the x, y and z coordinates of each point on the body surface, as well as the components corresponding to the displacement along z and to the rotations around x and y. The data must then be organized as: x, y, z, dz0, rx0, ry0, dz1, rx1, ry1, etc. The utility `interpolate_modes` uses SciPy RBF implementation to interpolate the modes on the aerodynamic grid. where `modes_path` is the path to the file containing the mode shapes. Each row of this file must contain the x, y and z coordinates of each point on the body surface, as well as the components corresponding to the displacement along z and to the rotations around x and y. The data must then be organized as: x, y, z, dz0, rx0, ry0, dz1, rx1, ry1, etc. The utility `interpolate_modes` uses SciPy RBF implementation to interpolate the modes on the aerodynamic grid.
A transonic correction can be added using A transonic correction can be added using
```python ```python
dcp = sdpm.utils.interpolate_pressure(pres_path, wing.getNodes()) dcp = sdpm.utils.interpolate_pressure(pres_path, wing.getNodes())
wing.setTransonicPressureGrad(dcp[:,0]) wing.setTransonicPressureGrad(dcp[:,0])
``` ```
where `pres_path` is the path to the file containing the reference pressure derivative with respect to the angle of attack. Each row of this file must contain the x, y and z coordinates of each point on the body surface, as well as its associated pressure derivative. The utility `interpolate_pressure` uses SciPy RBF implementation to interpolate the pressure derivative on the aerodynamic grid. where `pres_path` is the path to the file containing the reference pressure derivative with respect to the angle of attack. Each row of this file must contain the x, y and z coordinates of each point on the body surface, as well as its associated pressure derivative. The utility `interpolate_pressure` uses SciPy RBF implementation to interpolate the pressure derivative on the aerodynamic grid.
### Solving the flow ### Solving the flow
Once the problem has been fully defined, the solver can be instantiated, updated and run using Once the problem has been fully defined, the solver can be instantiated, updated and run using
```python ```python
sol = sdpm.Solver(pbl) sol = sdpm.Solver(pbl)
sol.update() sol.update()
sol.run() sol.run()
``` ```
If a transonic correction is to be applied, `SolverTransonic` can be used instead of `Solver`. The solution can be saved in Gmsh format using If a transonic correction is to be applied, `SolverTransonic` can be used instead of `Solver`. The solution can be saved in Gmsh format using
```python ```python
sol.save(sdpm.GmshExport(msh)) sol.save(sdpm.GmshExport(msh))
``` ```
If the gradients of the solution are needed, automatic differentiation must be enabled and the solver must be run using If the gradients of the solution are needed, automatic differentiation must be enabled and the solver must be run using
```python ```python
adj = sdpm.Adjoint(sol) adj = sdpm.Adjoint(sol)
adj.solve() adj.solve()
``` ```
### Post-processing the solution ### Post-processing the solution
The solution can be accessed in Python, mainly by using The solution can be accessed in Python, mainly by using
```python ```python
# Steady and unsteady pressure # Steady and unsteady pressure
sol.getPressure() sol.getPressure()
sol.getUnsteadyPressureReal(i_mode, i_freq) sol.getUnsteadyPressureReal(i_mode, i_freq)
sol.getUnsteadyPressureImag(i_mode, i_freq) sol.getUnsteadyPressureImag(i_mode, i_freq)
# Steady and unsteady aerodynamic load coefficients # Steady and unsteady aerodynamic load coefficients
sol.getLiftCoeff() sol.getLiftCoeff()
sol.getDragCoeff() sol.getDragCoeff()
sol.getSideforceCoeff() sol.getSideforceCoeff()
sol.getMomentCoeff() sol.getMomentCoeff()
sol.getUnsteadyLiftCoeffReal(i_mode, i_freq) sol.getUnsteadyLiftCoeffReal(i_mode, i_freq)
sol.getUnsteadyLiftCoeffImag(i_mode, i_freq) sol.getUnsteadyLiftCoeffImag(i_mode, i_freq)
sol.getUnsteadyDragCoeffReal(i_mode, i_freq) sol.getUnsteadyDragCoeffReal(i_mode, i_freq)
sol.getUnsteadyDragCoeffImag(i_mode, i_freq) sol.getUnsteadyDragCoeffImag(i_mode, i_freq)
sol.getUnsteadySideforceCoeffReal(i_mode, i_freq) sol.getUnsteadySideforceCoeffReal(i_mode, i_freq)
sol.getUnsteadySideforceCoeffImag(i_mode, i_freq) sol.getUnsteadySideforceCoeffImag(i_mode, i_freq)
sol.getUnsteadyMomentCoeffReal(i_mode, i_freq) sol.getUnsteadyMomentCoeffReal(i_mode, i_freq)
sol.getUnsteadyMomentCoeffImag(i_mode, i_freq) sol.getUnsteadyMomentCoeffImag(i_mode, i_freq)
# Generalized aerodynamic forces matrix # Generalized aerodynamic forces matrix
sol.getGafMatrixRe(i_freq) sol.getGafMatrixRe(i_freq)
sol.getGafMatrixIm(i_freq) sol.getGafMatrixIm(i_freq)
``` ```
where `i_mode`, `i_freq` and the indices of the requested mode and frequency. where `i_mode`, `i_freq` and the indices of the requested mode and frequency.
Additionally, chordwise pressure distributions can be extracted along the wing span by using Additionally, chordwise pressure distributions can be extracted along the wing span by using
```python ```python
sdpm.utils.create_slices(wing, vars, model_name, y_sec) sdpm.utils.create_slices(wing, vars, model_name, y_sec)
``` ```
where `vars` is a dictionary containing the name and its corresponding data array (e.g. `{'cp': sol.getPressure()}`), `model_name` is the name of the model (inherited from the model file) and `y_sec` is an array of y coordinates. where `vars` is a dictionary containing the name and its corresponding data array (e.g. `{'cp': sol.getPressure()}`), `model_name` is the name of the model (inherited from the model file) and `y_sec` is an array of y coordinates.
The sensitivities of the solution can be computed using The derivatives of the Generalized Aerodynamic Force (GAF) matrices can be calculated and accessed using
```python ```python
sens = adj.compute(of, seeds, wrt) grd = sdpm.Gradient(sol)
``` grd.run()
where `of` is the name of the variable for which the sensitivities are requested, `seeds` is an array of seed values and `wrt` is an array of names of the variables which will be differentiated with respect to. The sensitivities will be stored in the dictionary `sens` and can be accessed using `sens[wrt_var][i]` where `wrt_var` is one of the variable of `wrt` and `i` is the index of interest if the variable `wrt_var` is an array. grd.computePartialsGafRe(i_freq, i_mode, i_row, j_col, wrt)
grd.computePartialsGafIm(i_freq, i_mode, i_row, j_col, wrt)
_Note_ currently, the variables `cl`, `cd`, `cl1_re`, `cl1_im`, `cd1_re` and `cd1_im` can be differentiated with respect to `aoa` and `x_wing`. ```
_Note_ here, the sensitivity is the (adjoint) transposed matrix-vector product between the gradient and the seed. If `df` and `dx` are the seed and the sensitivity vectors, it can be expressed as The last two lines will return the derivative of the real and imaginary parts of the GAF matrix entry at row `i_row` and column `j_col` with respect to the component `wrt` of mode number `i_mode` at frequency `i_freq` for all body nodes. `wrt` can either be `'a'` (rigid pitch), `'h'` (rigid plunge), `'dz'` (modal displacements along z), `'rx'` (modal rotations about x) or `'ry'` (modal rotations about y).
```python If automatic differentiation has been enabled, the sensitivities of the solution can be computed using
[dx0] = [df0_dx0 df1_dx0] [df0] ```python
[dx1] [df0_dx1 df1_dx1] [df1] sens = adj.compute(of, seeds, wrt)
```
where `of` is the name of the variable for which the sensitivities are requested, `seeds` is an array of seed values and `wrt` is an array of names of the variables which will be differentiated with respect to. The sensitivities will be stored in the dictionary `sens` and can be accessed using `sens[wrt_var][i]` where `wrt_var` is one of the variable of `wrt` and `i` is the index of interest if the variable `wrt_var` is an array.
_Note_ currently, the variables `cl`, `cd`, `cl1_re`, `cl1_im`, `cd1_re` and `cd1_im` can be differentiated with respect to `aoa` and `x_wing`.
_Note_ here, the sensitivity is the (adjoint) transposed matrix-vector product between the gradient and the seed. If `df` and `dx` are the seed and the sensitivity vectors, it can be expressed as
```python
[dx0] = [df0_dx0 df1_dx0] [df0]
[dx1] [df0_dx1 df1_dx1] [df1]
``` ```
\ No newline at end of file