Create use_custom authored by Adrien Crovato's avatar Adrien Crovato
## Custom scripts
Users can customize their use of SDPM by writing their own scripts. Each script follows a similar pattern:
1. Creating the mesh
2. Defining the problem
3. Solving the flow
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).
### Creating the mesh
The mesh is created using
```python
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.
### Defining the problem
The problem is defined as follows
```python
pbl = sdpm.Problem(msh)
pbl.setGeometry(sref, cref, xref, yref, zref, ysym)
pbl.setFreestream(aoa, aos, minf)
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.
Lifting bodies are then added using
```python
wing = sdpm.Body(msh, body_name, te_name, cref, n_div, n_modes, n_kred)
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.
Unsteady rigid body motions can be added using
```python
wing.addMotion()
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.
Unsteady flexible motions (mode shapes) can be added using
```python
x_modes, amp_modes = sdpm.utils.interpolate_modes(modes_path, wing.getNodes())
for i in range(n_modes):
wing.addMotion()
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.
A transonic correction can be added using
```python
dcp = sdpm.utils.interpolate_pressure(pres_path, wing.getNodes())
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.
### Solving the flow
Once the problem has been fully defined, the solver can be instantiated, updated and run using
```python
sol = sdpm.Solver(pbl)
sol.update()
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
```python
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
```python
adj = sdpm.Adjoint(sol)
adj.solve()
```
### Post-processing the solution
The solution can be accessed in Python, mainly by using
```python
# Steady and unsteady pressure
sol.getPressure()
sol.getUnsteadyPressureReal(i_mode, i_freq)
sol.getUnsteadyPressureImag(i_mode, i_freq)
# Steady and unsteady aerodynamic load coefficients
sol.getLiftCoeff()
sol.getDragCoeff()
sol.getSideforceCoeff()
sol.getMomentCoeff()
sol.getUnsteadyLiftCoeffReal(i_mode, i_freq)
sol.getUnsteadyLiftCoeffImag(i_mode, i_freq)
sol.getUnsteadyDragCoeffReal(i_mode, i_freq)
sol.getUnsteadyDragCoeffImag(i_mode, i_freq)
sol.getUnsteadySideforceCoeffReal(i_mode, i_freq)
sol.getUnsteadySideforceCoeffImag(i_mode, i_freq)
sol.getUnsteadyMomentCoeffReal(i_mode, i_freq)
sol.getUnsteadyMomentCoeffImag(i_mode, i_freq)
# Generalized aerodynamic forces matrix
sol.getGafMatrixRe(i_freq)
sol.getGafMatrixIm(i_freq)
```
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
```python
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.
The sensitivities of the solution can be computed using
```python
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