Skip to content
Snippets Groups Projects
Verified Commit 31238be0 authored by Paul Dechamps's avatar Paul Dechamps :speech_balloon:
Browse files

(feat) Refactor optimzation problem setup

parent 86f5fd63
No related branches found
No related tags found
No related merge requests found
Pipeline #51054 passed
Hide whitespace changes
Inline Side-by-side
blast/api/mda_api.py
+ 107
− 2
View file @ 31238be0
......@@ -24,7 +24,6 @@ from fwk.coloring import ccolors
import openmdao.api as om
import blast
import blast.blUtils as vutils
from pygeo import DVGeometryCST
from scipy.optimize import least_squares
import fwk
......@@ -79,7 +78,7 @@ class BlasterSolver(om.ExplicitComponent):
self.coupler.run(write=False)
self.__write(sfx='baseline')
self.coupler.reset()
print('object created')
print('BLASTER objects created')
def setup(self):
self.add_input('aoa', val=0.0, desc='Flow angle of attack')
......@@ -205,6 +204,7 @@ class BlasterSolver(om.ExplicitComponent):
def getBoundaryFFDGeometry(self):
"""Airfoil FFD geometry
"""
from pygeo import DVGeometryCST
with open('surface_DVGeo.dat', 'w') as f:
for i, node_coord in enumerate(self.isol.vBnd[0][0].nodesCoord):
f.write(f"{node_coord[0]} {node_coord[1]}\n")
......@@ -498,3 +498,108 @@ class BlasterFFDGeometry(om.ExplicitComponent):
Ri = calc_R(*center)
R = Ri.mean()
return center, R
class groupMDA(om.Group):
def __init__(self, icfg, vcfg, ocfg):
super().__init__()
self.icfg = icfg
self.vcfg = vcfg
self.ocfg = ocfg
if 'aoa' not in self.ocfg['initialValues']:
raise ValueError('Initial value for aoa is missing')
if self.ocfg['geometryType'] == 'NACA':
if 'tau' not in self.ocfg['initialValues']:
raise ValueError('Initial value for tau is missing')
if 'eps' not in self.ocfg['initialValues']:
raise ValueError('Initial value for eps is missing')
if 'p' not in self.ocfg['initialValues']:
raise ValueError('Initial value for p is missing')
for key in self.ocfg['initialValues']:
if key not in ['aoa', 'tau', 'eps', 'p']:
raise KeyError('Invalid key in initialValues', key)
elif self.ocfg['geometryType'] == 'FFD':
if 'upper_shape' not in self.ocfg['initialValues']:
raise ValueError('Initial value for upper_shape is missing')
if 'lower_shape' not in self.ocfg['initialValues']:
raise ValueError('Initial value for lower_shape is missing')
if 'class_shape_n1' not in self.ocfg['initialValues']:
raise ValueError('Initial value for class_shape_n1 is missing')
if 'class_shape_n2' not in self.ocfg['initialValues']:
raise ValueError('Initial value for class_shape_n2 is missing')
for key in self.ocfg['initialValues']:
if key not in ['aoa', 'upper_shape', 'lower_shape', 'class_shape_n1', 'class_shape_n2']:
raise KeyError('Invalid key in initialValues', key)
def setup(self):
# Add the MDAO components
blaster = BlasterSolver(self.vcfg, self.icfg)
print('Optimization setup')
print('------------------')
# Design variables box
self.add_subsystem('dvs', om.IndepVarComp(), promotes=['*'])
self.add_subsystem('mod', BlasterModeUpdater(tolerance=10., force_mode=self.ocfg['mode'])) # Mode updater
self.add_subsystem('mesh', blaster.getBoundaryMesh())
# Geometry component
if self.ocfg['geometryType'] == 'NACA':
print('Adding NACA geometry')
self.add_subsystem('geometry', blaster.getBoundaryNACAGeometry())
elif self.ocfg['geometryType'] == 'FFD':
print('Adding FFD geometry')
self.add_subsystem('geometry', blaster.getBoundaryFFDGeometry())
# Add an intermediate variable for the constraint (such that upper_shape[0] = lower_shape[0])
self.add_subsystem('constraint_comp', om.ExecComp('constraint_var = lower_shape[0] + upper_shape[0]',
lower_shape=np.zeros(len(self.ocfg['initialValues']['upper_shape']['value'])), upper_shape=np.zeros(len(self.ocfg['initialValues']['lower_shape']['value']))), promotes=['*'])
else:
raise ValueError('Invalid geometryType', self.ocfg['geometryType'])
# Solver
self.add_subsystem('solver', blaster)
# Objectives
for obj in self.ocfg['objectives']:
print('Adding objective', obj['name'])
self.add_subsystem('obj_cmp_'+obj['name'], obj['object'], promotes=obj['promotions'])
# Connect meshes
self.connect('mesh.x_aero0', 'geometry.x_aero0')
self.connect('geometry.x_aero_out', 'solver.x_aero')
self.connect('mod.visc_mode', 'solver.visc_mode')
if 'nonlinear_solver' in self.ocfg:
self.nonlinear_solver = self.ocfg['nonlinearSolver']
if 'linear_solver' in self.ocfg:
self.linear_solver = self.ocfg['linearSolver']
print('')
def configure(self):
print('Optimization config')
print('------------------')
for variable, setup in self.ocfg['initialValues'].items():
print(f'- dvs output {variable} linked to {setup["link"]} with value {setup["value"]}')
self.dvs.add_output(variable, val=setup['value'])
self.connect(variable, setup['link']+'.'+variable)
if variable in self.ocfg['designVariables']:
design_config = self.ocfg['designVariables'][variable]
print(f'- design variable {variable} with bounds {design_config["lower_bound"]} {design_config["upper_bound"]} and scaler {design_config["scaler"]}')
self.add_design_var(variable, lower=design_config['lower_bound'], upper=design_config['upper_bound'], scaler=design_config['scaler'])
for const in self.ocfg['constraints']:
if const['type'] == 'equality':
print('- equality constraint:', const['name'], '=', const['value'])
self.add_constraint(const['name'], equals=const['value'], scaler=const['scaler'])
elif const['type'] == 'inequality':
print('- inequality constraint:', const['name'], end=' ')
if 'lower' in const:
print('-- lower bound:', const['lower'])
self.add_constraint(const['name'], lower=const['lower'], scaler=const['scaler'])
elif 'upper' in const:
print('-- upper bound:', const['upper'])
self.add_constraint(const['name'], upper=const['upper'], scaler=const['scaler'])
if self.ocfg['geometryType'] == 'FFD':
self.add_constraint('constraint_var', equals=0.0, linear=True)
for obj in self.ocfg['objectives']:
print('- objective:', obj['name'], 'linked to', obj['link'])
self.connect(obj['link']+'.'+obj['name'], obj['name'])
self.add_objective('obj', scaler=obj['scaler'])
print('')
blast/mdao/opti_NACA2.py 0 → 100644
+ 196
− 0
View file @ 31238be0
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# Copyright 2024 University of Liège
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
# @author Paul Dechamps
# @date 2024
# Optimization case using NACA parametrization
from fwk.wutils import parseargs
import fwk
from fwk.testing import *
from fwk.coloring import ccolors
from blast.api.mda_api import groupMDA
import numpy as np
import openmdao.api as om
def cfgInviscid(nthrds, verb):
import os.path
# Parameters
return {
# Options
'scenario' : 'aerodynamic',
'task' : 'optimization',
'Threads' : nthrds, # number of threads
'Verb' : verb, # verbosity
# Model (geometry or mesh)
'File' : os.path.dirname(os.path.dirname(os.path.abspath(__file__))) + '/models/dart/n0012.geo', # Input file containing the model
'Pars' : {'xLgt' : 50, 'yLgt' : 50, 'msF' : 30, 'msTe' : 0.01, 'msLe' : 0.0075}, # parameters for input file model
'Dim' : 2, # problem dimension
'Format' : 'gmsh', # save format (vtk or gmsh)
# Markers
'Fluid' : 'field', # name of physical group containing the fluid
'Farfield' : ['upstream', 'farfield', 'downstream'], # LIST of names of physical groups containing the farfield boundaries (upstream/downstream should be first/last element)
'Wing' : 'airfoil', # LIST of names of physical groups containing the lifting surface boundary
'Wake' : 'wake', # LIST of names of physical group containing the wake
'WakeTip' : 'wakeTip', # LIST of names of physical group containing the edge of the wake
'Te' : 'te', # LIST of names of physical group containing the trailing edge
'dbc' : False,
'Upstream' : 'upstream',
# Freestream
'M_inf' : 0.1, # freestream Mach number
'AoA' : 0., # freestream angle of attack
# Geometry
'S_ref' : 1., # reference surface length
'c_ref' : 1., # reference chord length
'x_ref' : .25, # reference point for moment computation (x)
'y_ref' : 0.0, # reference point for moment computation (y)
'z_ref' : 0.0, # reference point for moment computation (z)
# Numerical
'LSolver' : 'GMRES', # inner solver (Pardiso, MUMPS or GMRES)
'G_fill' : 2, # fill-in factor for GMRES preconditioner
'G_tol' : 1e-5, # tolerance for GMRES
'G_restart' : 50, # restart for GMRES
'Rel_tol' : 1e-6, # relative tolerance on solver residual
'Abs_tol' : 1e-8, # absolute tolerance on solver residual
'Max_it' : 20 # solver maximum number of iterations
}
def cfgBlast(verb):
return {
'Re' : 1e7, # Freestream Reynolds number
'Minf' : 0.1, # Freestream Mach number (used for the computation of the time step only)
'CFL0' : 1., # Inital CFL number of the calculation
'Verb': verb, # Verbosity level of the solver
'couplIter': 50, # Maximum number of iterations
'couplTol' : 1e-3, # Tolerance of the VII methodology
'iterPrint': 5, # int, number of iterations between outputs
'resetInv' : True, # bool, flag to reset the inviscid calculation at every iteration.
'sections' : [0], # List of sections for boundary layer calculation
'xtrF' : [0.05, 0.05], # Forced transition location
'interpolator' : 'Matching',# Interpolator for the coupling
'nSections' : 1, # Number of sections in the domain
'span' : 1.0 # Wing span
}
def cfgOpti():
return{
'geometryType': 'NACA',
'mode': 'visc',
'linearSolver': om.DirectSolver(),
'objectives': [{'name': 'cd',
'link': 'solver',
'object': om.ExecComp('obj = cd'),
'promotions': ['obj', 'cd'],
'scaler': 1000}
],
'initialValues': {'aoa':{'link':'solver', 'value': 0.0},
'tau':{'link':'geometry', 'value': 0.12},
'eps':{'link':'geometry', 'value': 0.0},
'p':{'link': 'geometry', 'value': 0.4}
},
'designVariables': {'aoa':{'lower_bound': -3.0, 'upper_bound': 3.0, 'scaler': 1},
'eps':{'lower_bound': 0.0, 'upper_bound': 0.05, 'scaler': 10},
'p':{'lower_bound': 0.35, 'upper_bound': 0.45, 'scaler': 10}
},
'constraints': [{'name': 'solver.cl', 'type': 'equality', 'value': 0.4, 'scaler': 10}]
}
def main():
# Timer.
tms = fwk.Timers()
tms['total'].start()
args = parseargs()
icfg = cfgInviscid(args.k, args.verb)
vcfg = cfgBlast(args.verb)
ocfg = cfgOpti()
prob = om.Problem()
prob.model = groupMDA(icfg, vcfg, ocfg)
prob.driver = om.ScipyOptimizeDriver()
prob.driver.options['optimizer'] = 'SLSQP'
prob.driver.options['tol'] = 1e-3
recorder = om.SqliteRecorder('reports/BlasterSolver_file.sql')
prob.driver.add_recorder(recorder)
prob.driver.recording_options['record_objectives'] = True
prob.driver.recording_options['record_desvars'] = True
prob.driver.recording_options['includes'] = ['aoa', 'solver.cl', 'solver.cd', 'mesh.x_aero0', 'geometry.x_aero_out']
prob.setup()
#prob.check_partials(includes=['cl', 'aoa'], compact_print=True)
prob.run_driver()
# Instantiate your CaseReader
cr = om.CaseReader(prob.get_reports_dir() / "../BlasterSolver_file.sql")
# Isolate "problem" as your source
driver_cases = cr.list_cases('driver', out_stream=None)
# Get evolution of f and x
aoa_evol = []
cl_evol = []
cd_evol = []
shape = []
shape0 = []
for c in driver_cases:
case = cr.get_case(c)
aoa_evol.append(case.get_val('aoa')[0]*180/np.pi)
cl_evol.append(case.get_val('solver.cl')[0])
cd_evol.append(case.get_val('solver.cd')[0])
shape0.append(case.get_val('mesh.x_aero0'))
shape.append(case.get_val('geometry.x_aero_out'))
print('success')
print('aoa:', prob['aoa']*np.pi/180)
print('cl:', prob['solver.cl'])
print('cd:', prob['solver.cd'])
print('obj:', prob['obj'])
# Recover selig format
shape0_selig = np.array(shape0[0])
shape1_selig = np.array(shape[-1])
for i, val in enumerate(prob.model.solver.isol.vBnd[0][0].connectListNodes):
for j in range(prob.model.solver.isol.getnDim()):
shape0_selig[i*prob.model.solver.isol.getnDim()+j] = shape0[0][val*prob.model.solver.isol.getnDim()+j]
shape1_selig[i*prob.model.solver.isol.getnDim()+j] = shape[-1][val*prob.model.solver.isol.getnDim()+j]
from matplotlib import pyplot as plt
plt.plot(shape0_selig[0::2], shape0_selig[1::2], '--', color='grey', label='Initial shape - cl = {:.2f} - cd = {:.5f}'.format(cl_evol[0], cd_evol[0]))
plt.plot(shape1_selig[0::2], shape1_selig[1::2], '-', color='blue', label='Final shape - cl = {:.2f} - cd = {:.5f}'.format(cl_evol[-1], cd_evol[-1]))
plt.legend()
plt.axis('equal')
plt.show()
quit()
plt.plot(shape0[0][::2], shape0[0][1::2], 'x', color='grey')
for xx in shape:
plt.plot(xx[::2], xx[1::2], 'x', color='blue')
plt.axis('equal')
plt.show()
quit()
if __name__ == "__main__":
main()
blast/mdao/opti_transonic.py
+ 1
− 1
View file @ 31238be0
......@@ -63,7 +63,7 @@ def cfgInviscid(nthrds, verb):
'y_ref' : 0.0, # reference point for moment computation (y)
'z_ref' : 0.0, # reference point for moment computation (z)
# Numerical
'LSolver' : 'SparseLu', # inner solver (Pardiso, MUMPS or GMRES)
'LSolver' : 'SparseLU', # inner solver (Pardiso, MUMPS or GMRES)
'G_fill' : 2, # fill-in factor for GMRES preconditioner
'G_tol' : 1e-5, # tolerance for GMRES
'G_restart' : 50, # restart for GMRES
......
blast/mdao/opti_transonic2.py 0 → 100644
+ 179
− 0
View file @ 31238be0
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# Copyright 2024 University of Liège
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
# @author Paul Dechamps
# @date 2024
# Optimization for a transonic case
from fwk.wutils import parseargs
import fwk
from fwk.testing import *
from fwk.coloring import ccolors
from blast.api.mda_api import groupMDA
import numpy as np
import openmdao.api as om
def cfgInviscid(nthrds, verb):
import os.path
# Parameters
return {
# Options
'scenario' : 'aerodynamic',
'task' : 'optimization',
'Threads' : nthrds, # number of threads
'Verb' : verb, # verbosity
# Model (geometry or mesh)
'File' : os.path.dirname(os.path.dirname(os.path.abspath(__file__))) + '/models/dart/n0012.geo', # Input file containing the model
'Pars' : {'xLgt' : 50, 'yLgt' : 50, 'msF' : 20, 'msTe' : 0.01, 'msLe' : 0.008}, # parameters for input file model
'Dim' : 2, # problem dimension
'Format' : 'gmsh', # save format (vtk or gmsh)
# Markers
'Fluid' : 'field', # name of physical group containing the fluid
'Farfield' : ['upstream', 'farfield', 'downstream'], # LIST of names of physical groups containing the farfield boundaries (upstream/downstream should be first/last element)
'Wing' : 'airfoil', # LIST of names of physical groups containing the lifting surface boundary
'Wake' : 'wake', # LIST of names of physical group containing the wake
'WakeTip' : 'wakeTip', # LIST of names of physical group containing the edge of the wake
'Te' : 'te', # LIST of names of physical group containing the trailing edge
'Upstream' : 'upstream',
# Freestream
'M_inf' : 0.78, # freestream Mach number
'AoA' : 0., # freestream angle of attack
# Geometry
'S_ref' : 1., # reference surface length
'c_ref' : 1., # reference chord length
'x_ref' : .25, # reference point for moment computation (x)
'y_ref' : 0.0, # reference point for moment computation (y)
'z_ref' : 0.0, # reference point for moment computation (z)
# Numerical
'LSolver' : 'SparseLU', # inner solver (Pardiso, MUMPS or GMRES)
'G_fill' : 2, # fill-in factor for GMRES preconditioner
'G_tol' : 1e-5, # tolerance for GMRES
'G_restart' : 50, # restart for GMRES
'Rel_tol' : 1e-6, # relative tolerance on solver residual
'Abs_tol' : 1e-8, # absolute tolerance on solver residual
'Max_it' : 50 # solver maximum number of iterations
}
def cfgBlast(verb):
return {
'Re' : 1e8, # Freestream Reynolds number
'Minf' : 0.78, # Freestream Mach number (used for the computation of the time step only)
'CFL0' : 1., # Inital CFL number of the calculation
'Verb': verb, # Verbosity level of the solver
'couplIter': 75, # Maximum number of iterations
'couplTol' : 1e-4, # Tolerance of the VII methodology
'iterPrint': 5, # int, number of iterations between outputs
'resetInv' : True, # bool, flag to reset the inviscid calculation at every iteration.
'sections' : [0], # List of sections for boundary layer calculation
'xtrF' : [0.05, 0.05], # Forced transition location
'interpolator' : 'Matching',# Interpolator for the coupling
'nSections' : 1, # Number of sections in the domain
'span' : 1.0 # Wing span
}
def cfgOpti():
return{
'geometryType': 'FFD',
'mode': 'visc',
'linearSolver': om.DirectSolver(),
'objectives': [{'name': 'cd',
'link': 'solver',
'object': om.ExecComp('obj = cd'),
'promotions': ['obj', 'cd'],
'scaler': 1000}
],
'initialValues': {'aoa':{'link':'solver', 'value': 0.0},
'upper_shape':{'link':'geometry', 'value': [0.17553542, 0.14723497, 0.14531762, 0.13826955]},
'lower_shape':{'link':'geometry', 'value': [-0.16935506, -0.15405552, -0.13912672, -0.14168947]},
'class_shape_n1':{'link': 'geometry', 'value': 0.5},
'class_shape_n2':{'link': 'geometry', 'value': 1.0}
},
'designVariables': {'aoa':{'lower_bound': -5.0, 'upper_bound': 5.0, 'scaler': 1.0},
'upper_shape':{'lower_bound': -0.15, 'upper_bound': 0.25, 'scaler': 10},
'lower_shape':{'lower_bound': -0.2, 'upper_bound': 0.2, 'scaler': 10}
},
'constraints': [{'name': 'solver.cl', 'type': 'equality', 'value': 0.4, 'scaler': 10, 'linear': False},
{'name': 'geometry.thickness', 'type': 'inequality', 'lower': 0.12, 'scaler':10, 'upper':None, 'linear': True},
{'name': 'geometry.LEradius', 'type': 'equality', 'value': 0.0224, 'scaler':100, 'linear': True}]
}
def main():
# Timer.
tms = fwk.Timers()
tms['total'].start()
args = parseargs()
icfg = cfgInviscid(args.k, args.verb)
vcfg = cfgBlast(args.verb)
ocfg = cfgOpti()
prob = om.Problem()
prob.model = groupMDA(icfg, vcfg, ocfg)
prob.driver = om.ScipyOptimizeDriver()
prob.driver.options['optimizer'] = 'SLSQP'
prob.driver.options['tol'] = 1e-4
recorder = om.SqliteRecorder('reports/BlasterSolver_file.sql')
prob.driver.add_recorder(recorder)
prob.driver.recording_options['record_objectives'] = True
prob.driver.recording_options['record_desvars'] = True
prob.driver.recording_options['includes'] = ['aoa', 'solver.cl', 'solver.cd', 'mesh.x_aero0', 'geometry.x_aero_out']
prob.setup()
# print('check')
# prob.check_partials(includes=['cl', 'cd'], compact_print=True)
# print('done')
prob.run_driver()
# problem timers
print(prob.model.solver.tms)
# Case reader
cr = om.CaseReader(prob.get_reports_dir() / "../BlasterSolver_file.sql")
driver_cases = cr.list_cases('driver', out_stream=None)
# Get evolution of f and x
aoa_evol = []
cl_evol = []
cd_evol = []
shape = []
shape0 = []
for c in driver_cases:
case = cr.get_case(c)
aoa_evol.append(case.get_val('aoa')[0])
cl_evol.append(case.get_val('solver.cl')[0])
cd_evol.append(case.get_val('solver.cd')[0])
shape0.append(case.get_val('mesh.x_aero0'))
shape.append(case.get_val('geometry.x_aero_out'))
print('aoa:', prob['aoa'])
print('cl:', prob['solver.cl'])
print('cd:', prob['solver.cd'])
print('obj:', prob['obj'])
# eof
tms['total'].stop()
print(tms)
if __name__ == "__main__":
main()
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment