From be59c21185c0fcd80d53f9443b207012458d87fe Mon Sep 17 00:00:00 2001
From: acrovato <a.crovato@uliege.be>
Date: Fri, 4 Oct 2024 11:51:38 +0200
Subject: [PATCH] Add partials to API. Fix bug in row assignement for
multibody.
---
sdpm/api/core.py | 14 ++++++-------
sdpm/api/om.py | 36 ++++++++++++++++++++------------
sdpm/src/sdpmGradient.cpp | 3 ++-
sdpm/src/sdpmSolverTransonic.cpp | 3 ++-
4 files changed, 34 insertions(+), 22 deletions(-)
diff --git a/sdpm/api/core.py b/sdpm/api/core.py
index 10a355e..de386ad 100644
--- a/sdpm/api/core.py
+++ b/sdpm/api/core.py
@@ -17,15 +17,13 @@
## Initialize SDPM
# Adrien Crovato
-def init_sdpm(cfg, use_ad=False):
+def init_sdpm(cfg):
"""Initialize SDPM components
Parameters
----------
cfg : dict
Dictionary of parameters to configure SDPM
- use_ad : bool, optional
- Whether to use AD within SDPM (default: False)
Returns
-------
@@ -44,8 +42,8 @@ def init_sdpm(cfg, use_ad=False):
Lifting body
sol : sdpm.Solver object
Direct solver
- adj : sdpm.Adjoint object
- Adjoint solver
+ grd : sdpm.Gradient object
+ Gradient calculator
"""
# Imports
import sdpm
@@ -87,7 +85,9 @@ def init_sdpm(cfg, use_ad=False):
# Solver
vrb = cfg.get('Verb_lvl', 1)
_sol = sdpm.Solver(_pbl, vrb) if 'Transonic_pressure_grad' not in cfg else sdpm.SolverTransonic(_pbl, vrb)
- _adj = sdpm.Adjoint(_sol) if use_ad else None
+
+ # Gradient
+ _grd = sdpm.Gradient(_sol)
# Return
return {
@@ -98,5 +98,5 @@ def init_sdpm(cfg, use_ad=False):
'pbl': _pbl,
'bdy': _bdy,
'sol': _sol,
- 'adj': _adj
+ 'grd': _grd
}
diff --git a/sdpm/api/om.py b/sdpm/api/om.py
index 9fe3a5c..d5edd72 100644
--- a/sdpm/api/om.py
+++ b/sdpm/api/om.py
@@ -50,20 +50,20 @@ class SdpmSolver(om.ExplicitComponent):
self.options.declare('nm', desc='number of modes')
self.options.declare('bdy', desc='lifting body', recordable=False)
self.options.declare('sol', desc='direct solver', recordable=False)
- self.options.declare('adj', desc='adjoint solver', recordable=False)
+ self.options.declare('grd', desc='gradient calculator', recordable=False)
def setup(self):
self.nf = self.options['nf']
self.nm = self.options['nm']
self.bdy = self.options['bdy']
self.sol = self.options['sol']
- self.adj = self.options['adj']
+ self.grd = self.options['grd']
self.add_input('x_aero0', shape_by_conn=True, desc='aerodynamic surface node coordinates')
self.add_input('q_aero', shape_by_conn=True, desc='modal displacements')
self.add_output('Q_re', val=np.ones((self.nf, self.nm, self.nm)), desc='real part of generalized aerodynamic forces matrices')
self.add_output('Q_im', val=np.ones((self.nf, self.nm, self.nm)), desc='imag part of generalized aerodynamic forces matrices')
# Partials
- # self.declare_partials(of=['Q_re', 'Q_im'], wrt=['x_aero', 'q_aero'], method='exact')
+ self.declare_partials(of=['Q_re', 'Q_im'], wrt=['q_aero'], method='exact')
def compute(self, inputs, outputs):
# Update coordinates
@@ -75,16 +75,28 @@ class SdpmSolver(om.ExplicitComponent):
x_modes[j, :] = inputs['q_aero'][3*j:3*j+3, i].T
self.bdy.setFlexibleMotion(i, 1, x_modes)
# Compute
- if self.adj is None:
- self.sol.update()
- self.sol.run()
- else:
- self.adj.solve()
+ self.sol.update()
+ self.sol.run()
# Set outputs
for i in range(self.nf):
outputs['Q_re'][i, :, :] = np.array(self.sol.getGafMatrixRe(i), dtype=float)
outputs['Q_im'][i, :, :] = np.array(self.sol.getGafMatrixIm(i), dtype=float)
+ def compute_partials(self, inputs, partials):
+ # Compute gradients and set partials
+ self.grd.run()
+ for ifrq in range(self.nf):
+ for imod in range(self.nm):
+ for jmod in range(self.nm):
+ dq_re = np.zeros((inputs['q_aero'].shape[0], self.nm))
+ dq_im = np.zeros((inputs['q_aero'].shape[0], self.nm))
+ for kmod in range(self.nm):
+ for idof, dof in enumerate(['dz', 'rx', 'ry']):
+ dq_re[idof::3, kmod] = np.array(self.grd.computePartialsGafRe(ifrq, kmod, imod, jmod, dof), dtype=float)
+ dq_im[idof::3, kmod] = np.array(self.grd.computePartialsGafIm(ifrq, kmod, imod, jmod, dof), dtype=float)
+ partials['Q_re', 'q_aero'][ifrq * self.nm * self.nm + imod * self.nm + jmod, :] = dq_re.flatten()
+ partials['Q_im', 'q_aero'][ifrq * self.nm * self.nm + imod * self.nm + jmod, :] = dq_im.flatten()
+
class SdpmBuilder(Builder):
"""SDPM builder for OpenMDAO
@@ -99,18 +111,16 @@ class SdpmBuilder(Builder):
'sol': direct solver
'adj': adjoint sovler
"""
- def __init__(self, cfg, use_ad=False):
+ def __init__(self, cfg):
"""Instantiate and initialize SDPM components
Parameters
----------
cfg : dict
Dictionary of parameters to configure SDPM
- use_ad : bool, optional
- Whether to use AD within SDPM (default: False)
"""
from .core import init_sdpm
- self.__sdpm = init_sdpm(cfg, use_ad)
+ self.__sdpm = init_sdpm(cfg)
def get_mesh(self):
"""Return OpenMDAO component to get the initial surface mesh coordinates
@@ -120,7 +130,7 @@ class SdpmBuilder(Builder):
def get_solver(self, scenario_name=''):
"""Return OpenMDAO component containing the solver
"""
- return SdpmSolver(nf=self.__sdpm['n_f'], nm=self.__sdpm['n_m'], bdy=self.__sdpm['bdy'], sol=self.__sdpm['sol'], adj=self.__sdpm['adj'])
+ return SdpmSolver(nf=self.__sdpm['n_f'], nm=self.__sdpm['n_m'], bdy=self.__sdpm['bdy'], sol=self.__sdpm['sol'], grd=self.__sdpm['grd'])
def get_number_of_nodes(self):
"""Return the number of surface nodes
diff --git a/sdpm/src/sdpmGradient.cpp b/sdpm/src/sdpmGradient.cpp
index 1603da9..5af4c58 100644
--- a/sdpm/src/sdpmGradient.cpp
+++ b/sdpm/src/sdpmGradient.cpp
@@ -30,7 +30,8 @@ Gradient::Gradient(Solver &sol) : _sol(sol),
int inod = 0;
for (auto body : _sol._pbl.getBodies())
for (auto n : body->getNodes())
- _rowsN.insert(std::pair<Node *, int>(n, inod++));
+ if (_rowsN.insert(std::pair<Node *, int>(n, inod)).second == true)
+ ++inod;
_nN = _rowsN.size();
// Resize displacement and load containers
diff --git a/sdpm/src/sdpmSolverTransonic.cpp b/sdpm/src/sdpmSolverTransonic.cpp
index 8783512..6033212 100644
--- a/sdpm/src/sdpmSolverTransonic.cpp
+++ b/sdpm/src/sdpmSolverTransonic.cpp
@@ -28,7 +28,8 @@ SolverTransonic::SolverTransonic(Problem &pbl, int vrb) : Solver(pbl, vrb), _dCo
int i = 0;
for (auto body : _pbl.getBodies())
for (auto n : body->getNodes())
- _rowsN.insert(std::pair<Node *, int>(n, i++));
+ if (_rowsN.insert(std::pair<Node *, int>(n, i)).second == true)
+ ++i;
}
/**
--
GitLab