Skip to content
Snippets Groups Projects

use_custom

Last edited by 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. 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

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

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

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

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

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

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

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

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

adj = sdpm.Adjoint(sol)
adj.solve()

Post-processing the solution

The solution can be accessed in Python, mainly by using

# 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

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 derivatives of the Generalized Aerodynamic Force (GAF) matrices can be calculated and accessed using

grd = sdpm.Gradient(sol)
grd.run()
grd.computePartialsGafRe(i_freq, i_mode, i_row, j_col, wrt)
grd.computePartialsGafIm(i_freq, i_mode, i_row, j_col, wrt)

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).
If automatic differentiation has been enabled, the sensitivities of the solution can be computed using

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

[dx0] = [df0_dx0 df1_dx0] [df0]  
[dx1]   [df0_dx1 df1_dx1] [df1]