Solutions: Debugging solvers – a hands-on set of problems
Contents
Solutions: Debugging solvers – a hands-on set of problems¶
Level: Advanced
Throughout this notebook I will mark specific important changes to code using a # Solution! comment.
Circuit example with a Newton solver¶
In this first case, we had a nonlinear solver added but no linear solver! Because we are using a Newton solver, you need a linear solver to compute the derivative information to be used during the nonlinear solve. Here I’ve added a DirectSolver.
That alone does not make the model converge; we also need to add a linesearch to the Newton method to help move through the solution space correctly. With these two changes we see convergence.
In the original problem setup, you might’ve noticed the Newton solver wasn’t able to reduce the residuals even from the beginning. This is sometimes a sign that the derivative information is incorrect (unsolved linear system) because the solver doesn’t know where to move in the solution space.
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.linear_solver = om.DirectSolver() # Solution!
self.nonlinear_solver.options['iprint'] = 2
self.nonlinear_solver.options['maxiter'] = 10
self.nonlinear_solver.options['solve_subsystems'] = True
self.nonlinear_solver.linesearch = om.ArmijoGoldsteinLS() # Solution!
self.nonlinear_solver.linesearch.options['maxiter'] = 10
self.nonlinear_solver.linesearch.options['iprint'] = 2
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
| LS: AG 1 ; 6.9707252e+152 1
| LS: AG 2 ; 8.51199079e+68 0.5
| LS: AG 3 ; 9.40605982e+26 0.25
| LS: AG 4 ; 988771.709 0.125
| LS: AG 5 ; 0.00130322117 0.0625
NL: Newton 1 ; 0.00130322117 0.921608691
| LS: AG 1 ; 5580826.07 1
| LS: AG 2 ; 13.3748871 0.5
| LS: AG 3 ; 0.0198015756 0.25
| LS: AG 4 ; 0.000828003297 0.125
NL: Newton 2 ; 0.000828003297 0.585545301
NL: Newton 3 ; 7.02986117e-06 0.00497135964
NL: Newton 4 ; 2.64550002e-08 1.87083809e-05
NL: Newton 5 ; 3.78431444e-13 2.67618202e-10
NL: Newton Converged
Simple implicit calculation doesn’t converge¶
In this problem the linear solver we were using (Linear Block Gauss-Seidel) doesn’t correctly solve the linear system. We saw some print-outs saying that the linear solver didn’t converge. By changing to a different linear solver, we’re able to converge the system.
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.ScipyKrylov() # Solution!
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
NL: Newton 1 ; 0.117134648 0.352649822
NL: Newton 2 ; 0.00208780933 0.00628563452
NL: Newton 3 ; 7.36738647e-07 2.2180521e-06
NL: Newton 4 ; 9.19264664e-14 2.76757155e-13
NL: Newton Converged
M = [1.48352986]
E = [2.02317564]
Propulsor Newton solver doesn’t converge¶
For this next case, it turns out we just needed to change the initial guess. The pyCycle model is quite touchy and heavily relies on good initial guesses to converge. By changing W (the mass flow rate) through the engine, we’re able to converge this system correctly.
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'] = 100. # Solution!
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.0863489 1
| LS: BCHK 0 ; 13215486.8 821534.263
NL: Newton 1 ; 0.480750635 0.0298856278
| LS: BCHK 0 ; 612567.884 1274190.48
NL: Newton 2 ; 0.0511605572 0.00318037098
| LS: BCHK 0 ; 38285.792 748345.877
NL: Newton 3 ; 0.00128938347 8.01538918e-05
| LS: BCHK 0 ; 187.538455 145448.161
NL: Newton 4 ; 0.000731094371 4.54481235e-05
| LS: BCHK 0 ; 22.2153148 30386.3847
NL: Newton 5 ; 1.23203982e-05 7.65891516e-07
NL: Newton Converged