Debugging solvers – a hands-on set of problems
Contents
Debugging solvers – a hands-on set of problems¶
Level: Advanced
Topics: Model construction and debugging solving
Let’s work through some solver issues together! Each one of these setups will have some sort of problem or potential problem and I want you to find a solution or fix.
Circuit example with a Newton solver¶
Here we have an example system entirely defined in the notebook. Then we try to converge it using a Newton solver but it fails. What’s going on here? How would you tackle this? Please modify this code so that the solver converges.
Make sure to check out the solver debugging checklist and to use the tips from there.
Click here to reveal a small hint
What solver type is converging the nonlinear system?
Click here to reveal a big hint
The model is missing a certain type of solver.
import numpy as np
import openmdao.api as om
class Resistor(om.ExplicitComponent):
"""Computes current across a resistor using Ohm's law."""
def initialize(self):
self.options.declare('R', default=1., desc='Resistance in Ohms')
def setup(self):
self.add_input('V_in', units='V')
self.add_input('V_out', units='V')
self.add_output('I', units='A')
# partial derivs are constant, so we can assign their values in setup
R = self.options['R']
self.declare_partials('I', 'V_in', val=1 / R)
self.declare_partials('I', 'V_out', val=-1 / R)
def compute(self, inputs, outputs):
deltaV = inputs['V_in'] - inputs['V_out']
outputs['I'] = deltaV / self.options['R']
class Diode(om.ExplicitComponent):
"""Computes current across a diode using the Shockley diode equation."""
def initialize(self):
self.options.declare('Is', default=1e-15, desc='Saturation current in Amps')
self.options.declare('Vt', default=.025875, desc='Thermal voltage in Volts')
def setup(self):
self.add_input('V_in', units='V')
self.add_input('V_out', units='V')
self.add_output('I', units='A')
# non-linear component, so we'll declare the partials here but compute them in compute_partials
self.declare_partials('I', 'V_in')
self.declare_partials('I', 'V_out')
def compute(self, inputs, outputs):
deltaV = inputs['V_in'] - inputs['V_out']
Is = self.options['Is']
Vt = self.options['Vt']
outputs['I'] = Is * (np.exp(deltaV / Vt) - 1)
def compute_partials(self, inputs, J):
deltaV = inputs['V_in'] - inputs['V_out']
Is = self.options['Is']
Vt = self.options['Vt']
I = Is * np.exp(deltaV / Vt)
J['I', 'V_in'] = I/Vt
J['I', 'V_out'] = -I/Vt
class Node(om.ImplicitComponent):
"""Computes voltage residual across a node based on incoming and outgoing current."""
def initialize(self):
self.options.declare('n_in', default=1, types=int, desc='number of connections with + assumed in')
self.options.declare('n_out', default=1, types=int, desc='number of current connections + assumed out')
def setup(self):
self.add_output('V', val=5., units='V')
for i in range(self.options['n_in']):
i_name = 'I_in:{}'.format(i)
self.add_input(i_name, units='A')
self.declare_partials('V', i_name, val=1)
for i in range(self.options['n_out']):
i_name = 'I_out:{}'.format(i)
self.add_input(i_name, units='A')
self.declare_partials('V', i_name, val=-1)
# note: we don't declare any partials wrt `V` here,
# because the residual doesn't directly depend on it
def apply_nonlinear(self, inputs, outputs, residuals):
residuals['V'] = 0.
for i_conn in range(self.options['n_in']):
residuals['V'] += inputs['I_in:{}'.format(i_conn)]
for i_conn in range(self.options['n_out']):
residuals['V'] -= inputs['I_out:{}'.format(i_conn)]
class Circuit(om.Group):
def setup(self):
self.add_subsystem('n1', Node(n_in=1, n_out=2), promotes_inputs=[('I_in:0', 'I_in')])
self.add_subsystem('n2', Node()) # leaving defaults
self.add_subsystem('R1', Resistor(R=100.), promotes_inputs=[('V_out', 'Vg')])
self.add_subsystem('R2', Resistor(R=10000.))
self.add_subsystem('D1', Diode(), promotes_inputs=[('V_out', 'Vg')])
self.connect('n1.V', ['R1.V_in', 'R2.V_in'])
self.connect('R1.I', 'n1.I_out:0')
self.connect('R2.I', 'n1.I_out:1')
self.connect('n2.V', ['R2.V_out', 'D1.V_in'])
self.connect('R2.I', 'n2.I_in:0')
self.connect('D1.I', 'n2.I_out:0')
self.nonlinear_solver = om.NewtonSolver()
self.nonlinear_solver.options['iprint'] = 2
self.nonlinear_solver.options['maxiter'] = 10
self.nonlinear_solver.options['solve_subsystems'] = True
p = om.Problem()
model = p.model
model.add_subsystem('circuit', Circuit())
p.setup()
p.set_val('circuit.I_in', 0.1)
p.set_val('circuit.Vg', 0.)
# set some initial guesses
p.set_val('circuit.n1.V', 10.)
p.set_val('circuit.n2.V', 1e-3)
p.run_model()
=======
circuit
=======
NL: Newton 0 ; 0.00141407214 1
NL: Newton 1 ; 0.0014214429 1.00521244
NL: Newton 2 ; 0.00142899707 1.01055457
NL: Newton 3 ; 0.00143666458 1.01597687
NL: Newton 4 ; 0.00144444604 1.02147974
NL: Newton 5 ; 0.00145234276 1.02706412
NL: Newton 6 ; 0.00146035606 1.03273095
NL: Newton 7 ; 0.0014684873 1.03848118
NL: Newton 8 ; 0.00147673781 1.04431575
NL: Newton 9 ; 0.00148510896 1.05023564
NL: Newton 10 ; 0.00149360212 1.05624181
NL: NewtonSolver 'NL: Newton' on system 'circuit' failed to converge in 10 iterations.
Simple implicit calculation doesn’t converge¶
Here we have a new system based on Kepler’s Equation. The solver isn’t able to converge with this setup. What can you change to get the model to converge?
Click here to reveal a small hint
Take a close look at the code output for any hints or warnings.
import numpy as np
import openmdao.api as om
prob = om.Problem()
bal = om.BalanceComp()
bal.add_balance(name='E', val=0.0, units='rad', eq_units='rad', rhs_name='M')
# Use M (mean anomaly) as the initial guess for E (eccentric anomaly)
def guess_function(inputs, outputs, residuals):
if np.abs(residuals['E']) > 1.0E-2:
outputs['E'] = inputs['M']
bal.options['guess_func'] = guess_function
# ExecComp used to compute the LHS of Kepler's equation.
lhs_comp = om.ExecComp('lhs=E - ecc * sin(E)',
lhs={'units': 'rad'},
E={'units': 'rad'},
ecc={'units': None})
prob.model.set_input_defaults('M', 85.0, units='deg')
prob.model.add_subsystem(name='lhs_comp', subsys=lhs_comp,
promotes_inputs=['E', 'ecc'])
prob.model.add_subsystem(name='balance', subsys=bal,
promotes_inputs=['M'],
promotes_outputs=['E'])
# Explicit connections
prob.model.connect('lhs_comp.lhs', 'balance.lhs:E')
# Set up solvers
prob.model.linear_solver = om.LinearBlockGS()
prob.model.nonlinear_solver = om.NewtonSolver(solve_subsystems=False, maxiter=10, iprint=2)
prob.setup()
prob.set_val('M', 85.0, units='deg')
prob.set_val('E', 85.0, units='deg')
prob.set_val('ecc', 0.6)
prob.run_model()
M = prob.get_val('M')
E = prob.get_val('E')
print(f'M = {M}')
print(f'E = {E}')
NL: Newton 0 ; 0.3321557 1
| LN: LNBGSSolver 'LN: LNBGS' on system '' failed to converge in 10 iterations.
NL: Newton 1 ; 0.237871778 0.716145404
| LN: LNBGSSolver 'LN: LNBGS' on system '' failed to converge in 10 iterations.
NL: Newton 2 ; 0.138431564 0.41676709
| LN: LNBGSSolver 'LN: LNBGS' on system '' failed to converge in 10 iterations.
NL: Newton 3 ; 0.078332712 0.235831305
| LN: LNBGSSolver 'LN: LNBGS' on system '' failed to converge in 10 iterations.
NL: Newton 4 ; 0.0434105989 0.130693524
| LN: LNBGSSolver 'LN: LNBGS' on system '' failed to converge in 10 iterations.
NL: Newton 5 ; 0.0233532448 0.070308126
| LN: LNBGSSolver 'LN: LNBGS' on system '' failed to converge in 10 iterations.
NL: Newton 6 ; 0.0121412104 0.0365527685
| LN: LNBGSSolver 'LN: LNBGS' on system '' failed to converge in 10 iterations.
NL: Newton 7 ; 0.00610725685 0.0183867291
| LN: LNBGSSolver 'LN: LNBGS' on system '' failed to converge in 10 iterations.
NL: Newton 8 ; 0.00298686556 0.00899236582
| LN: LNBGSSolver 'LN: LNBGS' on system '' failed to converge in 10 iterations.
NL: Newton 9 ; 0.00142908564 0.00430245707
| LN: LNBGSSolver 'LN: LNBGS' on system '' failed to converge in 10 iterations.
NL: Newton 10 ; 0.000672860219 0.00202573739
NL: NewtonSolver 'NL: Newton' on system '' failed to converge in 10 iterations.
M = [1.48352986]
E = [2.02265228]
Propulsor Newton solver doesn’t converge¶
This next problem uses the pyCycle engine modeling tool to create complex realistic problems. You will need to install pyCycle to work through this notebook. You can do so with pip install om-pycycle.
You don’t strictly need to understand pyCycle to find value from this notebook. Of course if you understand the physics, it might be easier to understand the solver hierarchy and why the model is set up that way. But if you don’t understand it, then you can view the states and residuals just as numbers without a particular physical meaning.
Let’s first start with a case where a Newton solver doesn’t converge. First we’ll define some pyCycle groups for a propulsor.
import openmdao.api as om
import pycycle.api as pyc
class Propulsor(pyc.Cycle):
def setup(self):
design = self.options['design']
USE_TABULAR = True
if USE_TABULAR:
self.options['thermo_method'] = 'TABULAR'
self.options['thermo_data'] = pyc.AIR_JETA_TAB_SPEC
else:
self.options['thermo_method'] = 'CEA'
self.options['thermo_data'] = pyc.species_data.janaf
FUEL_TYPE = 'JP-7'
self.add_subsystem('fc', pyc.FlightConditions())
self.add_subsystem('inlet', pyc.Inlet())
self.add_subsystem('fan', pyc.Compressor(map_data=pyc.FanMap, map_extrap=True))
self.add_subsystem('nozz', pyc.Nozzle())
self.add_subsystem('perf', pyc.Performance(num_nozzles=1, num_burners=0))
balance = om.BalanceComp()
if design:
self.add_subsystem('shaft', om.IndepVarComp('Nmech', 1., units='rpm'))
self.connect('shaft.Nmech', 'fan.Nmech')
balance.add_balance('W', units='lbm/s', eq_units='hp', val=50., lower=1., upper=500.)
self.add_subsystem('balance', balance,
promotes_inputs=[('rhs:W', 'pwr_target')])
self.connect('fan.power', 'balance.lhs:W')
else:
# vary mass flow till the nozzle area matches the design values
balance.add_balance('W', units='lbm/s', eq_units='inch**2', val=50, lower=1., upper=500.)
self.connect('nozz.Throat:stat:area', 'balance.lhs:W')
balance.add_balance('Nmech', val=1., units='rpm', lower=0.1, upper=2.0, eq_units='hp')
self.connect('balance.Nmech', 'fan.Nmech')
self.connect('fan.power', 'balance.lhs:Nmech')
self.add_subsystem('balance', balance,
promotes_inputs=[('rhs:Nmech', 'pwr_target')])
self.pyc_connect_flow('fc.Fl_O', 'inlet.Fl_I')
self.pyc_connect_flow('inlet.Fl_O', 'fan.Fl_I')
self.pyc_connect_flow('fan.Fl_O', 'nozz.Fl_I')
self.connect('fc.Fl_O:stat:P', 'nozz.Ps_exhaust')
self.connect('inlet.Fl_O:tot:P', 'perf.Pt2')
self.connect('fan.Fl_O:tot:P', 'perf.Pt3')
self.connect('inlet.F_ram', 'perf.ram_drag')
self.connect('nozz.Fg', 'perf.Fg_0')
self.connect('balance.W', 'fc.W')
newton = self.nonlinear_solver = om.NewtonSolver()
newton.options['atol'] = 1e-12
newton.options['rtol'] = 1e-12
newton.options['iprint'] = 2
newton.options['maxiter'] = 10
newton.options['solve_subsystems'] = True
newton.options['max_sub_solves'] = 10
newton.options['reraise_child_analysiserror'] = False
newton.linesearch = om.BoundsEnforceLS()
self.linear_solver = om.DirectSolver()
# base_class setup should be called as the last thing in your setup
super().setup()
class MPpropulsor(pyc.MPCycle):
def setup(self):
design = self.pyc_add_pnt('design', Propulsor(design=True, thermo_method='CEA'))
self.pyc_add_cycle_param('pwr_target', 100.)
# define the off-design conditions we want to run
self.od_pts = ['off_design']
self.od_MNs = [0.8,]
self.od_alts = [15000,]
self.od_Rlines = [2.2,]
for i, pt in enumerate(self.od_pts):
self.pyc_add_pnt(pt, Propulsor(design=False, thermo_method='CEA'))
self.set_input_defaults(pt+'.fc.MN', val=self.od_MNs[i])
self.set_input_defaults(pt+'.fc.alt', val=self.od_alts, units='m')
self.set_input_defaults(pt+'.fan.map.RlineMap', val=self.od_Rlines[i])
self.pyc_use_default_des_od_conns()
self.pyc_connect_des_od('nozz.Throat:stat:area', 'balance.rhs:W')
super().setup()
Next, let’s run the propulser through a design and off-design case to see what happens.
import numpy as np
prob = om.Problem()
prob.model = mp_propulsor = MPpropulsor()
prob.setup()
#Define the design point
prob.set_val('design.fc.alt', 10000, units='m')
prob.set_val('design.fc.MN', 0.8)
prob.set_val('design.inlet.MN', 0.6)
prob.set_val('design.fan.PR', 1.2)
prob.set_val('pwr_target', -3486.657, units='hp')
prob.set_val('design.fan.eff', 0.96)
# Set initial guesses for balances
prob['design.balance.W'] = 500.
for i, pt in enumerate(mp_propulsor.od_pts):
# initial guesses
prob['off_design.fan.PR'] = 3.
prob['off_design.balance.W'] = 10.
prob['off_design.balance.Nmech'] = 20. # normalized value
prob.set_solver_print(level=-1)
prob.set_solver_print(level=2, depth=2)
prob.model.design.nonlinear_solver.options['atol'] = 1e-6
prob.model.design.nonlinear_solver.options['rtol'] = 1e-6
prob.model.off_design.nonlinear_solver.options['atol'] = 1e-6
prob.model.off_design.nonlinear_solver.options['rtol'] = 1e-6
prob.run_model()
======
design
======
NL: Newton 0 ; 0.468905936 1
| LS: BCHK 0 ; 50.7119874 108.149596
NL: Newton 1 ; 0.0900132352 0.191964376
| LS: BCHK 0 ; 0.00611530407 0.0679378322
NL: Newton 2 ; 0.00611400388 0.0130388707
| LS: BCHK 0 ; 3.68697569e-05 0.00603037839
NL: Newton 3 ; 3.6869756e-05 7.86293225e-05
| LS: BCHK 0 ; 1.35988258e-09 3.68834169e-05
NL: Newton 4 ; 1.35927833e-09 2.89882943e-09
NL: Newton Converged
==========
off_design
==========
NL: Newton 0 ; 16.13985 1
| LS: BCHK 0 ; 98264889.5 6088339.71
NL: Newton 1 ; 2.66730484 0.165262059
| LS: BCHK 0 ; 3.21833742e+10 1.20658778e+10
NL: Newton 2 ; 64.0614724 3.96914918
| LS: BCHK 0 ; 3.48694839e+10 544312870
off_design.nozz.staticMN.ps_resid [-2189344.43307203] [7411947.00498394] [150.]
off_design.nozz.staticMN.ps_resid [-2189344.43307203] [7411947.00498394] [150.]
off_design.nozz.staticMN.ps_resid [-4.81875566] [287.05016442] [187.73045403]
off_design.nozz.staticMN.ps_resid [-4.81875566] [287.05016442] [187.73045403]
/home/john/anaconda3/envs/course/lib/python3.8/site-packages/openmdao/solvers/solver.py:651: SolverWarning:Your model has stalled three times and may be violating the bounds. In the future, turn on print_bound_enforce in your solver options here:
off_design.fan.ideal_flow.nonlinear_solver.linesearch.options['print_bound_enforce']=True.
The bound(s) being violated now are:
/home/john/anaconda3/envs/course/lib/python3.8/site-packages/openmdao/solvers/linesearch/backtracking.py:39: SolverWarning:'off_design.fan.ideal_flow.balance.T' exceeds lower bounds
Val: [-2.44706198e+11]
Lower: [150.]
/home/john/anaconda3/envs/course/lib/python3.8/site-packages/openmdao/solvers/solver.py:651: SolverWarning:Your model has stalled three times and may be violating the bounds. In the future, turn on print_bound_enforce in your solver options here:
off_design.fan.real_flow.nonlinear_solver.linesearch.options['print_bound_enforce']=True.
The bound(s) being violated now are:
/home/john/anaconda3/envs/course/lib/python3.8/site-packages/openmdao/solvers/linesearch/backtracking.py:35: SolverWarning:'off_design.fan.real_flow.balance.T' exceeds upper bounds
Val: [2613.1103812]
Upper: [2500.]
/mnt/c/Users/John/Dropbox/git/pyCycle/pycycle/elements/compressor.py:75: RuntimeWarning: invalid value encountered in log
outputs['eff_poly'] = Rt * np.log(PR) / ( Rt*np.log(PR) + S_out - S_in )
/mnt/c/Users/John/Dropbox/git/pyCycle/pycycle/thermo/static_ps_resid.py:128: RuntimeWarning: invalid value encountered in sqrt
Vsonic = (i['gamma']*i['R']*i['Ts'])**0.5
/mnt/c/Users/John/Dropbox/git/pyCycle/pycycle/thermo/static_ps_resid.py:128: RuntimeWarning: invalid value encountered in sqrt
Vsonic = (i['gamma']*i['R']*i['Ts'])**0.5
/mnt/c/Users/John/Dropbox/git/pyCycle/pycycle/thermo/static_ps_resid.py:128: RuntimeWarning: invalid value encountered in sqrt
Vsonic = (i['gamma']*i['R']*i['Ts'])**0.5
/mnt/c/Users/John/Dropbox/git/pyCycle/pycycle/thermo/static_ps_resid.py:128: RuntimeWarning: invalid value encountered in sqrt
Vsonic = (i['gamma']*i['R']*i['Ts'])**0.5
/mnt/c/Users/John/Dropbox/git/pyCycle/pycycle/thermo/static_ps_resid.py:128: RuntimeWarning: invalid value encountered in sqrt
Vsonic = (i['gamma']*i['R']*i['Ts'])**0.5
/mnt/c/Users/John/Dropbox/git/pyCycle/pycycle/thermo/static_ps_resid.py:128: RuntimeWarning: invalid value encountered in sqrt
Vsonic = (i['gamma']*i['R']*i['Ts'])**0.5
/home/john/anaconda3/envs/course/lib/python3.8/site-packages/openmdao/solvers/solver.py:651: SolverWarning:Your model has stalled three times and may be violating the bounds. In the future, turn on print_bound_enforce in your solver options here:
off_design.nozz.throat_total.nonlinear_solver.linesearch.options['print_bound_enforce']=True.
The bound(s) being violated now are:
/home/john/anaconda3/envs/course/lib/python3.8/site-packages/openmdao/solvers/linesearch/backtracking.py:35: SolverWarning:'off_design.nozz.throat_total.balance.T' exceeds upper bounds
Val: [2613.1103812]
Upper: [2500.]
/mnt/c/Users/John/Dropbox/git/pyCycle/pycycle/thermo/static_ps_resid.py:108: RuntimeWarning: invalid value encountered in sqrt
Vsonic = (i['gamma']*i['R']*i['Ts'])**0.5
/mnt/c/Users/John/Dropbox/git/pyCycle/pycycle/thermo/static_ps_resid.py:108: RuntimeWarning: invalid value encountered in sqrt
Vsonic = (i['gamma']*i['R']*i['Ts'])**0.5
/home/john/anaconda3/envs/course/lib/python3.8/site-packages/openmdao/solvers/solver.py:651: SolverWarning:Your model has stalled three times and may be violating the bounds. In the future, turn on print_bound_enforce in your solver options here:
off_design.nozz.staticMN.nonlinear_solver.linesearch.options['print_bound_enforce']=True.
The bound(s) being violated now are:
/home/john/anaconda3/envs/course/lib/python3.8/site-packages/openmdao/solvers/linesearch/backtracking.py:39: SolverWarning:'off_design.nozz.staticMN.ps_resid.Ps' exceeds lower bounds
Val: [-772646.27279157]
Lower: [-0.0009009]
/home/john/anaconda3/envs/course/lib/python3.8/site-packages/openmdao/solvers/linesearch/backtracking.py:39: SolverWarning:'off_design.nozz.staticMN.ps_resid.area' exceeds lower bounds
Val: [-1.19971121e-09]
Lower: [1.e-05]
NL: Newton 3 ; nan nan
NL: NewtonSolver 'NL: Newton' on system 'off_design': residuals contain 'inf' or 'NaN' after 3 iterations.
/mnt/c/Users/John/Dropbox/git/pyCycle/pycycle/elements/compressor.py:75: RuntimeWarning: invalid value encountered in log
outputs['eff_poly'] = Rt * np.log(PR) / ( Rt*np.log(PR) + S_out - S_in )
/mnt/c/Users/John/Dropbox/git/pyCycle/pycycle/thermo/static_ps_resid.py:128: RuntimeWarning: invalid value encountered in sqrt
Vsonic = (i['gamma']*i['R']*i['Ts'])**0.5
The design case converges well, but the off-design Newton solver fails after a few iterations. We get a whole bunch of warnings and some hints about what could be going on. How can you get this solver to converge?
Click here to reveal a big hint
Try some different initial guesses. How can you know which ones need changing?