diff --git a/appveyor.yml b/appveyor.yml index 7effaac6..1b5d47a8 100644 --- a/appveyor.yml +++ b/appveyor.yml @@ -7,7 +7,6 @@ platform: environment: matrix: - - PYTHON: 2.7 - PYTHON: 3.6 init: diff --git a/examples/B738.py b/examples/B738.py index fc0d630b..77f58a9a 100644 --- a/examples/B738.py +++ b/examples/B738.py @@ -75,8 +75,7 @@ def setup(self): # airplanes which consume fuel will need to integrate # fuel usage across the mission and subtract it from TOW - nn_simpson = int((nn - 1) / 2) - self.add_subsystem('intfuel', Integrator(num_intervals=nn_simpson, method='simpson', + self.add_subsystem('intfuel', Integrator(num_nodes=nn, method='simpson', quantity_units='kg', diff_units='s', time_setup='duration'), promotes_inputs=[('dqdt', 'fuel_flow'), 'duration', @@ -97,7 +96,7 @@ def setup(self): nn = 11 # Define a bunch of design varaiables and airplane-specific parameters - dv_comp = self.add_subsystem('dv_comp', DictIndepVarComp(acdata, seperator='|'), + dv_comp = self.add_subsystem('dv_comp', DictIndepVarComp(acdata), promotes_outputs=["*"]) dv_comp.add_output_from_dict('ac|aero|CLmax_TO') dv_comp.add_output_from_dict('ac|aero|polar|e') diff --git a/examples/Caravan.py b/examples/Caravan.py index bf70d997..9aaf7957 100644 --- a/examples/Caravan.py +++ b/examples/Caravan.py @@ -66,8 +66,7 @@ def setup(self): # airplanes which consume fuel will need to integrate # fuel usage across the mission and subtract it from TOW - nn_simpson = int((nn - 1) / 2) - self.add_subsystem('intfuel', Integrator(num_intervals=nn_simpson, method='simpson', + self.add_subsystem('intfuel', Integrator(num_nodes=nn, method='simpson', quantity_units='kg', diff_units='s', time_setup='duration'), promotes_inputs=[('dqdt', 'fuel_flow'), 'duration', @@ -89,7 +88,7 @@ def setup(self): nn = 11 # Define a bunch of design varaiables and airplane-specific parameters - dv_comp = self.add_subsystem('dv_comp', DictIndepVarComp(acdata, seperator='|'), + dv_comp = self.add_subsystem('dv_comp', DictIndepVarComp(acdata), promotes_outputs=["*"]) dv_comp.add_output_from_dict('ac|aero|CLmax_TO') dv_comp.add_output_from_dict('ac|aero|polar|e') diff --git a/examples/ElectricSinglewithThermal.py b/examples/ElectricSinglewithThermal.py index b5e262d0..82a8dd12 100644 --- a/examples/ElectricSinglewithThermal.py +++ b/examples/ElectricSinglewithThermal.py @@ -71,7 +71,7 @@ class ElectricTBMAnalysisGroup(Group): def setup(self): nn = 11 - dv_comp = self.add_subsystem('dv_comp',DictIndepVarComp(acdata,seperator='|'),promotes_outputs=["*"]) + dv_comp = self.add_subsystem('dv_comp',DictIndepVarComp(acdata),promotes_outputs=["*"]) #eventually replace the following aerodynamic parameters with an analysis module (maybe OpenAeroStruct) dv_comp.add_output_from_dict('ac|aero|CLmax_TO') dv_comp.add_output_from_dict('ac|aero|polar|e') diff --git a/examples/HybridTwin.py b/examples/HybridTwin.py index 9ffd7103..946e1726 100644 --- a/examples/HybridTwin.py +++ b/examples/HybridTwin.py @@ -88,8 +88,7 @@ def setup(self): self.connect('propmodel.eng1.component_weight', 'W_engine') self.connect('propmodel.gen1.component_weight', 'W_generator') self.connect('propmodel.motors_weight', 'W_motors') - nn_simpson = int((nn - 1) / 2) - self.add_subsystem('intfuel', Integrator(num_intervals=nn_simpson, method='simpson', + self.add_subsystem('intfuel', Integrator(num_nodes=nn, method='simpson', quantity_units='kg', diff_units='s', time_setup='duration'), promotes_inputs=[('dqdt', 'fuel_flow'), 'duration', @@ -110,7 +109,7 @@ def setup(self): nn = 11 # Define a bunch of design varaiables and airplane-specific parameters - dv_comp = self.add_subsystem('dv_comp', DictIndepVarComp(acdata, seperator='|'), + dv_comp = self.add_subsystem('dv_comp', DictIndepVarComp(acdata), promotes_outputs=["*"]) dv_comp.add_output_from_dict('ac|aero|CLmax_TO') dv_comp.add_output_from_dict('ac|aero|polar|e') diff --git a/examples/HybridTwin_thermal.py b/examples/HybridTwin_thermal.py index 6dac5e81..fefbd3d7 100644 --- a/examples/HybridTwin_thermal.py +++ b/examples/HybridTwin_thermal.py @@ -100,8 +100,7 @@ def setup(self): self.connect('propmodel.hx.frontal_area','hxadder.hx_frontal_area') self.connect('propmodel.area_nozzle','hxadder.nozzle_area') - nn_simpson = int((nn - 1) / 2) - self.add_subsystem('intfuel', Integrator(num_intervals=nn_simpson, method='simpson', + self.add_subsystem('intfuel', Integrator(num_nodes=nn, method='simpson', quantity_units='kg', diff_units='s', time_setup='duration'), promotes_inputs=[('dqdt', 'fuel_flow'), 'duration', @@ -122,7 +121,7 @@ def setup(self): nn = 11 # Define a bunch of design varaiables and airplane-specific parameters - dv_comp = self.add_subsystem('dv_comp', DictIndepVarComp(acdata, seperator='|'), + dv_comp = self.add_subsystem('dv_comp', DictIndepVarComp(acdata), promotes_outputs=["*"]) dv_comp.add_output_from_dict('ac|aero|CLmax_TO') dv_comp.add_output_from_dict('ac|aero|polar|e') diff --git a/examples/KingAirC90GT.py b/examples/KingAirC90GT.py index 27365aab..ac2a8145 100644 --- a/examples/KingAirC90GT.py +++ b/examples/KingAirC90GT.py @@ -66,9 +66,8 @@ def setup(self): self.connect('propmodel.engines_weight', 'W_engine') # airplanes which consume fuel will need to integrate - # fuel usage across the mission and subtract it from TOW - nn_simpson = int((nn - 1) / 2) - self.add_subsystem('intfuel', Integrator(num_intervals=nn_simpson, method='simpson', + # fuel usage across the mission and subtract it from TOW + self.add_subsystem('intfuel', Integrator(num_nodes=nn, method='simpson', quantity_units='kg', diff_units='s', time_setup='duration'), promotes_inputs=[('dqdt', 'fuel_flow'), 'duration', @@ -90,7 +89,7 @@ def setup(self): nn = 11 # Define a bunch of design varaiables and airplane-specific parameters - dv_comp = self.add_subsystem('dv_comp', DictIndepVarComp(acdata, seperator='|'), + dv_comp = self.add_subsystem('dv_comp', DictIndepVarComp(acdata), promotes_outputs=["*"]) dv_comp.add_output_from_dict('ac|aero|CLmax_TO') dv_comp.add_output_from_dict('ac|aero|polar|e') diff --git a/examples/TBM850.py b/examples/TBM850.py index 8b7bd817..336cd900 100644 --- a/examples/TBM850.py +++ b/examples/TBM850.py @@ -66,8 +66,7 @@ def setup(self): # airplanes which consume fuel will need to integrate # fuel usage across the mission and subtract it from TOW - nn_simpson = int((nn - 1) / 2) - self.add_subsystem('intfuel', Integrator(num_intervals=nn_simpson, method='simpson', + self.add_subsystem('intfuel', Integrator(num_nodes=nn, method='simpson', quantity_units='kg', diff_units='s', time_setup='duration'), promotes_inputs=[('dqdt', 'fuel_flow'), 'duration', @@ -89,7 +88,7 @@ def setup(self): nn = 11 # Define a bunch of design varaiables and airplane-specific parameters - dv_comp = self.add_subsystem('dv_comp', DictIndepVarComp(acdata, seperator='|'), + dv_comp = self.add_subsystem('dv_comp', DictIndepVarComp(acdata), promotes_outputs=["*"]) dv_comp.add_output_from_dict('ac|aero|CLmax_TO') dv_comp.add_output_from_dict('ac|aero|polar|e') diff --git a/openconcept/analysis/atmospherics/tests/test_atmospherics.py b/openconcept/analysis/atmospherics/tests/test_atmospherics.py index afb9af6b..7d4fe2d3 100644 --- a/openconcept/analysis/atmospherics/tests/test_atmospherics.py +++ b/openconcept/analysis/atmospherics/tests/test_atmospherics.py @@ -1,7 +1,7 @@ from __future__ import division import unittest import numpy as np -from openmdao.utils.assert_utils import assert_rel_error, assert_check_partials +from openmdao.utils.assert_utils import assert_near_equal, assert_check_partials from openmdao.api import IndepVarComp, Group, Problem from openconcept.analysis.atmospherics.compute_atmos_props import ComputeAtmosphericProperties @@ -27,22 +27,22 @@ def setUp(self): def test_sea_level_and_30kft(self): #check conditions at sea level - assert_rel_error(self, self.prob['fltcond|rho'][0],1.225,tolerance=1e-4) - assert_rel_error(self, self.prob['fltcond|p'][0],101325,tolerance=1e-3) - assert_rel_error(self, self.prob['fltcond|T'][0],288.15,tolerance=1e-3) - assert_rel_error(self, self.prob['fltcond|Utrue'][0],61.7333,tolerance=1e-3) - assert_rel_error(self, self.prob['fltcond|q'][0],2334.2398,tolerance=1e-3) - assert_rel_error(self, self.prob['fltcond|M'][0],0.1814,tolerance=1e-3) - assert_rel_error(self, self.prob['fltcond|a'][0],340.294,tolerance=1e-3) + assert_near_equal(self.prob['fltcond|rho'][0],1.225,tolerance=1e-4) + assert_near_equal(self.prob['fltcond|p'][0],101325,tolerance=1e-3) + assert_near_equal(self.prob['fltcond|T'][0],288.15,tolerance=1e-3) + assert_near_equal(self.prob['fltcond|Utrue'][0],61.7333,tolerance=1e-3) + assert_near_equal(self.prob['fltcond|q'][0],2334.2398,tolerance=1e-3) + assert_near_equal(self.prob['fltcond|M'][0],0.1814,tolerance=1e-3) + assert_near_equal(self.prob['fltcond|a'][0],340.294,tolerance=1e-3) #check conditions at 30kft (1976 standard atmosphere verified at https://www.digitaldutch.com/atmoscalc/) - assert_rel_error(self, self.prob['fltcond|rho'][-1],0.458312,tolerance=1e-4) - assert_rel_error(self, self.prob['fltcond|p'][-1],30089.6,tolerance=1e-3) - assert_rel_error(self, self.prob['fltcond|T'][-1],228.714,tolerance=1e-3) - assert_rel_error(self, self.prob['fltcond|Utrue'][-1],61.7333*np.sqrt(1.225/0.458312),tolerance=1e-3) - assert_rel_error(self, self.prob['fltcond|q'][-1],2334.2398,tolerance=1e-3) - assert_rel_error(self, self.prob['fltcond|M'][-1],0.3326,tolerance=1e-3) - assert_rel_error(self, self.prob['fltcond|a'][-1],303.2301,tolerance=1e-3) + assert_near_equal(self.prob['fltcond|rho'][-1],0.458312,tolerance=1e-4) + assert_near_equal(self.prob['fltcond|p'][-1],30089.6,tolerance=1e-3) + assert_near_equal(self.prob['fltcond|T'][-1],228.714,tolerance=1e-3) + assert_near_equal(self.prob['fltcond|Utrue'][-1],61.7333*np.sqrt(1.225/0.458312),tolerance=1e-3) + assert_near_equal(self.prob['fltcond|q'][-1],2334.2398,tolerance=1e-3) + assert_near_equal(self.prob['fltcond|M'][-1],0.3326,tolerance=1e-3) + assert_near_equal(self.prob['fltcond|a'][-1],303.2301,tolerance=1e-3) def test_partials(self): partials = self.prob.check_partials(method='cs', out_stream=None) @@ -55,13 +55,13 @@ def setUp(self): self.prob.run_model() def test_sea_level(self): - assert_rel_error(self, self.prob['fltcond|rho'][0],1.225,tolerance=1e-4) - assert_rel_error(self, self.prob['fltcond|p'][0],101325,tolerance=1e-3) - assert_rel_error(self, self.prob['fltcond|T'][0],288.15,tolerance=1e-3) - assert_rel_error(self, self.prob['fltcond|Utrue'][0],61.7333,tolerance=1e-3) - assert_rel_error(self, self.prob['fltcond|q'][0],2334.2398,tolerance=1e-3) - assert_rel_error(self, self.prob['fltcond|M'][0],0.1814,tolerance=1e-3) - assert_rel_error(self, self.prob['fltcond|a'][0],340.294,tolerance=1e-3) + assert_near_equal(self.prob['fltcond|rho'][0],1.225,tolerance=1e-4) + assert_near_equal(self.prob['fltcond|p'][0],101325,tolerance=1e-3) + assert_near_equal(self.prob['fltcond|T'][0],288.15,tolerance=1e-3) + assert_near_equal(self.prob['fltcond|Utrue'][0],61.7333,tolerance=1e-3) + assert_near_equal(self.prob['fltcond|q'][0],2334.2398,tolerance=1e-3) + assert_near_equal(self.prob['fltcond|M'][0],0.1814,tolerance=1e-3) + assert_near_equal(self.prob['fltcond|a'][0],340.294,tolerance=1e-3) def test_partials(self): partials = self.prob.check_partials(method='cs', out_stream=None) diff --git a/openconcept/analysis/performance/solver_phases.py b/openconcept/analysis/performance/solver_phases.py index 99c9abcb..d1832529 100644 --- a/openconcept/analysis/performance/solver_phases.py +++ b/openconcept/analysis/performance/solver_phases.py @@ -546,13 +546,13 @@ def setup(self): self.add_subsystem('haccel',HorizontalAcceleration(num_nodes=nn), promotes_inputs=['*'],promotes_outputs=['*']) - nn_simpson = int((nn-1)/2) + if flight_phase == 'v1v0': #unfortunately need to shoot backwards to avoid negative airspeeds #reverse the order of the accelerations so the last one is first (and make them negative) self.add_subsystem('flipaccel', FlipVectorComp(num_nodes=nn, units='m/s**2', negative=True), promotes_inputs=[('vec_in','accel_horiz')]) #integrate the timesteps in reverse from near zero speed. - self.add_subsystem('intvelocity',Integrator(num_intervals=nn_simpson, method='simpson', quantity_units='m/s', diff_units='s',time_setup='duration', lower=1.5), + self.add_subsystem('intvelocity',Integrator(num_nodes=nn, method='simpson', quantity_units='m/s', diff_units='s',time_setup='duration', lower=1.5), promotes_inputs=['duration',('q_initial','zero_speed')],promotes_outputs=[('q_final','fltcond|Utrue_initial')]) self.connect('flipaccel.vec_out','intvelocity.dqdt') #flip the result of the reverse integration again so the flight condition is forward and consistent with everythign else @@ -563,7 +563,7 @@ def setup(self): promotes_inputs=['*'],promotes_outputs=['duration']) else: # forward shooting for these acceleration segmentes - self.add_subsystem('intvelocity',Integrator(num_intervals=nn_simpson, method='simpson', quantity_units='m/s', diff_units='s', time_setup='duration', lower=1.5), + self.add_subsystem('intvelocity',Integrator(num_nodes=nn, method='simpson', quantity_units='m/s', diff_units='s', time_setup='duration', lower=1.5), promotes_inputs=[('dqdt','accel_horiz'),'duration',('q_initial','fltcond|Utrue_initial')],promotes_outputs=[('q','fltcond|Utrue'),('q_final','fltcond|Utrue_final')]) if flight_phase == 'v0v1': self.connect('zero_speed','fltcond|Utrue_initial') @@ -574,10 +574,10 @@ def setup(self): promotes_inputs=['*'],promotes_outputs=['duration']) if zero_start: - self.add_subsystem('intrange',Integrator(num_intervals=nn_simpson, method='simpson', quantity_units='m', diff_units='s',zero_start=zero_start, time_setup='duration'), + self.add_subsystem('intrange',Integrator(num_nodes=nn, method='simpson', quantity_units='m', diff_units='s',zero_start=zero_start, time_setup='duration'), promotes_inputs=[('dqdt','fltcond|groundspeed'),'duration'],promotes_outputs=[('q','range'),('q_final','range_final')]) else: - self.add_subsystem('intrange',Integrator(num_intervals=nn_simpson, method='simpson', quantity_units='m', diff_units='s',zero_start=zero_start, time_setup='duration'), + self.add_subsystem('intrange',Integrator(num_nodes=nn, method='simpson', quantity_units='m', diff_units='s',zero_start=zero_start, time_setup='duration'), promotes_inputs=[('dqdt','fltcond|groundspeed'),'duration',('q_initial','range_initial')],promotes_outputs=[('q','range'),('q_final','range_final')]) class RotationPhase(Group): @@ -673,17 +673,17 @@ def setup(self): self.add_subsystem('lift',Lift(num_nodes=nn), promotes_inputs=['*'],promotes_outputs=['*']) self.add_subsystem('haccel',HorizontalAcceleration(num_nodes=nn), promotes_inputs=['*'],promotes_outputs=['*']) self.add_subsystem('vaccel',VerticalAcceleration(num_nodes=nn), promotes_inputs=['*'],promotes_outputs=['*']) - nn_simpson = int((nn-1)/2) + # TODO always starts from zero altitude self.add_subsystem('clear_obstacle',BalanceComp(name='duration',units='s',val=1,eq_units='m',rhs_name='fltcond|h_final',lhs_name='h_obs',lower=0.1,upper=15), promotes_inputs=['*'],promotes_outputs=['duration']) - self.add_subsystem('intvelocity',Integrator(num_intervals=nn_simpson, method='simpson', quantity_units='m/s', diff_units='s',time_setup='duration',lower=0.1), + self.add_subsystem('intvelocity',Integrator(num_nodes=nn, method='simpson', quantity_units='m/s', diff_units='s',time_setup='duration',lower=0.1), promotes_inputs=[('dqdt','accel_horiz'),'duration',('q_initial','fltcond|Utrue_initial')],promotes_outputs=[('q','fltcond|Utrue'),('q_final','fltcond|Utrue_final')]) - self.add_subsystem('intrange',Integrator(num_intervals=nn_simpson, method='simpson', quantity_units='m', diff_units='s', time_setup='duration'), + self.add_subsystem('intrange',Integrator(num_nodes=nn, method='simpson', quantity_units='m', diff_units='s', time_setup='duration'), promotes_inputs=[('dqdt','fltcond|groundspeed'),'duration',('q_initial','range_initial')],promotes_outputs=[('q','range'),('q_final','range_final')]) - self.add_subsystem('intvs',Integrator(num_intervals=nn_simpson, method='simpson', quantity_units='m/s', diff_units='s', time_setup='duration',zero_start=True), + self.add_subsystem('intvs',Integrator(num_nodes=nn, method='simpson', quantity_units='m/s', diff_units='s', time_setup='duration',zero_start=True), promotes_inputs=[('dqdt','accel_vert'),'duration'],promotes_outputs=[('q','fltcond|vs'),('q_final','fltcond|vs_final')]) - self.add_subsystem('inth',Integrator(num_intervals=nn_simpson, method='simpson', quantity_units='m', diff_units='s', time_setup='duration',zero_start=True), + self.add_subsystem('inth',Integrator(num_nodes=nn, method='simpson', quantity_units='m', diff_units='s', time_setup='duration',zero_start=True), promotes_inputs=[('dqdt','fltcond|vs'),'duration'],promotes_outputs=[('q','fltcond|h'),('q_final','fltcond|h_final')]) class SteadyFlightPhase(Group): @@ -762,8 +762,8 @@ def setup(self): ivcomp.add_output('fltcond|Ueas',val=np.ones((nn,))*90, units='m/s') ivcomp.add_output('fltcond|vs',val=np.ones((nn,))*1, units='m/s') ivcomp.add_output('zero_accel',val=np.zeros((nn,)),units='m/s**2') - nn_simpson = int((nn-1)/2) - self.add_subsystem('inth',Integrator(num_intervals=nn_simpson, method='simpson', quantity_units='m', diff_units='s', time_setup='duration'), + + self.add_subsystem('inth',Integrator(num_nodes=nn, method='simpson', quantity_units='m', diff_units='s', time_setup='duration'), promotes_inputs=[('dqdt','fltcond|vs'),'duration',('q_initial','fltcond|h_initial')],promotes_outputs=[('q','fltcond|h'),('q_final','fltcond|h_final')]) self.add_subsystem('atmos', ComputeAtmosphericProperties(num_nodes=nn, true_airspeed_in=False), promotes_inputs=['*'], promotes_outputs=['*']) self.add_subsystem('gs',Groundspeeds(num_nodes=nn),promotes_inputs=['*'],promotes_outputs=['*']) @@ -773,7 +773,7 @@ def setup(self): self.add_subsystem('lift',Lift(num_nodes=nn), promotes_inputs=['*'],promotes_outputs=['*']) self.add_subsystem('haccel',HorizontalAcceleration(num_nodes=nn), promotes_inputs=['*'],promotes_outputs=['*']) - self.add_subsystem('intrange',Integrator(num_intervals=nn_simpson, method='simpson', quantity_units='m', diff_units='s', time_setup='duration'), + self.add_subsystem('intrange',Integrator(num_nodes=nn, method='simpson', quantity_units='m', diff_units='s', time_setup='duration'), promotes_inputs=[('dqdt','fltcond|groundspeed'),'duration',('q_initial','range_initial')],promotes_outputs=[('q','range'),('q_final','range_final')]) diff --git a/openconcept/analysis/performance/tests/test_solver_phase_helpers.py b/openconcept/analysis/performance/tests/test_solver_phase_helpers.py index af0bfa8a..1955ebf4 100644 --- a/openconcept/analysis/performance/tests/test_solver_phase_helpers.py +++ b/openconcept/analysis/performance/tests/test_solver_phase_helpers.py @@ -1,7 +1,7 @@ from __future__ import division import unittest import numpy as np -from openmdao.utils.assert_utils import assert_rel_error, assert_check_partials +from openmdao.utils.assert_utils import assert_near_equal, assert_check_partials from openmdao.api import IndepVarComp, Group, Problem from openconcept.analysis.performance.solver_phases import ClimbAngleComp, FlipVectorComp, Groundspeeds, HorizontalAcceleration, VerticalAcceleration, SteadyFlightCL, TakeoffTransition @@ -30,12 +30,12 @@ def setUp(self): self.prob.run_model() def test_level_flight(self): - assert_rel_error(self, self.prob['gamma'][0],0,tolerance=1e-10) + assert_near_equal(self.prob['gamma'][0],0,tolerance=1e-10) def test_climb_flight(self): self.prob['thrust'] = np.ones((1,))*1200 self.prob.run_model() - assert_rel_error(self, self.prob['gamma'][0], np.arcsin(200 / 1000 / g), tolerance=1e-10) + assert_near_equal(self.prob['gamma'][0], np.arcsin(200 / 1000 / g), tolerance=1e-10) def test_partials(self): partials = self.prob.check_partials(method='cs', out_stream=None) @@ -67,7 +67,7 @@ def setUp(self): self.prob.run_model() def test_flip_vec_order(self): - assert_rel_error(self, self.prob['vec_out'], np.linspace(100, 0, 11), tolerance=1e-10) + assert_near_equal(self.prob['vec_out'], np.linspace(100, 0, 11), tolerance=1e-10) def test_partials(self): partials = self.prob.check_partials(method='cs', out_stream=None) @@ -80,7 +80,7 @@ def setUp(self): self.prob.run_model() def test_flip_vec_order(self): - assert_rel_error(self, self.prob['vec_out'], np.zeros((1,)), tolerance=1e-10) + assert_near_equal(self.prob['vec_out'], np.zeros((1,)), tolerance=1e-10) def test_partials(self): partials = self.prob.check_partials(method='cs', out_stream=None) @@ -93,7 +93,7 @@ def setUp(self): self.prob.run_model() def test_flip_vec_order(self): - assert_rel_error(self, self.prob['vec_out'], np.linspace(-100, 0, 11), tolerance=1e-10) + assert_near_equal(self.prob['vec_out'], np.linspace(-100, 0, 11), tolerance=1e-10) def test_partials(self): partials = self.prob.check_partials(method='cs', out_stream=None) @@ -121,15 +121,15 @@ def setUp(self): self.prob.run_model() def test_level_flight(self): - assert_rel_error(self, self.prob['fltcond|groundspeed'][0],50,tolerance=1e-10) - assert_rel_error(self, self.prob['fltcond|cosgamma'][0],1.,tolerance=1e-10) - assert_rel_error(self, self.prob['fltcond|singamma'][0],0.,tolerance=1e-10) + assert_near_equal(self.prob['fltcond|groundspeed'][0],50,tolerance=1e-10) + assert_near_equal(self.prob['fltcond|cosgamma'][0],1.,tolerance=1e-10) + assert_near_equal(self.prob['fltcond|singamma'][0],0.,tolerance=1e-10) def test_climb_flight(self): gs = np.sqrt(50**2 - 3**2) - assert_rel_error(self, self.prob['fltcond|groundspeed'][-1],gs,tolerance=1e-10) - assert_rel_error(self, self.prob['fltcond|cosgamma'][-1],gs/50.,tolerance=1e-10) - assert_rel_error(self, self.prob['fltcond|singamma'][-1],3./50.,tolerance=1e-10) + assert_near_equal(self.prob['fltcond|groundspeed'][-1],gs,tolerance=1e-10) + assert_near_equal(self.prob['fltcond|cosgamma'][-1],gs/50.,tolerance=1e-10) + assert_near_equal(self.prob['fltcond|singamma'][-1],3./50.,tolerance=1e-10) def test_partials(self): partials = self.prob.check_partials(method='cs', out_stream=None) @@ -159,7 +159,7 @@ def setUp(self): self.prob.run_model() def test_steady_level_flights(self): - assert_rel_error(self, self.prob['accel_horiz'], np.zeros((9,)), tolerance=1e-10) + assert_near_equal(self.prob['accel_horiz'], np.zeros((9,)), tolerance=1e-10) def test_partials(self): partials = self.prob.check_partials(method='cs', out_stream=None) @@ -174,7 +174,7 @@ def setUp(self): self.prob.run_model() def test_steady_climb_flights(self): - assert_rel_error(self, self.prob['accel_horiz'], np.zeros((9,)), tolerance=1e-10) + assert_near_equal(self.prob['accel_horiz'], np.zeros((9,)), tolerance=1e-10) def test_partials(self): partials = self.prob.check_partials(method='cs', out_stream=None) @@ -189,7 +189,7 @@ def setUp(self): self.prob.run_model() def test_steady_climb_flights(self): - assert_rel_error(self, self.prob['accel_horiz'], np.zeros((9,)), tolerance=1e-10) + assert_near_equal(self.prob['accel_horiz'], np.zeros((9,)), tolerance=1e-10) def test_partials(self): partials = self.prob.check_partials(method='cs', out_stream=None) @@ -214,7 +214,7 @@ def test_accel_with_braking(self): brakeforce = 0.03 * (weight-lift) slopeforce = weight * singamma accel_horz_actual = (thrust - drag - brakeforce - slopeforce) / mass - assert_rel_error(self, self.prob['accel_horiz'][0], accel_horz_actual, tolerance=1e-10) + assert_near_equal(self.prob['accel_horiz'][0], accel_horz_actual, tolerance=1e-10) def test_accel_with_braking_and_lift(self): drag = 50.0 @@ -226,7 +226,7 @@ def test_accel_with_braking_and_lift(self): brakeforce = 0.03 * (weight-lift) slopeforce = weight * singamma accel_horz_actual = (thrust - drag - brakeforce - slopeforce) / mass - assert_rel_error(self, self.prob['accel_horiz'][4], accel_horz_actual, tolerance=1e-10) + assert_near_equal(self.prob['accel_horiz'][4], accel_horz_actual, tolerance=1e-10) def test_accel_lift_exceeds_weight(self): drag = 50.0 @@ -238,7 +238,7 @@ def test_accel_lift_exceeds_weight(self): brakeforce = 0.0 slopeforce = weight * singamma accel_horz_actual = (thrust - drag - brakeforce - slopeforce) / mass - assert_rel_error(self, self.prob['accel_horiz'][-1], accel_horz_actual, tolerance=1e-10) + assert_near_equal(self.prob['accel_horiz'][-1], accel_horz_actual, tolerance=1e-10) def test_partials(self): partials = self.prob.check_partials(method='cs', out_stream=None) @@ -268,7 +268,7 @@ def setUp(self): self.prob.run_model() def test_steady_level_flights(self): - assert_rel_error(self, self.prob['accel_vert'], np.zeros((9,)), tolerance=1e-10) + assert_near_equal(self.prob['accel_vert'], np.zeros((9,)), tolerance=1e-10) def test_partials(self): partials = self.prob.check_partials(method='cs', out_stream=None) @@ -285,7 +285,7 @@ def setUp(self): self.prob.run_model() def test_steady_climbing_flight(self): - assert_rel_error(self, self.prob['accel_vert'], np.zeros((9,)), tolerance=1e-10) + assert_near_equal(self.prob['accel_vert'], np.zeros((9,)), tolerance=1e-10) def test_partials(self): partials = self.prob.check_partials(method='cs', out_stream=None) @@ -299,7 +299,7 @@ def setUp(self): self.prob.run_model() def test_unsteady_pullup(self): - assert_rel_error(self, self.prob['accel_vert'], (100./100.)*np.ones((9,)), tolerance=1e-10) + assert_near_equal(self.prob['accel_vert'], (100./100.)*np.ones((9,)), tolerance=1e-10) def test_partials(self): partials = self.prob.check_partials(method='cs', out_stream=None) @@ -327,7 +327,7 @@ def setUp(self): self.prob.run_model() def test_steady_level_flights(self): - assert_rel_error(self, self.prob['fltcond|CL'], np.ones((9,))*100*g/1000./10./1.0, tolerance=1e-10) + assert_near_equal(self.prob['fltcond|CL'], np.ones((9,))*100*g/1000./10./1.0, tolerance=1e-10) def test_partials(self): partials = self.prob.check_partials(method='cs', out_stream=None) @@ -341,7 +341,7 @@ def setUp(self): self.prob.run_model() def test_steady_level_flights(self): - assert_rel_error(self, self.prob['fltcond|CL'], np.ones((9,))*100*g/1000./10.*0.98, tolerance=1e-10) + assert_near_equal(self.prob['fltcond|CL'], np.ones((9,))*100*g/1000./10.*0.98, tolerance=1e-10) def test_partials(self): partials = self.prob.check_partials(method='cs', out_stream=None) diff --git a/openconcept/analysis/tests/test_aerodynamics.py b/openconcept/analysis/tests/test_aerodynamics.py index beb071c9..c5e55d59 100644 --- a/openconcept/analysis/tests/test_aerodynamics.py +++ b/openconcept/analysis/tests/test_aerodynamics.py @@ -1,6 +1,6 @@ import unittest import numpy as np -from openmdao.utils.assert_utils import assert_rel_error, assert_check_partials +from openmdao.utils.assert_utils import assert_near_equal, assert_check_partials from openmdao.api import IndepVarComp, Group, Problem from openconcept.analysis.aerodynamics import PolarDrag, Lift, StallSpeed @@ -31,8 +31,8 @@ def setUp(self): def test_drag_vectorial(self): drag_cl0 = 0.5*1.225*70**2 * 30 *(0.02 + 0**2 / np.pi / 0.8 / 15) drag_cl1p5 = 0.5*1.225*70**2 * 30 *(0.02 + 1.5**2 / np.pi / 0.8 / 15) - assert_rel_error(self, self.prob['drag'][0], drag_cl0, tolerance=1e-8) - assert_rel_error(self, self.prob['drag'][-1], drag_cl1p5, tolerance=1e-8) + assert_near_equal(self.prob['drag'][0], drag_cl0, tolerance=1e-8) + assert_near_equal(self.prob['drag'][-1], drag_cl1p5, tolerance=1e-8) def test_partials(self): @@ -47,7 +47,7 @@ def setUp(self): def test_drag_scalar(self): drag_cl0 = 0.5*1.225*70**2 * 30 *(0.02 + 0**2 / np.pi / 0.8 / 15) - assert_rel_error(self, self.prob['drag'], drag_cl0, tolerance=1e-8) + assert_near_equal(self.prob['drag'], drag_cl0, tolerance=1e-8) def test_partials(self): partials = self.prob.check_partials(method='cs', out_stream=None) @@ -77,8 +77,8 @@ def setUp(self): def test_lift_vectorial(self): lift_cl0 = 0 lift_cl1p5 = 0.5*1.225*70**2 * 30 * 1.5 - assert_rel_error(self, self.prob['lift'][-1], lift_cl0, tolerance=1e-8) - assert_rel_error(self, self.prob['lift'][0], lift_cl1p5, tolerance=1e-8) + assert_near_equal(self.prob['lift'][-1], lift_cl0, tolerance=1e-8) + assert_near_equal(self.prob['lift'][0], lift_cl1p5, tolerance=1e-8) def test_partials(self): partials = self.prob.check_partials(method='cs', out_stream=None) @@ -92,7 +92,7 @@ def setUp(self): def test_lift_scalar(self): lift_cl1p5 = 0.5*1.225*70**2 * 30 * 1.5 - assert_rel_error(self, self.prob['lift'], lift_cl1p5, tolerance=1e-8) + assert_near_equal(self.prob['lift'], lift_cl1p5, tolerance=1e-8) def test_partials(self): partials = self.prob.check_partials(method='cs', out_stream=None) @@ -119,7 +119,7 @@ def setUp(self): def test_stall_speed(self): vstall = np.sqrt(2 * 1000 * 9.80665 / 1.225 / 30 / 2.5) - assert_rel_error(self, self.prob['Vstall_eas'], vstall, tolerance=1e-8) + assert_near_equal(self.prob['Vstall_eas'], vstall, tolerance=1e-8) def test_partials(self): partials = self.prob.check_partials(method='cs', out_stream=None) diff --git a/openconcept/components/battery.py b/openconcept/components/battery.py index d0dd3360..ffea8662 100644 --- a/openconcept/components/battery.py +++ b/openconcept/components/battery.py @@ -80,8 +80,8 @@ def setup(self): self.add_subsystem('divider',ElementMultiplyDivideComp(output_name='dSOCdt',input_names=['elec_load','max_energy'],vec_size=[nn,1],scaling_factor=-1,divide=[False,True],input_units=['W','kJ']), promotes_inputs=['*'],promotes_outputs=['*']) - nn_simpson = int((nn-1)/2) - self.add_subsystem('intload',Integrator(num_intervals=nn_simpson, method='simpson', quantity_units=None, diff_units='s', time_setup='duration'), + + self.add_subsystem('intload',Integrator(num_nodes=nn, method='simpson', quantity_units=None, diff_units='s', time_setup='duration'), promotes_inputs=[('q_initial','batt_SOC_initial'),'duration',('dqdt','dSOCdt')], promotes_outputs=[('q','SOC'),('q_final','SOC_final')]) diff --git a/openconcept/components/tests/test_heat_exchanger.py b/openconcept/components/tests/test_heat_exchanger.py index 7932dcba..b3ee142f 100644 --- a/openconcept/components/tests/test_heat_exchanger.py +++ b/openconcept/components/tests/test_heat_exchanger.py @@ -1,7 +1,7 @@ from __future__ import division import unittest import numpy as np -from openmdao.utils.assert_utils import assert_rel_error, assert_check_partials +from openmdao.utils.assert_utils import assert_near_equal, assert_check_partials from openmdao.api import IndepVarComp, Group, Problem from openconcept.components.heat_exchanger import OffsetStripFinGeometry, OffsetStripFinData, HydraulicDiameterReynoldsNumber, OutletTemperatures, PressureDrop from openconcept.components.heat_exchanger import NusseltFromColburnJ, ConvectiveCoefficient, FinEfficiency, UAOverall, NTUMethod, CrossFlowNTUEffectiveness, NTUEffectivenessActualHeatTransfer @@ -70,11 +70,11 @@ def test_default_settings(self): prob = Problem(OSFGeometryTestGroup(num_nodes=1)) prob.setup(check=True,force_alloc_complex=True) prob.run_model() - assert_rel_error(self, prob['osfgeometry.dh_cold'], 0.002462541, tolerance=1e-6) - assert_rel_error(self, prob['heat_transfer'], 10020.13126, tolerance=1e-6 ) - assert_rel_error(self, prob['delta_p_cold'], -131.9862069, tolerance=1e-6 ) - assert_rel_error(self, prob['delta_p_hot'], -9112.282754, tolerance=1e-6 ) - assert_rel_error(self, prob['component_weight'], 1.147605, tolerance=1e-5 ) + assert_near_equal(prob['osfgeometry.dh_cold'], 0.002462541, tolerance=1e-6) + assert_near_equal(prob['heat_transfer'], 10020.13126, tolerance=1e-6 ) + assert_near_equal(prob['delta_p_cold'], -131.9862069, tolerance=1e-6 ) + assert_near_equal(prob['delta_p_hot'], -9112.282754, tolerance=1e-6 ) + assert_near_equal(prob['component_weight'], 1.147605, tolerance=1e-5 ) partials = prob.check_partials(method='cs',compact_print=True) assert_check_partials(partials) @@ -89,10 +89,10 @@ def test_default_settings(self): # cost_base=0)) # prob.setup(check=True,force_alloc_complex=True) # prob.run_model() - # assert_rel_error(self, prob.get_val('battery.heat_out', units='kW'), np.ones(10)*100*0.05, tolerance=1e-15) - # assert_rel_error(self, prob['battery.component_sizing_margin'], np.ones(10)/3, tolerance=1e-15) - # assert_rel_error(self, prob['battery.component_cost'], 10000, tolerance=1e-15) - # assert_rel_error(self, prob.get_val('battery.max_energy', units='W*h'), 500*100, tolerance=1e-15) + # assert_near_equal(prob.get_val('battery.heat_out', units='kW'), np.ones(10)*100*0.05, tolerance=1e-15) + # assert_near_equal(prob['battery.component_sizing_margin'], np.ones(10)/3, tolerance=1e-15) + # assert_near_equal(prob['battery.component_cost'], 10000, tolerance=1e-15) + # assert_near_equal(prob.get_val('battery.max_energy', units='W*h'), 500*100, tolerance=1e-15) # partials = prob.check_partials(method='cs',compact_print=True) # assert_check_partials(partials) diff --git a/openconcept/components/tests/test_simple_comps.py b/openconcept/components/tests/test_simple_comps.py index cf3846b9..d69f13c1 100644 --- a/openconcept/components/tests/test_simple_comps.py +++ b/openconcept/components/tests/test_simple_comps.py @@ -1,7 +1,7 @@ from __future__ import division import unittest import numpy as np -from openmdao.utils.assert_utils import assert_rel_error, assert_check_partials +from openmdao.utils.assert_utils import assert_near_equal, assert_check_partials from openmdao.api import IndepVarComp, Group, Problem from openconcept.components import SimpleBattery, SimpleGenerator, SimpleMotor, SimplePropeller, SimpleTurboshaft, PowerSplit @@ -159,10 +159,10 @@ def test_default_settings(self): prob = Problem(BatteryTestGroup(vec_size=10, use_defaults=True)) prob.setup(check=True,force_alloc_complex=True) prob.run_model() - assert_rel_error(self, prob['battery.heat_out'], np.ones(10)*100*0.0, tolerance=1e-15) - assert_rel_error(self, prob['battery.component_sizing_margin'], np.ones(10)*0.20, tolerance=1e-15) - assert_rel_error(self, prob['battery.component_cost'], 5001, tolerance=1e-15) - assert_rel_error(self, prob.get_val('battery.max_energy', units='W*h'), 300*100, tolerance=1e-15) + assert_near_equal(prob['battery.heat_out'], np.ones(10)*100*0.0, tolerance=1e-15) + assert_near_equal(prob['battery.component_sizing_margin'], np.ones(10)*0.20, tolerance=1e-15) + assert_near_equal(prob['battery.component_cost'], 5001, tolerance=1e-15) + assert_near_equal(prob.get_val('battery.max_energy', units='W*h'), 300*100, tolerance=1e-15) partials = prob.check_partials(method='cs',compact_print=True) assert_check_partials(partials) @@ -177,10 +177,10 @@ def test_nondefault_settings(self): cost_base=0)) prob.setup(check=True,force_alloc_complex=True) prob.run_model() - assert_rel_error(self, prob.get_val('battery.heat_out', units='kW'), np.ones(10)*100*0.05, tolerance=1e-15) - assert_rel_error(self, prob['battery.component_sizing_margin'], np.ones(10)/3, tolerance=1e-15) - assert_rel_error(self, prob['battery.component_cost'], 10000, tolerance=1e-15) - assert_rel_error(self, prob.get_val('battery.max_energy', units='W*h'), 500*100, tolerance=1e-15) + assert_near_equal(prob.get_val('battery.heat_out', units='kW'), np.ones(10)*100*0.05, tolerance=1e-15) + assert_near_equal(prob['battery.component_sizing_margin'], np.ones(10)/3, tolerance=1e-15) + assert_near_equal(prob['battery.component_cost'], 10000, tolerance=1e-15) + assert_near_equal(prob.get_val('battery.max_energy', units='W*h'), 500*100, tolerance=1e-15) partials = prob.check_partials(method='cs',compact_print=True) assert_check_partials(partials) @@ -192,12 +192,12 @@ def test_default_settings(self): use_defaults=True)) prob.setup(check=True,force_alloc_complex=True) prob.run_model() - assert_rel_error(self, prob.get_val('motor.shaft_power_out', units='kW'), np.ones(10)*90, tolerance=1e-15) - assert_rel_error(self, prob.get_val('motor.elec_load', units='kW'), np.ones(10)*90, tolerance=1e-15) - assert_rel_error(self, prob.get_val('motor.heat_out', units='kW'), np.ones(10)*0.0, tolerance=1e-15) - assert_rel_error(self, prob.get_val('motor.component_sizing_margin'), np.ones(10)*0.90, tolerance=1e-15) - assert_rel_error(self, prob.get_val('motor.component_cost', units='USD'), 13423.818791946309, tolerance=1e-10) - assert_rel_error(self, prob.get_val('motor.component_weight', units='kg'), 20, tolerance=1e-15) + assert_near_equal(prob.get_val('motor.shaft_power_out', units='kW'), np.ones(10)*90, tolerance=1e-15) + assert_near_equal(prob.get_val('motor.elec_load', units='kW'), np.ones(10)*90, tolerance=1e-15) + assert_near_equal(prob.get_val('motor.heat_out', units='kW'), np.ones(10)*0.0, tolerance=1e-15) + assert_near_equal(prob.get_val('motor.component_sizing_margin'), np.ones(10)*0.90, tolerance=1e-15) + assert_near_equal(prob.get_val('motor.component_cost', units='USD'), 13423.818791946309, tolerance=1e-10) + assert_near_equal(prob.get_val('motor.component_weight', units='kg'), 20, tolerance=1e-15) partials = prob.check_partials(method='cs',compact_print=True) assert_check_partials(partials) @@ -212,12 +212,12 @@ def test_nondefault_settings(self): )) prob.setup(check=True,force_alloc_complex=True) prob.run_model() - assert_rel_error(self, prob.get_val('motor.shaft_power_out', units='kW'), np.ones(10)*90*0.95, tolerance=1e-15) - assert_rel_error(self, prob.get_val('motor.elec_load', units='kW'), np.ones(10)*90, tolerance=1e-15) - assert_rel_error(self, prob.get_val('motor.heat_out', units='kW'), np.ones(10)*90*0.05, tolerance=1e-15) - assert_rel_error(self, prob.get_val('motor.component_sizing_margin'), np.ones(10)*0.90, tolerance=1e-15) - assert_rel_error(self, prob.get_val('motor.component_cost', units='USD'), 203, tolerance=1e-15) - assert_rel_error(self, prob.get_val('motor.component_weight', units='kg'), 35.333333333333333333, tolerance=1e-10) + assert_near_equal(prob.get_val('motor.shaft_power_out', units='kW'), np.ones(10)*90*0.95, tolerance=1e-15) + assert_near_equal(prob.get_val('motor.elec_load', units='kW'), np.ones(10)*90, tolerance=1e-15) + assert_near_equal(prob.get_val('motor.heat_out', units='kW'), np.ones(10)*90*0.05, tolerance=1e-15) + assert_near_equal(prob.get_val('motor.component_sizing_margin'), np.ones(10)*0.90, tolerance=1e-15) + assert_near_equal(prob.get_val('motor.component_cost', units='USD'), 203, tolerance=1e-15) + assert_near_equal(prob.get_val('motor.component_weight', units='kg'), 35.333333333333333333, tolerance=1e-10) partials = prob.check_partials(method='cs',compact_print=True) assert_check_partials(partials) @@ -228,11 +228,11 @@ def test_default_settings(self): use_defaults=True)) prob.setup(check=True,force_alloc_complex=True) prob.run_model() - assert_rel_error(self, prob.get_val('generator.elec_power_out', units='kW'), np.ones(10)*90, tolerance=1e-15) - assert_rel_error(self, prob.get_val('generator.heat_out', units='kW'), np.ones(10)*0.0, tolerance=1e-15) - assert_rel_error(self, prob.get_val('generator.component_sizing_margin'), np.ones(10)*0.90, tolerance=1e-15) - assert_rel_error(self, prob.get_val('generator.component_cost', units='USD'), 13423.818791946309, tolerance=1e-10) - assert_rel_error(self, prob.get_val('generator.component_weight', units='kg'), 20, tolerance=1e-15) + assert_near_equal(prob.get_val('generator.elec_power_out', units='kW'), np.ones(10)*90, tolerance=1e-15) + assert_near_equal(prob.get_val('generator.heat_out', units='kW'), np.ones(10)*0.0, tolerance=1e-15) + assert_near_equal(prob.get_val('generator.component_sizing_margin'), np.ones(10)*0.90, tolerance=1e-15) + assert_near_equal(prob.get_val('generator.component_cost', units='USD'), 13423.818791946309, tolerance=1e-10) + assert_near_equal(prob.get_val('generator.component_weight', units='kg'), 20, tolerance=1e-15) partials = prob.check_partials(method='cs',compact_print=True) assert_check_partials(partials) @@ -247,11 +247,11 @@ def test_nondefault_settings(self): )) prob.setup(check=True,force_alloc_complex=True) prob.run_model() - assert_rel_error(self, prob.get_val('generator.elec_power_out', units='kW'), np.ones(10)*90*0.95, tolerance=1e-15) - assert_rel_error(self, prob.get_val('generator.heat_out', units='kW'), np.ones(10)*90*0.05, tolerance=1e-15) - assert_rel_error(self, prob.get_val('generator.component_sizing_margin'), np.ones(10)*0.90*0.95, tolerance=1e-15) - assert_rel_error(self, prob.get_val('generator.component_cost', units='USD'), 203, tolerance=1e-15) - assert_rel_error(self, prob.get_val('generator.component_weight', units='kg'), 35.333333333333333333, tolerance=1e-10) + assert_near_equal(prob.get_val('generator.elec_power_out', units='kW'), np.ones(10)*90*0.95, tolerance=1e-15) + assert_near_equal(prob.get_val('generator.heat_out', units='kW'), np.ones(10)*90*0.05, tolerance=1e-15) + assert_near_equal(prob.get_val('generator.component_sizing_margin'), np.ones(10)*0.90*0.95, tolerance=1e-15) + assert_near_equal(prob.get_val('generator.component_cost', units='USD'), 203, tolerance=1e-15) + assert_near_equal(prob.get_val('generator.component_weight', units='kg'), 35.333333333333333333, tolerance=1e-10) partials = prob.check_partials(method='cs',compact_print=True) assert_check_partials(partials) @@ -262,11 +262,11 @@ def test_default_settings(self): use_defaults=True)) prob.setup(check=True,force_alloc_complex=True) prob.run_model() - assert_rel_error(self, prob.get_val('turboshaft.shaft_power_out', units='hp'), np.ones(10)*1000*0.9, tolerance=1e-6) - assert_rel_error(self, prob.get_val('turboshaft.fuel_flow', units='lbm/h'), np.ones(10)*0.6*1000*0.9, tolerance=1e-6) - assert_rel_error(self, prob.get_val('turboshaft.component_sizing_margin'), np.ones(10)*0.90, tolerance=1e-15) - assert_rel_error(self, prob.get_val('turboshaft.component_cost', units='USD'), 775528., tolerance=1e-7) - assert_rel_error(self, prob.get_val('turboshaft.component_weight', units='kg'), 0, tolerance=1e-15) + assert_near_equal(prob.get_val('turboshaft.shaft_power_out', units='hp'), np.ones(10)*1000*0.9, tolerance=1e-6) + assert_near_equal(prob.get_val('turboshaft.fuel_flow', units='lbm/h'), np.ones(10)*0.6*1000*0.9, tolerance=1e-6) + assert_near_equal(prob.get_val('turboshaft.component_sizing_margin'), np.ones(10)*0.90, tolerance=1e-15) + assert_near_equal(prob.get_val('turboshaft.component_cost', units='USD'), 775528., tolerance=1e-7) + assert_near_equal(prob.get_val('turboshaft.component_weight', units='kg'), 0, tolerance=1e-15) partials = prob.check_partials(method='cs',compact_print=True) assert_check_partials(partials) @@ -280,10 +280,10 @@ def test_nondefault_settings(self): cost_base=50000)) prob.setup(check=True,force_alloc_complex=True) prob.run_model() - assert_rel_error(self, prob.get_val('turboshaft.shaft_power_out', units='hp'), np.ones(10)*1000*0.9, tolerance=1e-6) - assert_rel_error(self, prob.get_val('turboshaft.fuel_flow', units='lbm/h'), np.ones(10)*0.4*1000*0.9, tolerance=1e-6) - assert_rel_error(self, prob.get_val('turboshaft.component_sizing_margin'), np.ones(10)*0.90, tolerance=1e-15) - assert_rel_error(self, prob.get_val('turboshaft.component_cost', units='USD'), 1e6+50000., tolerance=1e-6) - assert_rel_error(self, prob.get_val('turboshaft.component_weight', units='kg'), 1050, tolerance=1e-10) + assert_near_equal(prob.get_val('turboshaft.shaft_power_out', units='hp'), np.ones(10)*1000*0.9, tolerance=1e-6) + assert_near_equal(prob.get_val('turboshaft.fuel_flow', units='lbm/h'), np.ones(10)*0.4*1000*0.9, tolerance=1e-6) + assert_near_equal(prob.get_val('turboshaft.component_sizing_margin'), np.ones(10)*0.90, tolerance=1e-15) + assert_near_equal(prob.get_val('turboshaft.component_cost', units='USD'), 1e6+50000., tolerance=1e-6) + assert_near_equal(prob.get_val('turboshaft.component_weight', units='kg'), 1050, tolerance=1e-10) partials = prob.check_partials(method='cs',compact_print=True) assert_check_partials(partials) \ No newline at end of file diff --git a/openconcept/components/thermal.py b/openconcept/components/thermal.py index 97a4a3f0..ac6a853e 100644 --- a/openconcept/components/thermal.py +++ b/openconcept/components/thermal.py @@ -305,7 +305,7 @@ def setup(self): num_nodes=nn), promotes_inputs=['q_in', 'mass']) self.add_subsystem('integratetemp', - Integrator(num_intervals=int((nn-1)/2), + Integrator(num_nodes=nn, quantity_units='K', diff_units='s', method='simpson', @@ -362,7 +362,7 @@ def setup(self): CoolantReservoirRate(num_nodes=nn), promotes_inputs=['T_in', 'T_out', 'mass', 'mdot_coolant']) self.add_subsystem('integratetemp', - Integrator(num_intervals=int((nn-1)/2), + Integrator(num_nodes=nn, quantity_units='K', diff_units='s', method='simpson', diff --git a/openconcept/tests/test_example_aircraft.py b/openconcept/tests/test_example_aircraft.py index cf8b17bc..bd923946 100644 --- a/openconcept/tests/test_example_aircraft.py +++ b/openconcept/tests/test_example_aircraft.py @@ -2,7 +2,7 @@ from __future__ import division import unittest import numpy as np -from openmdao.utils.assert_utils import assert_rel_error, assert_check_partials +from openmdao.utils.assert_utils import assert_near_equal, assert_check_partials from openmdao.api import IndepVarComp, Group, Problem from examples.B738 import run_738_analysis from examples.TBM850 import run_tbm_analysis @@ -19,10 +19,10 @@ def setUp(self): def test_values_TBM(self): prob = self.prob - assert_rel_error(self, prob.get_val('climb.OEW', units='lb'), 4756.772140709275, tolerance=1e-5) - assert_rel_error(self, prob.get_val('rotate.range_final', units='ft'), 2489.7142498884746, tolerance=1e-5) - assert_rel_error(self, prob.get_val('engineoutclimb.gamma',units='deg'), 8.78263, tolerance=1e-5) - assert_rel_error(self, prob.get_val('descent.fuel_used_final', units='lb'), 1605.32123542, tolerance=1e-5) + assert_near_equal(prob.get_val('climb.OEW', units='lb'), 4756.772140709275, tolerance=1e-5) + assert_near_equal(prob.get_val('rotate.range_final', units='ft'), 2489.7142498884746, tolerance=1e-5) + assert_near_equal(prob.get_val('engineoutclimb.gamma',units='deg'), 8.78263, tolerance=1e-5) + assert_near_equal(prob.get_val('descent.fuel_used_final', units='lb'), 1605.32123542, tolerance=1e-5) class CaravanAnalysisTestCase(unittest.TestCase): def setUp(self): @@ -30,8 +30,8 @@ def setUp(self): def test_values_TBM(self): prob = self.prob - assert_rel_error(self, prob.get_val('v1vr.range_final', units='ft'), 1375.61684, tolerance=1e-5) - assert_rel_error(self, prob.get_val('descent.fuel_used_final', units='lb'), 377.4998448764594, tolerance=1e-5) + assert_near_equal(prob.get_val('v1vr.range_final', units='ft'), 1375.61684, tolerance=1e-5) + assert_near_equal(prob.get_val('descent.fuel_used_final', units='lb'), 377.4998448764594, tolerance=1e-5) class HybridTwinThermalTestCase(unittest.TestCase): def setUp(self): @@ -39,14 +39,14 @@ def setUp(self): def test_values_thermalhybridtwin(self): prob = self.prob - assert_rel_error(self, prob.get_val('climb.OEW', units='lb'), 6673.001027260613, tolerance=1e-5) - assert_rel_error(self, prob.get_val('rotate.range_final', units='ft'), 4434.461454141321, tolerance=1e-5) - assert_rel_error(self, prob.get_val('engineoutclimb.gamma',units='deg'), 1.7508055214855194, tolerance=1e-5) - assert_rel_error(self, prob.get_val('descent.fuel_used_final', units='lb'), 862.667748103923, tolerance=1e-5) + assert_near_equal(prob.get_val('climb.OEW', units='lb'), 6673.001027260613, tolerance=1e-5) + assert_near_equal(prob.get_val('rotate.range_final', units='ft'), 4434.461454141321, tolerance=1e-5) + assert_near_equal(prob.get_val('engineoutclimb.gamma',units='deg'), 1.7508055214855194, tolerance=1e-5) + assert_near_equal(prob.get_val('descent.fuel_used_final', units='lb'), 862.667748103923, tolerance=1e-5) - assert_rel_error(self, prob.get_val('climb.propmodel.motorheatsink.T', units='degC')[-1], 76.29252153, tolerance=1e-5) - assert_rel_error(self, prob.get_val('climb.propmodel.batteryheatsink.T', units='degC')[-1], -0.12704621031633678, tolerance=1e-5) - assert_rel_error(self, prob.get_val('cruise.propmodel.duct.drag', units='lbf')[0], 7.935834730687241, tolerance=1e-5) + assert_near_equal(prob.get_val('climb.propmodel.motorheatsink.T', units='degC')[-1], 76.29252153, tolerance=1e-5) + assert_near_equal(prob.get_val('climb.propmodel.batteryheatsink.T', units='degC')[-1], -0.12704621031633678, tolerance=1e-5) + assert_near_equal(prob.get_val('cruise.propmodel.duct.drag', units='lbf')[0], 7.935834730687241, tolerance=1e-5) class HybridTwinTestCase(unittest.TestCase): @@ -55,10 +55,10 @@ def setUp(self): def test_values_hybridtwin(self): prob = self.prob - assert_rel_error(self, prob.get_val('climb.OEW', units='lb'), 6648.424765080086, tolerance=1e-5) - assert_rel_error(self, prob.get_val('rotate.range_final', units='ft'), 4383.871458066499, tolerance=1e-5) - assert_rel_error(self, prob.get_val('engineoutclimb.gamma',units='deg'), 1.7659046316724112, tolerance=1e-5) - assert_rel_error(self, prob.get_val('descent.fuel_used_final', units='lb'), 854.8937776195904, tolerance=1e-5) + assert_near_equal(prob.get_val('climb.OEW', units='lb'), 6648.424765080086, tolerance=1e-5) + assert_near_equal(prob.get_val('rotate.range_final', units='ft'), 4383.871458066499, tolerance=1e-5) + assert_near_equal(prob.get_val('engineoutclimb.gamma',units='deg'), 1.7659046316724112, tolerance=1e-5) + assert_near_equal(prob.get_val('descent.fuel_used_final', units='lb'), 854.8937776195904, tolerance=1e-5) class KingAirTestCase(unittest.TestCase): def setUp(self): @@ -66,9 +66,9 @@ def setUp(self): def test_values_kingair(self): prob = self.prob - assert_rel_error(self, prob.get_val('climb.OEW', units='lb'), 6471.539115423346, tolerance=1e-5) - assert_rel_error(self, prob.get_val('rotate.range_final', units='ft'), 3056.9443135452075, tolerance=1e-5) - assert_rel_error(self, prob.get_val('descent.fuel_used_final', units='lb'), 1663.490303796347, tolerance=1e-5) + assert_near_equal(prob.get_val('climb.OEW', units='lb'), 6471.539115423346, tolerance=1e-5) + assert_near_equal(prob.get_val('rotate.range_final', units='ft'), 3056.9443135452075, tolerance=1e-5) + assert_near_equal(prob.get_val('descent.fuel_used_final', units='lb'), 1663.490303796347, tolerance=1e-5) class ElectricSingleTestCase(unittest.TestCase): @@ -77,9 +77,9 @@ def setUp(self): def test_values_kingair(self): prob = self.prob - assert_rel_error(self, prob.get_val('rotate.range_final', units='ft'), 2419.111568458725, tolerance=1e-5) - assert_rel_error(self, prob.get_val('descent.propmodel.batt1.SOC')[-1], 0.1663373102614198, tolerance=1e-5) - assert_rel_error(self, prob.get_val('descent.propmodel.motorheatsink.T', units='degC')[-1], 14.918329533221709, tolerance=1e-5) + assert_near_equal(prob.get_val('rotate.range_final', units='ft'), 2419.111568458725, tolerance=1e-5) + assert_near_equal(prob.get_val('descent.propmodel.batt1.SOC')[-1], 0.1663373102614198, tolerance=1e-5) + assert_near_equal(prob.get_val('descent.propmodel.motorheatsink.T', units='degC')[-1], 14.918329533221709, tolerance=1e-5) class B738TestCase(unittest.TestCase): def setUp(self): @@ -88,6 +88,6 @@ def setUp(self): def test_values_B738(self): prob = self.prob # block fuel - assert_rel_error(self, prob.get_val('descent.fuel_used_final', units='lbm'), 28688.32933661591, tolerance=2e-5) + assert_near_equal(prob.get_val('descent.fuel_used_final', units='lbm'), 28688.32933661591, tolerance=2e-5) # total fuel - assert_rel_error(self, prob.get_val('loiter.fuel_used_final', units='lbm'), 34555.31347454542, tolerance=2e-5) + assert_near_equal(prob.get_val('loiter.fuel_used_final', units='lbm'), 34555.31347454542, tolerance=2e-5) diff --git a/openconcept/utilities/math/add_subtract_comp.py b/openconcept/utilities/math/add_subtract_comp.py index af2d1b09..e9d8346c 100644 --- a/openconcept/utilities/math/add_subtract_comp.py +++ b/openconcept/utilities/math/add_subtract_comp.py @@ -1,10 +1,9 @@ """Definition of the Add/Subtract Component.""" -import collections import numpy as np from scipy import sparse as sp from six import string_types -from collections import Iterable +from collections.abc import Iterable from openmdao.core.explicitcomponent import ExplicitComponent @@ -81,7 +80,7 @@ def __init__(self, output_name=None, input_names=None, vec_size=1, length=1, if isinstance(output_name, string_types): self._add_systems.append((output_name, input_names, vec_size, length, val, scaling_factors, kwargs)) - elif isinstance(output_name, collections.Iterable): + elif isinstance(output_name, Iterable): raise NotImplementedError('Declaring multiple addition systems ' 'on initiation is not implemented.' 'Use a string to name a single addition relationship or use ' diff --git a/openconcept/utilities/math/combine_split_comp.py b/openconcept/utilities/math/combine_split_comp.py index fa044de9..48b08277 100644 --- a/openconcept/utilities/math/combine_split_comp.py +++ b/openconcept/utilities/math/combine_split_comp.py @@ -1,6 +1,6 @@ """Definition of the Vector Combiner/Splitter Component.""" -import collections +from collections.abc import Iterable import numpy as np from scipy import sparse as sp from six import string_types @@ -57,14 +57,14 @@ def __init__(self, output_name=None, input_names=None, vec_sizes=None, length=1, self._add_systems = [] if isinstance(output_name, string_types): - if (not isinstance(input_names, collections.Iterable) or - not isinstance(vec_sizes, collections.Iterable)): + if (not isinstance(input_names, Iterable) or + not isinstance(vec_sizes, Iterable)): raise ValueError('User must provide list of input name(s)' 'and list of vec_sizes for each input') self._add_systems.append((output_name, input_names, vec_sizes, length, val, kwargs)) - elif isinstance(output_name, collections.Iterable): + elif isinstance(output_name, Iterable): raise NotImplementedError('Declaring multiple relations ' 'on initiation is not implemented.' 'Use a string to name a single addition relationship or use ' @@ -140,8 +140,8 @@ def add_relation(self, output_name, input_names, vec_sizes, length=1, val=1.0, 'lower': lower, 'upper': upper, 'ref': ref, 'ref0': ref0, 'res_ref': res_ref} - if (not isinstance(input_names, collections.Iterable) or - not isinstance(vec_sizes, collections.Iterable)): + if (not isinstance(input_names, Iterable) or + not isinstance(vec_sizes, Iterable)): raise ValueError('User must provide list of input name(s)' 'and list of vec_sizes for each input') @@ -279,14 +279,14 @@ def __init__(self, output_names=None, input_name=None, vec_sizes=None, length=1, self._add_systems = [] if isinstance(input_name, string_types): - if (not isinstance(output_names, collections.Iterable) or - not isinstance(vec_sizes, collections.Iterable)): + if (not isinstance(output_names, Iterable) or + not isinstance(vec_sizes, Iterable)): raise ValueError('User must provide list of output name(s)' 'and list of vec_sizes for each input') self._add_systems.append((output_names, input_name, vec_sizes, length, val, kwargs)) - elif isinstance(input_name, collections.Iterable): + elif isinstance(input_name, Iterable): raise NotImplementedError('Declaring multiple relations ' 'on initiation is not implemented.' 'Use a string to name a single addition relationship or use ' @@ -362,8 +362,8 @@ def add_relation(self, output_names, input_name, vec_sizes, length=1, val=1.0, 'lower': lower, 'upper': upper, 'ref': ref, 'ref0': ref0, 'res_ref': res_ref} - if (not isinstance(output_names, collections.Iterable) or - not isinstance(vec_sizes, collections.Iterable)): + if (not isinstance(output_names, Iterable) or + not isinstance(vec_sizes, Iterable)): raise ValueError('User must provide list of output name(s)' 'and list of vec_sizes for each input') diff --git a/openconcept/utilities/math/integrals.py b/openconcept/utilities/math/integrals.py index 0d7876e6..31860fac 100644 --- a/openconcept/utilities/math/integrals.py +++ b/openconcept/utilities/math/integrals.py @@ -2,6 +2,7 @@ import numpy as np import scipy.sparse as sp from openmdao.api import ExplicitComponent +import warnings def bdf3_cache_matrix(n,all_bdf=False): """ @@ -472,22 +473,21 @@ def integrator_partials_wrt_deltas(num_segments, num_intervals): partials_q_wrt_deltas = sp.csr_matrix(jacmat) return partials_q_wrt_deltas -class IntegrateQuantityEveryNode(ExplicitComponent): +class NewIntegrator(ExplicitComponent): """ - This component integrates a vector using a 3rd order Lagrange polynomial - equivalent to Simpson's rule, but with quadrature between every subinterval + Integrates rate variables implicitly. + Add new integrated quantities by using the add_integrand method. + "q" inputs here are illustrative only. Inputs ------ - dqdt : float - The vector quantity to integrate. - Length of the vector = (2 * num_intervals + 1) * num_segments - segment|dt : float - The timestep of "segment" (scalar) - 1 per segment + duration : float + The duration of the integration interval (can also use dt) (scalar) + dq_dt : float + Rate to integrate (vector) q_initial : float - Starting value of the quantity (scalar) - + Starting value of quantity (scalar) + Outputs ------- q : float @@ -499,190 +499,289 @@ class IntegrateQuantityEveryNode(ExplicitComponent): Options ------- - segment_names : list - A list of str with the names of the individual segments - By default, if no segment_names are provided, one segment will be assumed and segment|dt will just be named "dt" - num_intervals : int - The number of Simpson intervals per segment - The total number of points per segment is 2N + 1 where N = num_intervals - The total length of the vector q is n_segments x (2N + 1) - quantity_units : str - The units of quantity being integrated (not the rate) + num_nodes : int + num_nodes = 2N + 1 where N = num_intervals + The total length of the vector q is 2N + 1 diff_units : str The units of the integrand (none by default) method : str - Numerical method (default 'simpson', optionally 'trap') - final_only : bool - Returns only the final value q_final and none of the full state. + Numerical method (default 'bdf3'; alternatively, 'simpson) + time_setup : str + Time configuration (default 'dt') + 'dt' creates input 'dt' + 'duration' creates input 'duration' + 'bounds' creates inputs 't_initial', 't_final' """ + def __init__(self, **kwargs): + super(NewIntegrator, self).__init__(**kwargs) + self._state_vars = {} + num_nodes = self.options['num_nodes'] + method = self.options['method'] + + # check to make sure num nodes is OK + if (num_nodes - 1) % 2 > 0: + raise ValueError('num_nodes must be odd') + + if num_nodes > 1: + if method == 'bdf3': + self.tri_mat, self.repeat_mat = bdf3_cache_matrix(num_nodes) + elif method == 'simpson': + self.tri_mat, self.repeat_mat = simpson_cache_matrix(num_nodes) def initialize(self): - self.options.declare('segment_names', default=None, desc="Names of differentiation segments") - self.options.declare('quantity_units',default=None, desc="Units of the quantity being differentiated") self.options.declare('diff_units',default=None, desc="Units of the differential") - self.options.declare('num_intervals',default=5, desc="Number of Simpsons rule intervals per segment") - self.options.declare('method',default='simpson', desc="Numerical method to use.") - self.options.declare('final_only',default=False) + self.options.declare('num_nodes',default=11, desc="Analysis points per segment") + self.options.declare('method',default='bdf3', desc="Numerical method to use.") + self.options.declare('time_setup',default='dt') + + def add_integrand(self, name, rate_name=None, start_name=None, end_name=None, + units=None, rate_units=None, zero_start=False, final_only=False, lower=-1e30, upper=1e30): + """ + Add a new integrated variable q = integrate(dqdt) + q0 + This will add an output with the integrated quantity, an output with the final value, + an input with the rate source, and an input for the initial quantity. + + Parameters + ---------- + name : str + The name of the integrated variable to be created. + rate_name : str + The name of the input rate (default name"_rate") + start_name : str + The name of the initial value input (default value name"_initial") + end_name : str + The name of the end value output (default value name"_final") + units : str or None + Units for the integrated quantity (or inferred automatically from rate_units) + rate_units : str or None + Units of the rate (can be inferred automatically from units) + zero_start : bool + If true, eliminates start value input and always begins from zero (default False) + final_only : bool + If true, only integrates final quantity, not all the intermediate points (default False) + upper : float + Upper bound on integrated quantity + lower : float + Lower bound on integrated quantity + """ + + num_nodes = self.options['num_nodes'] + diff_units = self.options['diff_units'] + time_setup = self.options['time_setup'] + # check to make sure num nodes is OK + if (num_nodes - 1) % 2 > 0: + raise ValueError('num_nodes must be odd') + + if units and rate_units: + raise ValueError('Specify either quantity units or rate units, but not both') + if units: + # infer rate units from diff units and quantity units + if not diff_units: + rate_units = units + warnings.warn('You have specified a integral with respect to a unitless integrand. Be aware of this.') + else: + rate_units = '('+units+') / (' + diff_units +')' + elif rate_units: + # infer quantity units from rate units and diff units + if not diff_units: + units = rate_units + warnings.warn('You have specified a integral with respect to a unitless integrand. Be aware of this.') + else: + units = '('+rate_units+') * (' + diff_units + ')' + elif diff_units: + # neither quantity nor rate units specified + rate_units = '(' + diff_units +')** -1' + + if not rate_name: + rate_name = name + '_rate' + if not start_name: + start_name = name + '_initial' + if not end_name: + end_name = name + '_final' + + options = {'name': name, + 'rate_name': rate_name, + 'start_name': start_name, + 'end_name': end_name, + 'units': units, + 'rate_units': rate_units, + 'zero_start': zero_start, + 'final_only': final_only, + 'upper': upper, + 'lower': lower} + + # TODO maybe later can pass kwargs + self._state_vars[name] = options + self.add_input(rate_name, val=0, shape=(num_nodes), units=rate_units) + self.add_output(end_name, units=units, upper=options['upper'],lower=options['lower']) + if not final_only: + self.add_output(name, shape=(num_nodes), units=units, upper=options['upper'],lower=options['lower']) + if not zero_start: + self.add_input(start_name, val=0, units=units) + if not final_only: + self.declare_partials([name], [start_name], rows=np.arange(num_nodes), cols=np.zeros((num_nodes,)), val=np.ones((num_nodes,))) + self.declare_partials([end_name], [start_name], val=1) + + # set up sparse partial structure + if num_nodes > 1: + # single point analysis has no dqdt dependency since the outputs are equal to the inputs + dQdrate, dQddtlist = multistep_integrator(0, np.ones((num_nodes,)), np.ones((1,)), self.tri_mat, self.repeat_mat, + segment_names=None, segments_to_count=None, partials=True) + dQdrate_indices = dQdrate.nonzero() + dQfdrate_indices = dQdrate.getrow(-1).nonzero() + if not final_only: + self.declare_partials([name], [rate_name], rows=dQdrate_indices[0], cols=dQdrate_indices[1]) + self.declare_partials([end_name], [rate_name], rows=dQfdrate_indices[0], cols=dQfdrate_indices[1]) # rows are zeros + + dQddt_seg = dQddtlist[0] + dQddt_indices = dQddt_seg.nonzero() + dQfddt_indices = dQddt_seg.getrow(-1).nonzero() + + if time_setup == 'dt': + if not final_only: + self.declare_partials([name], ['dt'], rows=dQddt_indices[0], cols=dQddt_indices[1]) + self.declare_partials([end_name], ['dt'], rows=dQfddt_indices[0], cols=dQfddt_indices[1]) + elif time_setup == 'duration': + if not final_only: + self.declare_partials([name], ['duration'], rows=dQddt_indices[0], cols=dQddt_indices[1]) + self.declare_partials([end_name], ['duration'], rows=dQfddt_indices[0], cols=dQfddt_indices[1]) + elif time_setup == 'bounds': + if not final_only: + self.declare_partials([name], ['t_initial','t_final'], rows=dQddt_indices[0], cols=dQddt_indices[1]) + self.declare_partials([end_name], ['t_initial','t_final'], rows=dQfddt_indices[0], cols=dQfddt_indices[1]) + else: + raise ValueError('Only dt, duration, and bounds are allowable values of time_setup') + def setup(self): - segment_names = self.options['segment_names'] - quantity_units = self.options['quantity_units'] diff_units = self.options['diff_units'] - n_int_per_seg = self.options['num_intervals'] + num_nodes = self.options['num_nodes'] method = self.options['method'] - nn_seg = (n_int_per_seg*2 + 1) - final_only = self.options['final_only'] + time_setup = self.options['time_setup'] - if segment_names is None: - n_segments = 1 - else: - n_segments = len(segment_names) - nn_tot = nn_seg * n_segments + # check to make sure num nodes is OK + if (num_nodes - 1) % 2 > 0: + raise ValueError('num_nodes must be odd') - if quantity_units is None and diff_units is None: - rate_units = None - elif quantity_units is None: - rate_units = '(' + diff_units +')** -1' - elif diff_units is None: - rate_units = quantity_units - warnings.warn('You have specified a integral with respect to a unitless integrand. Be aware of this.') + # branch logic here for the corner case of 0 segments + # so point analysis can be run without breaking everything + if num_nodes == 1: + single_point = True else: - rate_units = '('+quantity_units+') / (' + diff_units +')' - # the output of this function is of length nn - 1. NO partial for first row (initial value) - # get the partials of the delta quantities WRT the rates dDelta / drate - if method == "simpson": - delta_q, dDelta_ddqdt, dDelta_dts = three_point_lagrange_integration(np.ones((nn_tot,)), np.ones((n_segments, )), - num_segments=n_segments, num_intervals=n_int_per_seg) - elif method == "trap": - delta_q, dDelta_ddqdt, dDelta_dts = trapezoid_integration(np.ones((nn_tot,)), np.ones((n_segments, )), - num_segments=n_segments, num_intervals=n_int_per_seg) - elif method == "backward_euler": - delta_q, dDelta_ddqdt, dDelta_dts = backward_euler(np.ones((nn_tot,)), np.ones((n_segments, )), - num_segments=n_segments, num_intervals=n_int_per_seg) - dq_dDelta = integrator_partials_wrt_deltas(n_segments, n_int_per_seg) - # we need the partial of q with respect to the rates - # which is dq / dDelta * dDelta / d{parameter} - dq_ddqdt = dq_dDelta.dot(dDelta_ddqdt) + single_point = False + if not single_point: + if method == 'bdf3': + self.tri_mat, self.repeat_mat = bdf3_cache_matrix(num_nodes) + elif method == 'simpson': + self.tri_mat, self.repeat_mat = simpson_cache_matrix(num_nodes) - self.add_input('dqdt', val=0, units=rate_units, desc='Quantity to integrate',shape=(nn_tot,)) - self.add_input('q_initial', val=0, units=quantity_units, desc='Initial value') - if not final_only: - self.add_output('q', units=quantity_units, desc='Integral of dqdt', shape=(nn_tot,)) - self.declare_partials(['q'], ['q_initial'], rows=np.arange(nn_tot), cols=np.zeros((nn_tot,)), val=np.ones((nn_tot,))) + if time_setup == 'dt': + self.add_input('dt', units=diff_units, desc='Time step') + elif time_setup == 'duration': + self.add_input('duration', units=diff_units, desc='Time duration') + elif time_setup == 'bounds': + self.add_input('t_initial', units=diff_units, desc='Initial time') + self.add_input('t_final', units=diff_units, desc='Initial time') + else: + raise ValueError('Only dt, duration, and bounds are allowable values of time_setup') - self.add_output('q_final', units=quantity_units, desc='Final value of q') - self.declare_partials(['q_final'], ['q_initial'], val=1) - dq_ddqdt_indices = dq_ddqdt.nonzero() - dqfinal_ddqdt_indices = dq_ddqdt.getrow(-1).nonzero() - if not final_only: - self.declare_partials(['q'], ['dqdt'], rows=dq_ddqdt_indices[0], cols=dq_ddqdt_indices[1]) - self.declare_partials(['q_final'], ['dqdt'], rows=dqfinal_ddqdt_indices[0], cols=dqfinal_ddqdt_indices[1]) # rows are zeros + def compute(self, inputs, outputs): + num_nodes = self.options['num_nodes'] + time_setup=self.options['time_setup'] - if segment_names is None: - self.add_input('dt', units=diff_units, desc='Time step') - dq_ddt = dq_dDelta.dot(dDelta_dts[0]) - dq_ddt_indices = dq_ddt.nonzero() - if not final_only: - self.declare_partials(['q'], ['dt'], rows=dq_ddt_indices[0], cols=dq_ddt_indices[1]) - dqfinal_ddt_indices = dq_ddt.getrow(-1).nonzero() - self.declare_partials(['q_final'], ['dt'], rows=dqfinal_ddt_indices[0], cols=dqfinal_ddt_indices[1]) + if num_nodes == 1: + single_point = True else: - for i_seg, segment_name in enumerate(segment_names): - self.add_input(segment_name +'|dt', units=diff_units, desc='Time step') - dq_ddt = dq_dDelta.dot(dDelta_dts[i_seg]) - dq_ddt_indices = dq_ddt.nonzero() - if not final_only: - self.declare_partials(['q'], [segment_name +'|dt'], rows=dq_ddt_indices[0], cols=dq_ddt_indices[1]) - dqfinal_ddt_indices = dq_ddt.getrow(-1).nonzero() - self.declare_partials(['q_final'], [segment_name +'|dt'], rows=dqfinal_ddt_indices[0], cols=dqfinal_ddt_indices[1]) + single_point = False - def compute(self, inputs, outputs): - segment_names = self.options['segment_names'] - n_int_per_seg = self.options['num_intervals'] - method = self.options['method'] - final_only = self.options['final_only'] - nn_seg = (n_int_per_seg*2 + 1) - if segment_names is None: - n_segments = 1 + if time_setup == 'dt': dts = [inputs['dt'][0]] - else: - n_segments = len(segment_names) - dts = [] - for i_seg, segment_name in enumerate(segment_names): - input_name = segment_name+'|dt' - dts.append(inputs[input_name][0]) - if method == 'simpson': - delta_q, dDelta_ddqdt, dDelta_dts = three_point_lagrange_integration(inputs['dqdt'], dts, - num_segments=n_segments, num_intervals=n_int_per_seg) - elif method == 'trap': - delta_q, dDelta_ddqdt, dDelta_dts = trapezoid_integration(inputs['dqdt'], dts, - num_segments=n_segments, num_intervals=n_int_per_seg) - elif method == "backward_euler": - delta_q, dDelta_ddqdt, dDelta_dts = backward_euler(inputs['dqdt'], dts, - num_segments=n_segments, num_intervals=n_int_per_seg) - cumsum = np.cumsum(delta_q) - cumsum = np.insert(cumsum, 0, 0) - cumsum = cumsum + inputs['q_initial'] - if n_segments > 1: - for i in range(1,n_segments): - duplicate_row = cumsum[i*nn_seg-1] - cumsum = np.insert(cumsum,i*nn_seg,duplicate_row) - if not final_only: - outputs['q'] = cumsum - outputs['q_final'] = cumsum[-1] + elif time_setup == 'duration': + if num_nodes == 1: + dts = [inputs['duration'][0]] + else: + dts = [inputs['duration'][0]/(num_nodes-1)] + elif time_setup == 'bounds': + delta_t = inputs['t_final'] - inputs['t_initial'] + dts = [delta_t[0]/(num_nodes-1)] + + for name, options in self._state_vars.items(): + if options['zero_start']: + q0 = 0 + else: + q0 = inputs[options['start_name']] + if not single_point: + Q = multistep_integrator(q0, inputs[options['rate_name']], dts, self.tri_mat, self.repeat_mat, + segment_names=None, segments_to_count=None, partials=False) + else: + # single point case, no change, no dependence on time + Q = q0 - def compute_partials(self, inputs, J): - segment_names = self.options['segment_names'] - quantity_units = self.options['quantity_units'] - diff_units = self.options['diff_units'] - n_int_per_seg = self.options['num_intervals'] - method = self.options['method'] - final_only = self.options['final_only'] + if not options['final_only']: + outputs[options['name']] = Q + outputs[options['end_name']] = Q[-1] - nn_seg = (n_int_per_seg*2 + 1) - if segment_names is None: - n_segments = 1 - else: - n_segments = len(segment_names) - nn_tot = nn_seg * n_segments + def compute_partials(self, inputs, J): + num_nodes = self.options['num_nodes'] + time_setup = self.options['time_setup'] - if segment_names is None: - n_segments = 1 - dts = [inputs['dt'][0]] + if num_nodes == 1: + single_point = True else: - n_segments = len(segment_names) - dts = [] - for i_seg, segment_name in enumerate(segment_names): - input_name = segment_name+'|dt' - dts.append(inputs[input_name][0]) - if method == 'simpson': - delta_q, dDelta_ddqdt, dDelta_dts = three_point_lagrange_integration(inputs['dqdt'], dts, - num_segments=n_segments, num_intervals=n_int_per_seg) - elif method == 'trap': - delta_q, dDelta_ddqdt, dDelta_dts = trapezoid_integration(inputs['dqdt'], dts, - num_segments=n_segments, num_intervals=n_int_per_seg) - elif method == "backward_euler": - delta_q, dDelta_ddqdt, dDelta_dts = backward_euler(inputs['dqdt'], dts, - num_segments=n_segments, num_intervals=n_int_per_seg) - dq_dDelta = integrator_partials_wrt_deltas(n_segments, n_int_per_seg) - - dq_ddqdt = dq_dDelta.dot(dDelta_ddqdt) - if not final_only: - J['q','dqdt'] = dq_ddqdt.data - J['q_final', 'dqdt'] = dq_ddqdt.getrow(-1).data + single_point = False + if not single_point: + if time_setup == 'dt': + dts = [inputs['dt'][0]] + elif time_setup == 'duration': + dts = [inputs['duration'][0]/(num_nodes-1)] + elif time_setup == 'bounds': + delta_t = inputs['t_final'] - inputs['t_initial'] + dts = [delta_t[0]/(num_nodes-1)] + + for name, options in self._state_vars.items(): + start_name = options['start_name'] + end_name = options['end_name'] + qty_name = options['name'] + rate_name = options['rate_name'] + final_only = options['final_only'] + if options['zero_start']: + q0 = 0 + else: + q0 = inputs[start_name] + dQdrate, dQddtlist = multistep_integrator(q0, inputs[rate_name], dts, self.tri_mat, self.repeat_mat, + segment_names=None, segments_to_count=None, partials=True) - if segment_names is None: - dq_ddt = dq_dDelta.dot(dDelta_dts[0]) - if not final_only: - J['q','dt'] = dq_ddt.data - J['q_final','dt'] = dq_ddt.getrow(-1).data - else: - for i_seg, segment_name in enumerate(segment_names): - dq_ddt = dq_dDelta.dot(dDelta_dts[i_seg]) if not final_only: - J['q',segment_name+'|dt'] = dq_ddt.data - J['q_final',segment_name+'|dt'] = dq_ddt.getrow(-1).data + J[qty_name, rate_name] = dQdrate.data + J[end_name, rate_name] = dQdrate.getrow(-1).data + + if time_setup == 'dt': + if not final_only: + J[qty_name, 'dt'] = np.squeeze(dQddtlist[0].toarray()[1:]) + J[end_name, 'dt'] = np.squeeze(dQddtlist[0].getrow(-1).toarray()) + + elif time_setup == 'duration': + if not final_only: + J[qty_name, 'duration'] = np.squeeze(dQddtlist[0].toarray()[1:] / (num_nodes - 1)) + J[end_name, 'duration'] = np.squeeze(dQddtlist[0].getrow(-1).toarray() / (num_nodes - 1)) + + elif time_setup == 'bounds': + if not final_only: + if len(dQddtlist[0].data) == 0: + J[qty_name, 't_initial'] = np.zeros(J[qty_name, 't_initial'].shape) + J[qty_name, 't_final'] = np.zeros(J[qty_name, 't_final'].shape) + else: + J[qty_name, 't_initial'] = -dQddtlist[0].data / (num_nodes - 1) + J[qty_name, 't_final'] = dQddtlist[0].data / (num_nodes - 1) + if len(dQddtlist[0].getrow(-1).data) == 0: + J[end_name, 't_initial'] = 0 + J[end_name, 't_final'] = 0 + else: + J[end_name, 't_initial'] = -dQddtlist[0].getrow(-1).data / (num_nodes - 1) + J[end_name, 't_final'] = dQddtlist[0].getrow(-1).data / (num_nodes - 1) + + class Integrator(ExplicitComponent): """ @@ -717,14 +816,15 @@ class Integrator(ExplicitComponent): segments_to_count : list A list of str with the names of segments to be included in the integration. By default, ALL segments will be included. - num_intervals : int - The number of Simpson intervals per segment - The total number of points per segment is 2N + 1 where N = num_intervals + num_nodes : int + num_nodes = 2N + 1 where N = num_intervals The total length of the vector q is n_segments x (2N + 1) quantity_units : str The units of quantity being integrated (not the rate) diff_units : str The units of the integrand (none by default) + rate_units : str + The units of the rate being integrated method : str Numerical method (default 'bdf3'; alternatively, 'simpson) zero_start : bool @@ -743,7 +843,8 @@ def initialize(self): self.options.declare('segments_to_count', default=None, desc="Names of differentiation segments") self.options.declare('quantity_units',default=None, desc="Units of the quantity being differentiated") self.options.declare('diff_units',default=None, desc="Units of the differential") - self.options.declare('num_intervals',default=5, desc="Number of Simpsons rule intervals per segment") + self.options.declare('rate_units',default=None, desc="Units of the rate being integrated") + self.options.declare('num_nodes',default=11, desc="Analysis points per segment") self.options.declare('method',default='bdf3', desc="Numerical method to use.") self.options.declare('zero_start',default=False) self.options.declare('final_only',default=False) @@ -756,30 +857,35 @@ def setup(self): segments_to_count = self.options['segments_to_count'] quantity_units = self.options['quantity_units'] diff_units = self.options['diff_units'] - n_int_per_seg = self.options['num_intervals'] + num_nodes = self.options['num_nodes'] method = self.options['method'] zero_start = self.options['zero_start'] final_only = self.options['final_only'] time_setup = self.options['time_setup'] - nn_seg = (n_int_per_seg*2 + 1) + + # check to make sure num nodes is OK + if (num_nodes - 1) % 2 > 0: + raise ValueError('num_nodes must be odd') # branch logic here for the corner case of 0 segments # so point analysis can be run without breaking everything - if nn_seg == 1: + if num_nodes == 1: single_point = True else: single_point = False if not single_point: if method == 'bdf3': - self.tri_mat, self.repeat_mat = bdf3_cache_matrix(nn_seg) + self.tri_mat, self.repeat_mat = bdf3_cache_matrix(num_nodes) elif method == 'simpson': - self.tri_mat, self.repeat_mat = simpson_cache_matrix(nn_seg) + self.tri_mat, self.repeat_mat = simpson_cache_matrix(num_nodes) if segment_names is None: n_segments = 1 else: n_segments = len(segment_names) - nn_tot = nn_seg * n_segments + nn_tot = num_nodes * n_segments + + # TODO enable specifying rate units if quantity_units is None and diff_units is None: rate_units = None @@ -863,14 +969,13 @@ def setup(self): def compute(self, inputs, outputs): segment_names = self.options['segment_names'] - n_int_per_seg = self.options['num_intervals'] + num_nodes = self.options['num_nodes'] segments_to_count = self.options['segments_to_count'] zero_start = self.options['zero_start'] final_only = self.options['final_only'] time_setup=self.options['time_setup'] - nn_seg = (n_int_per_seg*2 + 1) - if nn_seg == 1: + if num_nodes == 1: single_point = True else: single_point = False @@ -880,13 +985,13 @@ def compute(self, inputs, outputs): if time_setup == 'dt': dts = [inputs['dt'][0]] elif time_setup == 'duration': - if nn_seg == 1: + if num_nodes == 1: dts = [inputs['duration'][0]] else: - dts = [inputs['duration'][0]/(nn_seg-1)] + dts = [inputs['duration'][0]/(num_nodes-1)] elif time_setup == 'bounds': delta_t = inputs['t_final'] - inputs['t_initial'] - dts = [delta_t[0]/(nn_seg-1)] + dts = [delta_t[0]/(num_nodes-1)] else: n_segments = len(segment_names) dts = [] @@ -915,14 +1020,13 @@ def compute_partials(self, inputs, J): segment_names = self.options['segment_names'] quantity_units = self.options['quantity_units'] diff_units = self.options['diff_units'] - n_int_per_seg = self.options['num_intervals'] + num_nodes = self.options['num_nodes'] segments_to_count = self.options['segments_to_count'] zero_start = self.options['zero_start'] final_only = self.options['final_only'] time_setup = self.options['time_setup'] - nn_seg = (n_int_per_seg*2 + 1) - if nn_seg == 1: + if num_nodes == 1: single_point = True else: single_point = False @@ -931,17 +1035,17 @@ def compute_partials(self, inputs, J): n_segments = 1 else: n_segments = len(segment_names) - nn_tot = nn_seg * n_segments + nn_tot = num_nodes * n_segments if segment_names is None: n_segments = 1 if time_setup == 'dt': dts = [inputs['dt'][0]] elif time_setup == 'duration': - dts = [inputs['duration'][0]/(nn_seg-1)] + dts = [inputs['duration'][0]/(num_nodes-1)] elif time_setup == 'bounds': delta_t = inputs['t_final'] - inputs['t_initial'] - dts = [delta_t[0]/(nn_seg-1)] + dts = [delta_t[0]/(num_nodes-1)] else: n_segments = len(segment_names) dts = [] @@ -979,13 +1083,13 @@ def compute_partials(self, inputs, J): # if len(dQddtlist[0].data) == 0: # J['q','duration'] = np.zeros(J['q','duration'].shape) # else: - # J['q','duration'] = dQddtlist[0].data / (nn_seg - 1) - J['q','duration'] = np.squeeze(dQddtlist[0].toarray()[1:] / (nn_seg - 1)) + # J['q','duration'] = dQddtlist[0].data / (num_nodes - 1) + J['q','duration'] = np.squeeze(dQddtlist[0].toarray()[1:] / (num_nodes - 1)) # if len(dQddtlist[0].getrow(-1).data) == 0: # J['q_final','duration'] = 0 # else: - # J['q_final','duration'] = dQddtlist[0].getrow(-1).data / (nn_seg - 1) - J['q_final','duration'] = np.squeeze(dQddtlist[0].getrow(-1).toarray() / (nn_seg - 1)) + # J['q_final','duration'] = dQddtlist[0].getrow(-1).data / (num_nodes - 1) + J['q_final','duration'] = np.squeeze(dQddtlist[0].getrow(-1).toarray() / (num_nodes - 1)) elif time_setup == 'bounds': if not final_only: @@ -993,14 +1097,14 @@ def compute_partials(self, inputs, J): J['q','t_initial'] = np.zeros(J['q','t_initial'].shape) J['q','t_final'] = np.zeros(J['q','t_final'].shape) else: - J['q','t_initial'] = -dQddtlist[0].data / (nn_seg - 1) - J['q','t_final'] = dQddtlist[0].data / (nn_seg - 1) + J['q','t_initial'] = -dQddtlist[0].data / (num_nodes - 1) + J['q','t_final'] = dQddtlist[0].data / (num_nodes - 1) if len(dQddtlist[0].getrow(-1).data) == 0: J['q_final','t_initial'] = 0 J['q_final','t_final'] = 0 else: - J['q_final','t_initial'] = -dQddtlist[0].getrow(-1).data / (nn_seg - 1) - J['q_final','t_final'] = dQddtlist[0].getrow(-1).data / (nn_seg - 1) + J['q_final','t_initial'] = -dQddtlist[0].getrow(-1).data / (num_nodes - 1) + J['q_final','t_final'] = dQddtlist[0].getrow(-1).data / (num_nodes - 1) else: for i_seg, segment_name in enumerate(segment_names): if not final_only: diff --git a/openconcept/utilities/math/multiply_divide_comp.py b/openconcept/utilities/math/multiply_divide_comp.py index a267fd42..e848b6d5 100644 --- a/openconcept/utilities/math/multiply_divide_comp.py +++ b/openconcept/utilities/math/multiply_divide_comp.py @@ -1,10 +1,9 @@ """Definition of the Element Multiply Component.""" -import collections import numpy as np from scipy import sparse as sp from six import string_types -from collections import Iterable +from collections.abc import Iterable from openmdao.core.explicitcomponent import ExplicitComponent @@ -84,7 +83,7 @@ def __init__(self, output_name=None, input_names=None, vec_size=1, length=1, if isinstance(output_name, string_types): self._add_systems.append((output_name, input_names, vec_size, length, val, scaling_factor, divide, input_units, kwargs)) - elif isinstance(output_name, collections.Iterable): + elif isinstance(output_name, Iterable): raise NotImplementedError('Declaring multiple systems ' 'on initiation is not implemented.' 'Use a string to name a single addition relationship or use ' diff --git a/openconcept/utilities/math/sum_comp.py b/openconcept/utilities/math/sum_comp.py index 1a3ee0ef..7159d408 100644 --- a/openconcept/utilities/math/sum_comp.py +++ b/openconcept/utilities/math/sum_comp.py @@ -1,6 +1,6 @@ """Definition of the Element Summation Component.""" -import collections +from collections.abc import Iterable import numpy as np from scipy import sparse as sp from six import string_types @@ -72,7 +72,7 @@ def __init__(self, output_name=None, input_name=None, vec_size=1, length=1, if isinstance(output_name, string_types): self._add_systems.append((output_name, input_name, vec_size, length, val, scaling_factor, kwargs)) - elif isinstance(output_name, collections.Iterable): + elif isinstance(output_name, Iterable): raise NotImplementedError('Declaring multiple systems ' 'on initiation is not implemented.' 'Use a string to name a single addition relationship or use ' diff --git a/openconcept/utilities/math/tests/test_add_subtract_comp.py b/openconcept/utilities/math/tests/test_add_subtract_comp.py index 237a7ceb..b784c4e4 100644 --- a/openconcept/utilities/math/tests/test_add_subtract_comp.py +++ b/openconcept/utilities/math/tests/test_add_subtract_comp.py @@ -7,7 +7,7 @@ from openmdao.api import Problem, Group, IndepVarComp #from openmdao.components.add_subtract_comp import AddSubtractComp from openconcept.utilities.math.add_subtract_comp import AddSubtractComp -from openmdao.utils.assert_utils import assert_rel_error, assert_check_partials +from openmdao.utils.assert_utils import assert_near_equal, assert_check_partials class TestAddSubtractCompScalars(unittest.TestCase): @@ -41,7 +41,7 @@ def test_results(self): b = self.p['b'] out = self.p['add_subtract_comp.adder_output'] expected = a + b - assert_rel_error(self, out, expected,1e-16) + assert_near_equal(out, expected,1e-16) def test_partials(self): partials = self.p.check_partials(method='fd', out_stream=None) @@ -81,7 +81,7 @@ def test_results(self): b = self.p['b'] out = self.p['add_subtract_comp.adder_output'] expected = a + b - assert_rel_error(self, out, expected,1e-16) + assert_near_equal(out, expected,1e-16) def test_partials(self): partials = self.p.check_partials(method='fd', out_stream=None) @@ -120,7 +120,7 @@ def test_results(self): b = self.p['b'] out = self.p['add_subtract_comp.adder_output'] expected = a + b - assert_rel_error(self, out, expected,1e-16) + assert_near_equal(out, expected,1e-16) def test_partials(self): partials = self.p.check_partials(method='fd', out_stream=None) @@ -160,7 +160,7 @@ def test_results(self): b = self.p['b'] out = self.p['add_subtract_comp.adder_output'] expected = a + b - assert_rel_error(self, out, expected,1e-16) + assert_near_equal(out, expected,1e-16) def test_partials(self): partials = self.p.check_partials(method='fd', out_stream=None) @@ -204,7 +204,7 @@ def test_results(self): c = self.p['c'] out = self.p['add_subtract_comp.adder_output'] expected = a + b + c - assert_rel_error(self, out, expected,1e-16) + assert_near_equal(out, expected,1e-16) def test_partials(self): partials = self.p.check_partials(method='fd', out_stream=None) @@ -248,7 +248,7 @@ def test_results(self): c = self.p['c'] out = self.p['add_subtract_comp.adder_output'] expected = 2*a + b - c - assert_rel_error(self, out, expected,1e-16) + assert_near_equal(out, expected,1e-16) def test_partials(self): partials = self.p.check_partials(method='fd', out_stream=None) @@ -293,7 +293,7 @@ def test_results(self): out = self.p['add_subtract_comp.adder_output'] m_to_ft = 3.280839895 expected = a + b*m_to_ft + c*m_to_ft - assert_rel_error(self, out, expected,1e-8) + assert_near_equal(out, expected,1e-8) def test_partials(self): partials = self.p.check_partials(method='fd', out_stream=None) @@ -335,7 +335,7 @@ def test(self): import numpy as np #from openmdao.api import Problem, Group, IndepVarComp from openconcept.utilities.math.add_subtract_comp import AddSubtractComp - from openmdao.utils.assert_utils import assert_rel_error + from openmdao.utils.assert_utils import assert_near_equal n = 3 @@ -377,7 +377,7 @@ def test(self): # Verify the results expected_i = np.array([[100, 200, 300], [0, -1, -2]]).T - assert_rel_error(self, p.get_val('totalforcecomp.total_force', units='kN'), expected_i) + assert_near_equal(p.get_val('totalforcecomp.total_force', units='kN'), expected_i) if __name__ == '__main__': diff --git a/openconcept/utilities/math/tests/test_combine_split.py b/openconcept/utilities/math/tests/test_combine_split.py index 1590cd06..55954124 100644 --- a/openconcept/utilities/math/tests/test_combine_split.py +++ b/openconcept/utilities/math/tests/test_combine_split.py @@ -6,7 +6,7 @@ from openmdao.api import Problem, Group, IndepVarComp from openconcept.utilities.math.combine_split_comp import VectorConcatenateComp, VectorSplitComp -from openmdao.utils.assert_utils import assert_rel_error, assert_check_partials +from openmdao.utils.assert_utils import assert_near_equal, assert_check_partials class TestConcatenateScalars(unittest.TestCase): @@ -40,7 +40,7 @@ def test_results(self): b = self.p['b'] out = self.p['vector_concat_comp.concat_output'] expected = np.concatenate((a,b)) - assert_rel_error(self, out, expected,1e-16) + assert_near_equal(out, expected,1e-16) def test_partials(self): partials = self.p.check_partials(method='fd', out_stream=None) @@ -77,7 +77,7 @@ def test_results(self): b = self.p['b'] out = self.p['vector_concat_comp.concat_output'] expected = np.concatenate((a,b)) - assert_rel_error(self, out, expected,1e-16) + assert_near_equal(out, expected,1e-16) def test_partials(self): partials = self.p.check_partials(method='fd', out_stream=None) @@ -115,7 +115,7 @@ def test_results(self): b = self.p['b'] out = self.p['vector_concat_comp.concat_output'] expected = np.concatenate((a,b)) - assert_rel_error(self, out, expected,1e-16) + assert_near_equal(out, expected,1e-16) def test_partials(self): partials = self.p.check_partials(method='fd', out_stream=None) @@ -151,7 +151,7 @@ def test_results(self): b = self.p['b'] out = self.p['vector_concat_comp.concat_output'] expected = np.concatenate((a,b)) - assert_rel_error(self, out, expected,1e-16) + assert_near_equal(out, expected,1e-16) def test_partials(self): partials = self.p.check_partials(method='fd', out_stream=None) @@ -204,8 +204,8 @@ def test_results(self): expected1 = np.concatenate((a,b)) expected2 = np.concatenate((c,d)) - assert_rel_error(self, out1, expected1,1e-16) - assert_rel_error(self, out2, expected2,1e-16) + assert_near_equal(out1, expected1,1e-16) + assert_near_equal(out2, expected2,1e-16) def test_partials(self): partials = self.p.check_partials(method='fd', out_stream=None) @@ -244,7 +244,7 @@ def test_results(self): b = b*1000. out = self.p['vector_concat_comp.concat_output'] expected = np.concatenate((a,b)) - assert_rel_error(self, out, expected,1e-16) + assert_near_equal(out, expected,1e-16) class TestConcatenate3InputsDiffSizesNx1(unittest.TestCase): def setUp(self): @@ -283,7 +283,7 @@ def test_results(self): out = self.p['vector_concat_comp.concat_output'] expected = np.concatenate((a,b,c)) - assert_rel_error(self, out, expected,1e-16) + assert_near_equal(out, expected,1e-16) def test_partials(self): partials = self.p.check_partials(method='fd', out_stream=None) @@ -327,7 +327,7 @@ def test_results(self): out = self.p['vector_concat_comp.concat_output'] expected = np.concatenate((a,b,c)) - assert_rel_error(self, out, expected,1e-16) + assert_near_equal(out, expected,1e-16) def test_partials(self): partials = self.p.check_partials(method='fd', out_stream=None) @@ -385,8 +385,8 @@ def test_results(self): expected_a = input_to_split[0] expected_b = input_to_split[1] - assert_rel_error(self, out_a, expected_a,1e-16) - assert_rel_error(self, out_b, expected_b,1e-16) + assert_near_equal(out_a, expected_a,1e-16) + assert_near_equal(out_b, expected_b,1e-16) def test_partials(self): partials = self.p.check_partials(method='fd', out_stream=None) @@ -420,8 +420,8 @@ def test_results(self): expected_a = input_to_split[0:self.nn] expected_b = input_to_split[self.nn:2*self.nn] - assert_rel_error(self, out_a, expected_a,1e-16) - assert_rel_error(self, out_b, expected_b,1e-16) + assert_near_equal(out_a, expected_a,1e-16) + assert_near_equal(out_b, expected_b,1e-16) def test_partials(self): partials = self.p.check_partials(method='fd', out_stream=None) @@ -455,8 +455,8 @@ def test_results(self): expected_a = input_to_split[0:self.nn,:] expected_b = input_to_split[self.nn:2*self.nn,:] - assert_rel_error(self, out_a, expected_a,1e-16) - assert_rel_error(self, out_b, expected_b,1e-16) + assert_near_equal(out_a, expected_a,1e-16) + assert_near_equal(out_b, expected_b,1e-16) def test_partials(self): partials = self.p.check_partials(method='fd', out_stream=None) @@ -488,8 +488,8 @@ def test_results(self): expected_a = input_to_split[0:self.nn,:] expected_b = input_to_split[self.nn:2*self.nn,:] - assert_rel_error(self, out_a, expected_a,1e-16) - assert_rel_error(self, out_b, expected_b,1e-16) + assert_near_equal(out_a, expected_a,1e-16) + assert_near_equal(out_b, expected_b,1e-16) def test_partials(self): partials = self.p.check_partials(method='fd', out_stream=None) @@ -526,9 +526,9 @@ def test_results(self): expected_b = input_to_split[self.nn:2*self.nn,:] expected_c = input_to_split[2*self.nn:2*self.nn+2,:] - assert_rel_error(self, out_a, expected_a,1e-16) - assert_rel_error(self, out_b, expected_b,1e-16) - assert_rel_error(self, out_c, expected_c,1e-16) + assert_near_equal(out_a, expected_a,1e-16) + assert_near_equal(out_b, expected_b,1e-16) + assert_near_equal(out_c, expected_c,1e-16) def test_partials(self): partials = self.p.check_partials(method='fd', out_stream=None) @@ -562,8 +562,8 @@ def test_results(self): expected_a = input_to_split[0:self.nn,:] expected_b = input_to_split[self.nn:2*self.nn,:] - assert_rel_error(self, out_a, expected_a,1e-16) - assert_rel_error(self, out_b, expected_b,1e-16) + assert_near_equal(out_a, expected_a,1e-16) + assert_near_equal(out_b, expected_b,1e-16) def test_partials(self): partials = self.p.check_partials(method='fd', out_stream=None) @@ -605,7 +605,7 @@ def test(self): """ import numpy as np from openmdao.api import Problem, Group, IndepVarComp - from openmdao.utils.assert_utils import assert_rel_error + from openmdao.utils.assert_utils import assert_near_equal from openconcept.utilities.math.combine_split_comp import VectorConcatenateComp, VectorSplitComp n_takeoff_pts = 3 @@ -662,9 +662,9 @@ def test(self): expected_vel = np.array([[30, 0], [40, 0], [50, 0], [60, 5], [60, 0], [60, 0], [60, 0], [60, -5]]) expected_alt = np.array([0, 0, 0, 6000, 7500, 8000, 8500, 5000]) expected_split_vel = np.array([[60, 5], [60, 0], [60, 0], [60, 0], [60, -5]]) - assert_rel_error(self, p.get_val('combiner.velocity', units='m/s'), expected_vel) - assert_rel_error(self, p.get_val('combiner.altitude', units='m'), expected_alt) - assert_rel_error(self, p.get_val('divider.cruise_vel', units='m/s'), expected_split_vel) + assert_near_equal(p.get_val('combiner.velocity', units='m/s'), expected_vel) + assert_near_equal(p.get_val('combiner.altitude', units='m'), expected_alt) + assert_near_equal(p.get_val('divider.cruise_vel', units='m/s'), expected_split_vel) if __name__ == '__main__': unittest.main() \ No newline at end of file diff --git a/openconcept/utilities/math/tests/test_integrals.py b/openconcept/utilities/math/tests/test_integrals.py index 98ad132d..a4de7c59 100644 --- a/openconcept/utilities/math/tests/test_integrals.py +++ b/openconcept/utilities/math/tests/test_integrals.py @@ -2,7 +2,7 @@ import unittest import pytest import numpy as np -from openmdao.utils.assert_utils import assert_rel_error, assert_check_partials +from openmdao.utils.assert_utils import assert_near_equal, assert_check_partials from openmdao.api import IndepVarComp, Group, Problem from openconcept.utilities.math.integrals import Integrator @@ -14,10 +14,8 @@ class MultiPhaseIntegratorTestGroup(Group): segment_names : list A list of str with the names of the individual segments By default, if no segment_names are provided, one segment will be assumed and segment|dt will just be named "dt" - num_intervals : int - The number of Simpson intervals per segment - The total number of points per segment is 2N + 1 where N = num_intervals - The total length of the vector q is n_segments * (2N + 1) + num_nodes : int + Number of analysis points per segment quantity_units : str The units of the integral quantity q (NOT the rate) diff_units : str @@ -31,7 +29,7 @@ def initialize(self): self.options.declare('segments_to_count', default=None, desc="Names of differentiation segments") self.options.declare('quantity_units',default=None, desc="Units of the quantity being differentiated") self.options.declare('diff_units',default=None, desc="Units of the differential") - self.options.declare('num_intervals',default=5, desc="Number of Simpsons rule intervals per segment") + self.options.declare('num_nodes',default=11, desc="Number of nodes per segment") self.options.declare('integrator',default='simpson', desc="Which simpson integrator to use") self.options.declare('time_setup',default='dt') @@ -40,18 +38,18 @@ def setup(self): segments_to_count = self.options['segments_to_count'] quantity_units = self.options['quantity_units'] diff_units = self.options['diff_units'] - n_int_per_seg = self.options['num_intervals'] + num_nodes = self.options['num_nodes'] integrator_option = self.options['integrator'] time_setup = self.options['time_setup'] if segment_names is None: - nn_tot = (2*n_int_per_seg + 1) + nn_tot = num_nodes else: - nn_tot = (2*n_int_per_seg + 1) * len(segment_names) + nn_tot = num_nodes * len(segment_names) iv = self.add_subsystem('iv', IndepVarComp()) self.add_subsystem('integral', Integrator(segment_names=segment_names, segments_to_count=segments_to_count, quantity_units=quantity_units, - diff_units=diff_units, num_intervals=n_int_per_seg, method=integrator_option, time_setup=time_setup)) + diff_units=diff_units, num_nodes=num_nodes, method=integrator_option, time_setup=time_setup)) if quantity_units is None and diff_units is None: rate_units = None elif quantity_units is None: @@ -76,7 +74,7 @@ def setup(self): self.connect('iv.duration', 'integral.duration') elif time_setup == 'bounds': iv.add_output('t_initial', val=2, units=diff_units) - iv.add_output('t_final', val=2 + 1*(nn_tot-1)) + iv.add_output('t_final', val=2 + 1*(nn_tot-1), units=diff_units) self.connect('iv.t_initial','integral.t_initial') self.connect('iv.t_final','integral.t_final') else: @@ -90,120 +88,120 @@ class IntegratorEveryNodeCommonTestCases(object): """ def test_uniform_single_phase_no_units(self): - prob = Problem(MultiPhaseIntegratorTestGroup(num_intervals=self.num_intervals, integrator=self.integrator)) + prob = Problem(MultiPhaseIntegratorTestGroup(num_nodes=self.num_nodes, integrator=self.integrator)) prob.setup(check=True, force_alloc_complex=True) prob.run_model() - n_int_per_seg = self.num_intervals - nn_tot = (n_int_per_seg*2 + 1) - assert_rel_error(self, prob['integral.q'], np.linspace(0, nn_tot-1, nn_tot), tolerance=1e-14) - assert_rel_error(self, prob.get_val('integral.q_final', units=None), nn_tot-1, tolerance=1e-14) + num_nodes = self.num_nodes + nn_tot = num_nodes + assert_near_equal(prob['integral.q'], np.linspace(0, nn_tot-1, nn_tot), tolerance=1e-14) + assert_near_equal(prob.get_val('integral.q_final', units=None), nn_tot-1, tolerance=1e-14) partials = prob.check_partials(method='cs',compact_print=True) assert_check_partials(partials, atol=1e-8, rtol=1e0) def test_linear_single_phase_no_units(self): - n_int_per_seg = self.num_intervals - nn_tot = (n_int_per_seg*2 + 1) + num_nodes = self.num_nodes + nn_tot = num_nodes x = np.linspace(0, nn_tot-1, nn_tot) fprime = x f = x ** 2 / 2 - prob = Problem(MultiPhaseIntegratorTestGroup(num_intervals=self.num_intervals, integrator=self.integrator)) + prob = Problem(MultiPhaseIntegratorTestGroup(num_nodes=self.num_nodes, integrator=self.integrator)) prob.setup(check=True, force_alloc_complex=True) prob['iv.rate_to_integrate'] = fprime prob.run_model() - assert_rel_error(self, prob['integral.q'], f, tolerance=1e-14) - assert_rel_error(self, prob.get_val('integral.q_final', units=None), f[-1], tolerance=1e-14) + assert_near_equal(prob['integral.q'], f, tolerance=1e-14) + assert_near_equal(prob.get_val('integral.q_final', units=None), f[-1], tolerance=1e-14) partials = prob.check_partials(method='cs',compact_print=True) assert_check_partials(partials, atol=1e-8, rtol=1e0) def test_quadratic_single_phase_no_units(self): - n_int_per_seg = self.num_intervals - nn_tot = (n_int_per_seg*2 + 1) + num_nodes = self.num_nodes + nn_tot = num_nodes x = np.linspace(0, nn_tot-1, nn_tot) fprime = 4 * x **2 - 8*x + 5 f = 4 * x ** 3 / 3 - 8 * x ** 2 / 2 + 5*x - prob = Problem(MultiPhaseIntegratorTestGroup(num_intervals=self.num_intervals, integrator=self.integrator)) + prob = Problem(MultiPhaseIntegratorTestGroup(num_nodes=self.num_nodes, integrator=self.integrator)) prob.setup(check=True, force_alloc_complex=True) prob['iv.rate_to_integrate'] = fprime prob.run_model() - assert_rel_error(self, prob['integral.q'], f, tolerance=1e-14) - assert_rel_error(self, prob.get_val('integral.q_final', units=None), f[-1], tolerance=1e-14) + assert_near_equal(prob['integral.q'], f, tolerance=1e-14) + assert_near_equal(prob.get_val('integral.q_final', units=None), f[-1], tolerance=1e-14) partials = prob.check_partials(method='cs',compact_print=True) assert_check_partials(partials, atol=1e-8, rtol=1e0) def test_quadratic_single_phase_units(self): - n_int_per_seg = self.num_intervals - nn_tot = (n_int_per_seg*2 + 1) + num_nodes = self.num_nodes + nn_tot = num_nodes x = np.linspace(0, nn_tot-1, nn_tot) fprime = 4 * x **2 - 8*x + 5 f = 4 * x ** 3 / 3 - 8 * x ** 2 / 2 + 5*x - prob = Problem(MultiPhaseIntegratorTestGroup(num_intervals=self.num_intervals, integrator=self.integrator, + prob = Problem(MultiPhaseIntegratorTestGroup(num_nodes=self.num_nodes, integrator=self.integrator, quantity_units='kg', diff_units='s')) prob.setup(check=True, force_alloc_complex=True) prob['iv.rate_to_integrate'] = fprime prob.run_model() - assert_rel_error(self, prob.get_val('integral.q', units='kg'), f, tolerance=1e-14) - assert_rel_error(self, prob.get_val('integral.q_final', units='kg'), f[-1], tolerance=1e-14) + assert_near_equal(prob.get_val('integral.q', units='kg'), f, tolerance=1e-14) + assert_near_equal(prob.get_val('integral.q_final', units='kg'), f[-1], tolerance=1e-14) partials = prob.check_partials(method='cs',compact_print=True) assert_check_partials(partials, atol=1e-8, rtol=1e0) def test_quadratic_duration(self): - n_int_per_seg = self.num_intervals - nn_tot = (n_int_per_seg*2 + 1) + num_nodes = self.num_nodes + nn_tot = num_nodes x = np.linspace(0, nn_tot-1, nn_tot) fprime = 4 * x **2 - 8*x + 5 f = 4 * x ** 3 / 3 - 8 * x ** 2 / 2 + 5*x - prob = Problem(MultiPhaseIntegratorTestGroup(num_intervals=self.num_intervals, integrator=self.integrator, + prob = Problem(MultiPhaseIntegratorTestGroup(num_nodes=self.num_nodes, integrator=self.integrator, quantity_units='kg', diff_units='s',time_setup='duration')) prob.setup(check=True, force_alloc_complex=True) prob['iv.rate_to_integrate'] = fprime prob.run_model() - assert_rel_error(self, prob.get_val('integral.q', units='kg'), f, tolerance=1e-14) - assert_rel_error(self, prob.get_val('integral.q_final', units='kg'), f[-1], tolerance=1e-14) + assert_near_equal(prob.get_val('integral.q', units='kg'), f, tolerance=1e-14) + assert_near_equal(prob.get_val('integral.q_final', units='kg'), f[-1], tolerance=1e-14) partials = prob.check_partials(method='cs',compact_print=True) assert_check_partials(partials, atol=1e-8, rtol=1e0) @pytest.mark.filterwarnings("ignore:Input*:UserWarning") def test_quadratic_bounds(self): - n_int_per_seg = self.num_intervals - nn_tot = (n_int_per_seg*2 + 1) + num_nodes = self.num_nodes + nn_tot = num_nodes x = np.linspace(0, nn_tot-1, nn_tot) fprime = 4 * x **2 - 8*x + 5 f = 4 * x ** 3 / 3 - 8 * x ** 2 / 2 + 5*x - prob = Problem(MultiPhaseIntegratorTestGroup(num_intervals=self.num_intervals, integrator=self.integrator, + prob = Problem(MultiPhaseIntegratorTestGroup(num_nodes=self.num_nodes, integrator=self.integrator, quantity_units='kg', diff_units='s',time_setup='bounds')) prob.setup(check=True, force_alloc_complex=True) prob['iv.rate_to_integrate'] = fprime prob.run_model() - assert_rel_error(self, prob.get_val('integral.q', units='kg'), f, tolerance=1e-14) - assert_rel_error(self, prob.get_val('integral.q_final', units='kg'), f[-1], tolerance=1e-14) + assert_near_equal(prob.get_val('integral.q', units='kg'), f, tolerance=1e-14) + assert_near_equal(prob.get_val('integral.q_final', units='kg'), f[-1], tolerance=1e-14) partials = prob.check_partials(method='cs',compact_print=True) assert_check_partials(partials, atol=1e-8, rtol=1e0) def test_quadratic_single_phase_no_rate_units(self): - n_int_per_seg = self.num_intervals - nn_tot = (n_int_per_seg*2 + 1) + num_nodes = self.num_nodes + nn_tot = num_nodes x = np.linspace(0, nn_tot-1, nn_tot) fprime = 4 * x **2 - 8*x + 5 f = 4 * x ** 3 / 3 - 8 * x ** 2 / 2 + 5*x - prob = Problem(MultiPhaseIntegratorTestGroup(num_intervals=self.num_intervals, integrator=self.integrator, + prob = Problem(MultiPhaseIntegratorTestGroup(num_nodes=self.num_nodes, integrator=self.integrator, diff_units='s')) prob.setup(check=True, force_alloc_complex=True) prob['iv.rate_to_integrate'] = fprime prob.run_model() - assert_rel_error(self, prob.get_val('integral.q', units=None), f, tolerance=1e-14) - assert_rel_error(self, prob.get_val('integral.q_final', units=None), f[-1], tolerance=1e-14) + assert_near_equal(prob.get_val('integral.q', units=None), f, tolerance=1e-14) + assert_near_equal(prob.get_val('integral.q_final', units=None), f[-1], tolerance=1e-14) partials = prob.check_partials(method='cs',compact_print=True) assert_check_partials(partials, atol=1e-8, rtol=1e0) def test_quadratic_three_phase_units_equal_dt(self): - n_int_per_seg = self.num_intervals - nn_tot = (n_int_per_seg*2 + 1) + num_nodes = self.num_nodes + nn_tot = num_nodes x1 = np.linspace(0, nn_tot-1, nn_tot) x2 = np.linspace(nn_tot-1, 2*(nn_tot-1), nn_tot) x3 = np.linspace(2*(nn_tot-1), 3*(nn_tot-1), nn_tot) @@ -211,19 +209,19 @@ def test_quadratic_three_phase_units_equal_dt(self): fprime = 4 * x **2 - 8*x + 5 f = 4 * x ** 3 / 3 - 8 * x ** 2 / 2 + 5*x prob = Problem(MultiPhaseIntegratorTestGroup(segment_names=['climb','cruise','descent'], - num_intervals=self.num_intervals, integrator=self.integrator, + num_nodes=self.num_nodes, integrator=self.integrator, quantity_units='kg', diff_units='s')) prob.setup(check=True, force_alloc_complex=True) prob['iv.rate_to_integrate'] = fprime prob.run_model() - assert_rel_error(self, prob.get_val('integral.q', units='kg'), f, tolerance=1e-14) - assert_rel_error(self, prob.get_val('integral.q_final', units='kg'), f[-1], tolerance=1e-14) + assert_near_equal(prob.get_val('integral.q', units='kg'), f, tolerance=1e-14) + assert_near_equal(prob.get_val('integral.q_final', units='kg'), f[-1], tolerance=1e-14) partials = prob.check_partials(method='cs',compact_print=True) assert_check_partials(partials, atol=1e-8, rtol=1e0) def test_quadratic_three_phase_units_equal_dt_count_all(self): - n_int_per_seg = self.num_intervals - nn_tot = (n_int_per_seg*2 + 1) + num_nodes = self.num_nodes + nn_tot = num_nodes x1 = np.linspace(0, nn_tot-1, nn_tot) x2 = np.linspace(nn_tot-1, 2*(nn_tot-1), nn_tot) x3 = np.linspace(2*(nn_tot-1), 3*(nn_tot-1), nn_tot) @@ -232,19 +230,19 @@ def test_quadratic_three_phase_units_equal_dt_count_all(self): f = 4 * x ** 3 / 3 - 8 * x ** 2 / 2 + 5*x prob = Problem(MultiPhaseIntegratorTestGroup(segment_names=['climb','cruise','descent'], segments_to_count=['climb','cruise','descent'], - num_intervals=self.num_intervals, integrator=self.integrator, + num_nodes=self.num_nodes, integrator=self.integrator, quantity_units='kg', diff_units='s')) prob.setup(check=True, force_alloc_complex=True) prob['iv.rate_to_integrate'] = fprime prob.run_model() - assert_rel_error(self, prob.get_val('integral.q', units='kg'), f, tolerance=1e-14) - assert_rel_error(self, prob.get_val('integral.q_final', units='kg'), f[-1], tolerance=1e-14) + assert_near_equal(prob.get_val('integral.q', units='kg'), f, tolerance=1e-14) + assert_near_equal(prob.get_val('integral.q_final', units='kg'), f[-1], tolerance=1e-14) partials = prob.check_partials(method='cs',compact_print=True) assert_check_partials(partials, atol=1e-8, rtol=1e0) def test_quadratic_three_phase_units_equal_dt_skip_middle(self): - n_int_per_seg = self.num_intervals - nn_tot = (n_int_per_seg*2 + 1) + num_nodes = self.num_nodes + nn_tot = num_nodes x1 = np.linspace(0, nn_tot-1, nn_tot) x2 = np.linspace(nn_tot-1, 2*(nn_tot-1), nn_tot) x3 = np.linspace(2*(nn_tot-1), 3*(nn_tot-1), nn_tot) @@ -256,19 +254,19 @@ def test_quadratic_three_phase_units_equal_dt_skip_middle(self): fint[2*nn_tot:] - np.ones((nn_tot,))*(fint[2*nn_tot]-fint[nn_tot])]) prob = Problem(MultiPhaseIntegratorTestGroup(segment_names=['climb','cruise','descent'], segments_to_count=['climb','descent'], - num_intervals=self.num_intervals, integrator=self.integrator, + num_nodes=self.num_nodes, integrator=self.integrator, quantity_units='kg', diff_units='s')) prob.setup(check=True, force_alloc_complex=True) prob['iv.rate_to_integrate'] = fprime prob.run_model() - assert_rel_error(self, prob.get_val('integral.q', units='kg'), f, tolerance=1e-14) - assert_rel_error(self, prob.get_val('integral.q_final', units='kg'), f[-1], tolerance=1e-14) + assert_near_equal(prob.get_val('integral.q', units='kg'), f, tolerance=1e-14) + assert_near_equal(prob.get_val('integral.q_final', units='kg'), f[-1], tolerance=1e-14) partials = prob.check_partials(method='cs',compact_print=True) assert_check_partials(partials, atol=1e-8, rtol=1e0) def test_quadratic_three_phase_units_unequal_dt(self): - n_int_per_seg = self.num_intervals - nn_tot = (n_int_per_seg*2 + 1) + num_nodes = self.num_nodes + nn_tot = num_nodes x1 = np.linspace(0, nn_tot-1, nn_tot) x2 = np.linspace(nn_tot-1, 3*(nn_tot-1), nn_tot) x3 = np.linspace(3*(nn_tot-1), 6*(nn_tot-1), nn_tot) @@ -276,7 +274,7 @@ def test_quadratic_three_phase_units_unequal_dt(self): fprime = 4 * x **2 - 8*x + 5 f = 4 * x ** 3 / 3 - 8 * x ** 2 / 2 + 5*x prob = Problem(MultiPhaseIntegratorTestGroup(segment_names=['climb','cruise','descent'], - num_intervals=self.num_intervals, integrator=self.integrator, + num_nodes=self.num_nodes, integrator=self.integrator, quantity_units='kg', diff_units='s')) prob.setup(check=True, force_alloc_complex=True) prob['iv.rate_to_integrate'] = fprime @@ -284,14 +282,14 @@ def test_quadratic_three_phase_units_unequal_dt(self): prob['iv.cruise|dt'] = 2 prob['iv.descent|dt'] = 3 prob.run_model() - assert_rel_error(self, prob.get_val('integral.q', units='kg'), f, tolerance=1e-14) - assert_rel_error(self, prob.get_val('integral.q_final', units='kg'), f[-1], tolerance=1e-14) + assert_near_equal(prob.get_val('integral.q', units='kg'), f, tolerance=1e-14) + assert_near_equal(prob.get_val('integral.q_final', units='kg'), f[-1], tolerance=1e-14) partials = prob.check_partials(method='cs',compact_print=True) assert_check_partials(partials, atol=1e-8, rtol=1e0) def test_quadratic_three_phase_units_unequal_dt_initial_val(self): - n_int_per_seg = self.num_intervals - nn_tot = (n_int_per_seg*2 + 1) + num_nodes = self.num_nodes + nn_tot = num_nodes x1 = np.linspace(0, nn_tot-1, nn_tot) x2 = np.linspace(nn_tot-1, 3*(nn_tot-1), nn_tot) x3 = np.linspace(3*(nn_tot-1), 6*(nn_tot-1), nn_tot) @@ -300,7 +298,7 @@ def test_quadratic_three_phase_units_unequal_dt_initial_val(self): fprime = 4 * x **2 - 8*x + 5 f = 4 * x ** 3 / 3 - 8 * x ** 2 / 2 + 5*x + C prob = Problem(MultiPhaseIntegratorTestGroup(segment_names=['climb','cruise','descent'], - num_intervals=self.num_intervals, integrator=self.integrator, + num_nodes=self.num_nodes, integrator=self.integrator, quantity_units='kg', diff_units='s')) prob.setup(check=True, force_alloc_complex=True) prob['iv.rate_to_integrate'] = fprime @@ -309,8 +307,8 @@ def test_quadratic_three_phase_units_unequal_dt_initial_val(self): prob['iv.descent|dt'] = 3 prob['iv.initial_value'] = C prob.run_model() - assert_rel_error(self, prob.get_val('integral.q', units='kg'), f, tolerance=1e-14) - assert_rel_error(self, prob.get_val('integral.q_final', units='kg'), f[-1], tolerance=1e-14) + assert_near_equal(prob.get_val('integral.q', units='kg'), f, tolerance=1e-14) + assert_near_equal(prob.get_val('integral.q_final', units='kg'), f[-1], tolerance=1e-14) partials = prob.check_partials(method='cs',compact_print=True) assert_check_partials(partials, atol=1e-8, rtol=1e0) @@ -320,7 +318,7 @@ class SimpsonIntegratorEveryNode5PtTestCases(unittest.TestCase, IntegratorEveryN it cannot differentiate the quartic accurately """ def __init__(self, *args, **kwargs): - self.num_intervals = 5 + self.num_nodes = 11 self.integrator = 'simpson' super(SimpsonIntegratorEveryNode5PtTestCases, self).__init__(*args, **kwargs) @@ -330,7 +328,7 @@ class SimpsonIntegratorEveryNode3PtTestCases(unittest.TestCase, IntegratorEveryN it cannot differentiate the quartic accurately """ def __init__(self, *args, **kwargs): - self.num_intervals = 3 + self.num_nodes = 7 self.integrator = 'simpson' super(SimpsonIntegratorEveryNode3PtTestCases, self).__init__(*args, **kwargs) @@ -340,7 +338,7 @@ class BDFIntegratorEveryNode5PtTestCases(unittest.TestCase, IntegratorEveryNodeC it cannot differentiate the quartic accurately """ def __init__(self, *args, **kwargs): - self.num_intervals = 5 + self.num_nodes = 11 self.integrator = 'bdf3' super(BDFIntegratorEveryNode5PtTestCases, self).__init__(*args, **kwargs) @@ -350,6 +348,6 @@ class BDFIntegratorEveryNode3PtTestCases(unittest.TestCase, IntegratorEveryNodeC it cannot differentiate the quartic accurately """ def __init__(self, *args, **kwargs): - self.num_intervals = 3 + self.num_nodes = 7 self.integrator = 'bdf3' super(BDFIntegratorEveryNode3PtTestCases, self).__init__(*args, **kwargs) \ No newline at end of file diff --git a/openconcept/utilities/math/tests/test_math.py b/openconcept/utilities/math/tests/test_math.py index dffbb769..ddef0c44 100644 --- a/openconcept/utilities/math/tests/test_math.py +++ b/openconcept/utilities/math/tests/test_math.py @@ -2,7 +2,7 @@ import unittest import pytest import numpy as np -from openmdao.utils.assert_utils import assert_rel_error, assert_check_partials +from openmdao.utils.assert_utils import assert_near_equal, assert_check_partials from openmdao.api import IndepVarComp, Group, Problem from openconcept.utilities.math.simpson_integration import IntegrateQuantity from openconcept.utilities.math.derivatives import FirstDerivative @@ -29,7 +29,7 @@ def test_uniform_single(self): prob = Problem(SimpsonTestGroup(n_simp_intervals=1)) prob.setup(check=True) prob.run_model() - assert_rel_error(self,prob['delta_quantity'],1.,tolerance=1e-15) + assert_near_equal(prob['delta_quantity'],1.,tolerance=1e-15) partials = prob.check_partials(method='fd',compact_print=True) assert_check_partials(partials) @@ -37,7 +37,7 @@ def test_uniform(self): prob = Problem(SimpsonTestGroup(n_simp_intervals=5)) prob.setup(check=True) prob.run_model() - assert_rel_error(self,prob['delta_quantity'],1.,tolerance=1e-15) + assert_near_equal(prob['delta_quantity'],1.,tolerance=1e-15) partials = prob.check_partials(method='fd',compact_print=True) assert_check_partials(partials) @@ -46,14 +46,14 @@ def test_endpoints(self): prob.setup(check=False) prob['iv.end_pt'] = 2 prob.run_model() - assert_rel_error(self,prob['delta_quantity'],2.,tolerance=1e-15) + assert_near_equal(prob['delta_quantity'],2.,tolerance=1e-15) def test_function_level(self): prob = Problem(SimpsonTestGroup(n_simp_intervals=5)) prob.setup(check=False) prob['iv.function'] = 3*np.ones(2*5+1) prob.run_model() - assert_rel_error(self,prob['delta_quantity'],3.,tolerance=1e-15) + assert_near_equal(prob['delta_quantity'],3.,tolerance=1e-15) def test_trig_function_approx(self): prob = Problem(SimpsonTestGroup(n_simp_intervals=5)) @@ -62,7 +62,7 @@ def test_trig_function_approx(self): prob['iv.end_pt'] = np.pi prob['iv.function'] = np.sin(x) prob.run_model() - assert_rel_error(self,prob['delta_quantity'],2.,tolerance=1e-4) + assert_near_equal(prob['delta_quantity'],2.,tolerance=1e-4) partials = prob.check_partials(method='fd',compact_print=True) assert_check_partials(partials) @@ -72,7 +72,7 @@ def test_cubic_polynomial_exact(self): x = np.linspace(0,1,2*5+1) prob['iv.function'] = x**3-2*x**2+2.5*x-4 prob.run_model() - assert_rel_error(self,prob['delta_quantity'],-3.1666666666666666666666666666666666666,tolerance=1e-16) + assert_near_equal(prob['delta_quantity'],-3.1666666666666666666666666666666666666,tolerance=1e-16) partials = prob.check_partials(method='fd',compact_print=True) assert_check_partials(partials) @@ -139,7 +139,7 @@ def test_uniform_single_phase_no_units(self): prob.run_model() n_int_per_seg = 5 nn_tot = (n_int_per_seg*2 + 1) - assert_rel_error(self, prob['derivative.dqdt'], np.zeros((nn_tot,)),tolerance=1e-14) + assert_near_equal(prob['derivative.dqdt'], np.zeros((nn_tot,)),tolerance=1e-14) partials = prob.check_partials(method='cs',compact_print=True) assert_check_partials(partials, atol=1e-8, rtol=1e2) @@ -151,7 +151,7 @@ def test_linear_single_phase_no_units(self): prob['iv.quant_to_diff'] = np.linspace(0,1,nn_tot) prob['iv.dt'] = 1 / (nn_tot - 1) prob.run_model() - assert_rel_error(self, prob['derivative.dqdt'], np.ones((nn_tot,)),tolerance=1e-15) + assert_near_equal(prob['derivative.dqdt'], np.ones((nn_tot,)),tolerance=1e-15) partials = prob.check_partials(method='cs',compact_print=True) assert_check_partials(partials, atol=1e-8, rtol=1e2) @@ -166,7 +166,7 @@ def test_quadratic_single_phase_no_units(self): prob['iv.quant_to_diff'] = f_test prob['iv.dt'] = 2 / (nn_tot - 1) prob.run_model() - assert_rel_error(self, prob['derivative.dqdt'], fp_exact, tolerance=1e-14) + assert_near_equal(prob['derivative.dqdt'], fp_exact, tolerance=1e-14) partials = prob.check_partials(method='cs',compact_print=True) assert_check_partials(partials, atol=1e-8, rtol=1e-8) @@ -181,7 +181,7 @@ def test_quadratic_single_phase_units(self): prob['iv.quant_to_diff'] = f_test prob['iv.dt'] = 2 / (nn_tot - 1) prob.run_model() - assert_rel_error(self, prob.get_val('derivative.dqdt','m/s'), fp_exact, tolerance=1e-14) + assert_near_equal(prob.get_val('derivative.dqdt','m/s'), fp_exact, tolerance=1e-14) partials = prob.check_partials(method='cs',compact_print=True) assert_check_partials(partials, atol=1e-8, rtol=1e-8) @@ -196,7 +196,7 @@ def test_quadratic_single_phase_diff_units_only(self): prob['iv.quant_to_diff'] = f_test prob['iv.dt'] = 2 / (nn_tot - 1) prob.run_model() - assert_rel_error(self, prob.get_val('derivative.dqdt','s ** -1'), fp_exact, tolerance=1e-14) + assert_near_equal(prob.get_val('derivative.dqdt','s ** -1'), fp_exact, tolerance=1e-14) partials = prob.check_partials(method='cs',compact_print=True) assert_check_partials(partials, atol=1e-8, rtol=1e-8) @@ -211,7 +211,7 @@ def test_quadratic_single_named_phase_units(self): prob['iv.quant_to_diff'] = f_test prob['iv.cruise|dt'] = 2 / (nn_tot - 1) prob.run_model() - assert_rel_error(self, prob.get_val('derivative.dqdt','m/s'), fp_exact, tolerance=1e-14) + assert_near_equal(prob.get_val('derivative.dqdt','m/s'), fp_exact, tolerance=1e-14) partials = prob.check_partials(method='cs',compact_print=True) assert_check_partials(partials, atol=1e-8, rtol=1e-8) @@ -230,7 +230,7 @@ def test_quadratic_multi_phase_units(self): prob['iv.cruise|dt'] = 1 / (nn_seg - 1) prob['iv.descent|dt'] = 3 / (nn_seg - 1) prob.run_model() - assert_rel_error(self, prob.get_val('derivative.dqdt','m/s'), fp_exact, tolerance=1e-12) + assert_near_equal(prob.get_val('derivative.dqdt','m/s'), fp_exact, tolerance=1e-12) partials = prob.check_partials(method='cs',compact_print=True) assert_check_partials(partials, atol=1e-6, rtol=1e-6) @@ -249,7 +249,7 @@ def test_quadratic_multi_phase_units_7int(self): prob['iv.cruise|dt'] = 1 / (nn_seg - 1) prob['iv.descent|dt'] = 3 / (nn_seg - 1) prob.run_model() - assert_rel_error(self, prob.get_val('derivative.dqdt','m/s'), fp_exact, tolerance=1e-12) + assert_near_equal(prob.get_val('derivative.dqdt','m/s'), fp_exact, tolerance=1e-12) partials = prob.check_partials(method='cs',compact_print=True) assert_check_partials(partials, atol=1e-6, rtol=1e-6) @@ -281,7 +281,7 @@ def test_quartic_single_phase_no_units(self): prob['iv.quant_to_diff'] = f_test prob['iv.dt'] = 2 / (nn_tot - 1) prob.run_model() - assert_rel_error(self, prob['derivative.dqdt'], fp_exact, tolerance=1e-14) + assert_near_equal(prob['derivative.dqdt'], fp_exact, tolerance=1e-14) partials = prob.check_partials(method='cs',compact_print=True) assert_check_partials(partials, atol=1e-8, rtol=1e-8) @@ -296,7 +296,7 @@ def test_quartic_negative_single_phase_no_units(self): prob['iv.quant_to_diff'] = f_test prob['iv.dt'] = 2 / (nn_tot - 1) prob.run_model() - assert_rel_error(self, prob['derivative.dqdt'], fp_exact, tolerance=1e-14) + assert_near_equal(prob['derivative.dqdt'], fp_exact, tolerance=1e-14) partials = prob.check_partials(method='cs',compact_print=True) assert_check_partials(partials, atol=1e-6, rtol=1e-6) @@ -311,7 +311,7 @@ def test_quartic_single_phase_units(self): prob['iv.quant_to_diff'] = f_test prob['iv.dt'] = 2 / (nn_tot - 1) prob.run_model() - assert_rel_error(self, prob.get_val('derivative.dqdt',units='m/s'), fp_exact, tolerance=1e-14) + assert_near_equal(prob.get_val('derivative.dqdt',units='m/s'), fp_exact, tolerance=1e-14) partials = prob.check_partials(method='cs',compact_print=True) assert_check_partials(partials, atol=1e-8, rtol=1e-8) @@ -326,7 +326,7 @@ def test_quartic_single_phase_diff_units_only(self): prob['iv.quant_to_diff'] = f_test prob['iv.dt'] = 2 / (nn_tot - 1) prob.run_model() - assert_rel_error(self, prob.get_val('derivative.dqdt',units='s ** -1'), fp_exact, tolerance=1e-14) + assert_near_equal(prob.get_val('derivative.dqdt',units='s ** -1'), fp_exact, tolerance=1e-14) partials = prob.check_partials(method='cs',compact_print=True) assert_check_partials(partials, atol=1e-8, rtol=1e-8) @@ -342,7 +342,7 @@ def test_quartic_single_phase_qty_units_only(self): prob['iv.quant_to_diff'] = f_test prob['iv.dt'] = 2 / (nn_tot - 1) prob.run_model() - assert_rel_error(self, prob.get_val('derivative.dqdt',units='m'), fp_exact, tolerance=1e-14) + assert_near_equal(prob.get_val('derivative.dqdt',units='m'), fp_exact, tolerance=1e-14) partials = prob.check_partials(method='cs',compact_print=True) assert_check_partials(partials, atol=1e-8, rtol=1e-8) diff --git a/openconcept/utilities/math/tests/test_multiply_divide_comp.py b/openconcept/utilities/math/tests/test_multiply_divide_comp.py index e67c1596..159ffe93 100644 --- a/openconcept/utilities/math/tests/test_multiply_divide_comp.py +++ b/openconcept/utilities/math/tests/test_multiply_divide_comp.py @@ -7,7 +7,7 @@ from openmdao.api import Problem, Group, IndepVarComp #from openmdao.components.multiply_divide_comp import ElementMultiplyDivideComp from openconcept.utilities.math.multiply_divide_comp import ElementMultiplyDivideComp -from openmdao.utils.assert_utils import assert_rel_error, assert_check_partials +from openmdao.utils.assert_utils import assert_near_equal, assert_check_partials class TestElementMultiplyDivideCompScalars(unittest.TestCase): @@ -41,7 +41,7 @@ def test_results(self): b = self.p['b'] out = self.p['multiply_divide_comp.multdiv_output'] expected = a * b - assert_rel_error(self, out, expected,1e-16) + assert_near_equal(out, expected,1e-16) def test_partials(self): partials = self.p.check_partials(method='fd', out_stream=None) @@ -81,7 +81,7 @@ def test_results(self): b = self.p['b'] out = self.p['multiply_divide_comp.multdiv_output'] expected = a * b - assert_rel_error(self, out, expected,1e-16) + assert_near_equal(out, expected,1e-16) def test_partials(self): partials = self.p.check_partials(method='fd', out_stream=None) @@ -120,7 +120,7 @@ def test_results(self): b = self.p['b'] out = self.p['multiply_divide_comp.multdiv_output'] expected = a * b - assert_rel_error(self, out, expected,1e-16) + assert_near_equal(out, expected,1e-16) def test_partials(self): partials = self.p.check_partials(method='fd', out_stream=None) @@ -160,7 +160,7 @@ def test_results(self): b = self.p['b'] out = self.p['multiply_divide_comp.multdiv_output'] expected = a * b - assert_rel_error(self, out, expected,1e-16) + assert_near_equal(out, expected,1e-16) def test_partials(self): partials = self.p.check_partials(method='fd', out_stream=None) @@ -204,7 +204,7 @@ def test_results(self): c = self.p['c'] out = self.p['multiply_divide_comp.multdiv_output'] expected = a * b * c - assert_rel_error(self, out, expected,1e-16) + assert_near_equal(out, expected,1e-16) def test_partials(self): partials = self.p.check_partials(method='fd', out_stream=None) @@ -248,7 +248,7 @@ def test_results(self): c = self.p['c'] out = self.p['multiply_divide_comp.multdiv_output'] expected = 1 / a / b * c - assert_rel_error(self, out, expected,1e-15) + assert_near_equal(out, expected,1e-15) @pytest.mark.filterwarnings("ignore:Casting*") def test_partials(self): @@ -293,7 +293,7 @@ def test_results(self): c = self.p['c'] out = self.p['multiply_divide_comp.multdiv_output'] expected = 2 * a * b * c - assert_rel_error(self, out, expected,1e-16) + assert_near_equal(out, expected,1e-16) def test_partials(self): partials = self.p.check_partials(method='fd', out_stream=None) @@ -339,7 +339,7 @@ def test_results(self): c = self.p['c'] out = self.p.get_val('multiply_divide_comp.multdiv_output', units='kN') expected = a * b / c / 1000 - assert_rel_error(self, out, expected,1e-8) + assert_near_equal(out, expected,1e-8) def test_partials(self): partials = self.p.check_partials(method='cs', out_stream=None) @@ -385,7 +385,7 @@ def test_results(self): c = self.p['c'] out = self.p.get_val('multiply_divide_comp.multdiv_output', units='kN') expected = a * b / c / 1000 - assert_rel_error(self, out, expected,1e-8) + assert_near_equal(out, expected,1e-8) def test_partials(self): partials = self.p.check_partials(method='cs', out_stream=None) @@ -455,7 +455,7 @@ def test(self): import numpy as np #from openmdao.api import Problem, Group, IndepVarComp from openconcept.utilities.math.multiply_divide_comp import ElementMultiplyDivideComp - from openmdao.utils.assert_utils import assert_rel_error + from openmdao.utils.assert_utils import assert_near_equal n = 5 length = 4 @@ -489,7 +489,7 @@ def test(self): # Verify the results expected_i = - p['mass'] * p['acceleration'] / 1000 - assert_rel_error(self, p.get_val('inertialforcecomp.inertial_force', units='kN'), expected_i) + assert_near_equal(p.get_val('inertialforcecomp.inertial_force', units='kN'), expected_i) if __name__ == '__main__': diff --git a/openconcept/utilities/math/tests/test_new_integrals.py b/openconcept/utilities/math/tests/test_new_integrals.py new file mode 100644 index 00000000..0f21497b --- /dev/null +++ b/openconcept/utilities/math/tests/test_new_integrals.py @@ -0,0 +1,227 @@ +from __future__ import division +import unittest +import pytest +import numpy as np +from openmdao.utils.assert_utils import assert_near_equal, assert_check_partials +from openmdao.api import IndepVarComp, Group, Problem +from openconcept.utilities.math.integrals import NewIntegrator + +class MultiPhaseIntegratorTestGroup(Group): + """An OpenMDAO group to test the every-node integrator component + + Options + ------- + num_nodes : int + Number of analysis points per segment + quantity_units : str + The units of the integral quantity q (NOT the rate) + diff_units : str + The units of the integrand (none by default) + integrator : str + Which integration scheme to use (default 'simpson') + """ + + def initialize(self): + self.options.declare('quantity_units',default=None, desc="Units of the quantity being differentiated") + self.options.declare('diff_units',default=None, desc="Units of the differential") + self.options.declare('num_nodes',default=11, desc="Number of nodes per segment") + self.options.declare('integrator',default='simpson', desc="Which simpson integrator to use") + self.options.declare('time_setup',default='dt') + + def setup(self): + quantity_units = self.options['quantity_units'] + diff_units = self.options['diff_units'] + num_nodes = self.options['num_nodes'] + integrator_option = self.options['integrator'] + time_setup = self.options['time_setup'] + + num_nodes = num_nodes + + iv = self.add_subsystem('iv', IndepVarComp()) + integ = NewIntegrator(diff_units=diff_units, num_nodes=num_nodes, method=integrator_option, time_setup=time_setup) + integ.add_integrand('q', rate_name='dqdt', start_name='q_initial', end_name='q_final', + units=quantity_units) + self.add_subsystem('integral', integ) + if quantity_units is None and diff_units is None: + rate_units = None + elif quantity_units is None: + rate_units = '(' + diff_units +')** -1' + elif diff_units is None: + rate_units = quantity_units + else: + rate_units = '('+quantity_units+') / (' + diff_units +')' + + iv.add_output('rate_to_integrate', val=np.ones((num_nodes,)), units=rate_units) + iv.add_output('initial_value', val=0, units=quantity_units) + + self.connect('iv.rate_to_integrate','integral.dqdt') + self.connect('iv.initial_value', 'integral.q_initial') + + if time_setup == 'dt': + iv.add_output('dt', val=1, units=diff_units) + self.connect('iv.dt', 'integral.dt') + elif time_setup == 'duration': + iv.add_output('duration', val=1*(num_nodes-1), units=diff_units) + self.connect('iv.duration', 'integral.duration') + elif time_setup == 'bounds': + iv.add_output('t_initial', val=2, units=diff_units) + iv.add_output('t_final', val=2 + 1*(num_nodes-1), units=diff_units) + self.connect('iv.t_initial','integral.t_initial') + self.connect('iv.t_final','integral.t_final') + +class IntegratorEveryNodeCommonTestCases(object): + """ + A common set of test cases for the integrator component + """ + + def test_uniform_single_phase_no_units(self): + prob = Problem(MultiPhaseIntegratorTestGroup(num_nodes=self.num_nodes, integrator=self.integrator)) + prob.setup(check=True, force_alloc_complex=True) + prob.run_model() + num_nodes = self.num_nodes + nn_tot = num_nodes + assert_near_equal(prob['integral.q'], np.linspace(0, nn_tot-1, nn_tot), tolerance=1e-14) + assert_near_equal(prob.get_val('integral.q_final', units=None), nn_tot-1, tolerance=1e-14) + partials = prob.check_partials(method='cs',compact_print=True) + assert_check_partials(partials, atol=1e-8, rtol=1e0) + + def test_linear_single_phase_no_units(self): + num_nodes = self.num_nodes + nn_tot = num_nodes + x = np.linspace(0, nn_tot-1, nn_tot) + fprime = x + f = x ** 2 / 2 + + prob = Problem(MultiPhaseIntegratorTestGroup(num_nodes=self.num_nodes, integrator=self.integrator)) + prob.setup(check=True, force_alloc_complex=True) + prob['iv.rate_to_integrate'] = fprime + prob.run_model() + assert_near_equal(prob['integral.q'], f, tolerance=1e-14) + assert_near_equal(prob.get_val('integral.q_final', units=None), f[-1], tolerance=1e-14) + partials = prob.check_partials(method='cs',compact_print=True) + assert_check_partials(partials, atol=1e-8, rtol=1e0) + + def test_quadratic_single_phase_no_units(self): + num_nodes = self.num_nodes + nn_tot = num_nodes + x = np.linspace(0, nn_tot-1, nn_tot) + fprime = 4 * x **2 - 8*x + 5 + f = 4 * x ** 3 / 3 - 8 * x ** 2 / 2 + 5*x + + prob = Problem(MultiPhaseIntegratorTestGroup(num_nodes=self.num_nodes, integrator=self.integrator)) + prob.setup(check=True, force_alloc_complex=True) + prob['iv.rate_to_integrate'] = fprime + prob.run_model() + assert_near_equal(prob['integral.q'], f, tolerance=1e-14) + assert_near_equal(prob.get_val('integral.q_final', units=None), f[-1], tolerance=1e-14) + partials = prob.check_partials(method='cs',compact_print=True) + assert_check_partials(partials, atol=1e-8, rtol=1e0) + + def test_quadratic_single_phase_units(self): + num_nodes = self.num_nodes + nn_tot = num_nodes + x = np.linspace(0, nn_tot-1, nn_tot) + fprime = 4 * x **2 - 8*x + 5 + f = 4 * x ** 3 / 3 - 8 * x ** 2 / 2 + 5*x + + prob = Problem(MultiPhaseIntegratorTestGroup(num_nodes=self.num_nodes, integrator=self.integrator, + quantity_units='kg', diff_units='s')) + prob.setup(check=True, force_alloc_complex=True) + prob['iv.rate_to_integrate'] = fprime + prob.run_model() + assert_near_equal(prob.get_val('integral.q', units='kg'), f, tolerance=1e-14) + assert_near_equal(prob.get_val('integral.q_final', units='kg'), f[-1], tolerance=1e-14) + partials = prob.check_partials(method='cs',compact_print=True) + assert_check_partials(partials, atol=1e-8, rtol=1e0) + + def test_quadratic_duration(self): + num_nodes = self.num_nodes + nn_tot = num_nodes + x = np.linspace(0, nn_tot-1, nn_tot) + fprime = 4 * x **2 - 8*x + 5 + f = 4 * x ** 3 / 3 - 8 * x ** 2 / 2 + 5*x + + prob = Problem(MultiPhaseIntegratorTestGroup(num_nodes=self.num_nodes, integrator=self.integrator, + quantity_units='kg', diff_units='s',time_setup='duration')) + prob.setup(check=True, force_alloc_complex=True) + prob['iv.rate_to_integrate'] = fprime + prob.run_model() + assert_near_equal(prob.get_val('integral.q', units='kg'), f, tolerance=1e-14) + assert_near_equal(prob.get_val('integral.q_final', units='kg'), f[-1], tolerance=1e-14) + partials = prob.check_partials(method='cs',compact_print=True) + assert_check_partials(partials, atol=1e-8, rtol=1e0) + + @pytest.mark.filterwarnings("ignore:Input*:UserWarning") + def test_quadratic_bounds(self): + num_nodes = self.num_nodes + nn_tot = num_nodes + x = np.linspace(0, nn_tot-1, nn_tot) + fprime = 4 * x **2 - 8*x + 5 + f = 4 * x ** 3 / 3 - 8 * x ** 2 / 2 + 5*x + + prob = Problem(MultiPhaseIntegratorTestGroup(num_nodes=self.num_nodes, integrator=self.integrator, + quantity_units='kg', diff_units='s',time_setup='bounds')) + prob.setup(check=True, force_alloc_complex=True) + prob['iv.rate_to_integrate'] = fprime + prob.run_model() + assert_near_equal(prob.get_val('integral.q', units='kg'), f, tolerance=1e-14) + assert_near_equal(prob.get_val('integral.q_final', units='kg'), f[-1], tolerance=1e-14) + partials = prob.check_partials(method='cs',compact_print=True) + assert_check_partials(partials, atol=1e-8, rtol=1e0) + + def test_quadratic_single_phase_no_rate_units(self): + num_nodes = self.num_nodes + nn_tot = num_nodes + x = np.linspace(0, nn_tot-1, nn_tot) + fprime = 4 * x **2 - 8*x + 5 + f = 4 * x ** 3 / 3 - 8 * x ** 2 / 2 + 5*x + + prob = Problem(MultiPhaseIntegratorTestGroup(num_nodes=self.num_nodes, integrator=self.integrator, + diff_units='s')) + prob.setup(check=True, force_alloc_complex=True) + prob['iv.rate_to_integrate'] = fprime + prob.run_model() + assert_near_equal(prob.get_val('integral.q', units=None), f, tolerance=1e-14) + assert_near_equal(prob.get_val('integral.q_final', units=None), f[-1], tolerance=1e-14) + partials = prob.check_partials(method='cs',compact_print=True) + assert_check_partials(partials, atol=1e-8, rtol=1e0) + +class SimpsonIntegratorEveryNode5PtTestCases(unittest.TestCase, IntegratorEveryNodeCommonTestCases): + """ + Only run the common test cases using second order accuracy because + it cannot differentiate the quartic accurately + """ + def __init__(self, *args, **kwargs): + self.num_nodes = 11 + self.integrator = 'simpson' + super(SimpsonIntegratorEveryNode5PtTestCases, self).__init__(*args, **kwargs) + +class SimpsonIntegratorEveryNode3PtTestCases(unittest.TestCase, IntegratorEveryNodeCommonTestCases): + """ + Only run the common test cases using second order accuracy because + it cannot differentiate the quartic accurately + """ + def __init__(self, *args, **kwargs): + self.num_nodes = 7 + self.integrator = 'simpson' + super(SimpsonIntegratorEveryNode3PtTestCases, self).__init__(*args, **kwargs) + +class BDFIntegratorEveryNode5PtTestCases(unittest.TestCase, IntegratorEveryNodeCommonTestCases): + """ + Only run the common test cases using second order accuracy because + it cannot differentiate the quartic accurately + """ + def __init__(self, *args, **kwargs): + self.num_nodes = 11 + self.integrator = 'bdf3' + super(BDFIntegratorEveryNode5PtTestCases, self).__init__(*args, **kwargs) + +class BDFIntegratorEveryNode3PtTestCases(unittest.TestCase, IntegratorEveryNodeCommonTestCases): + """ + Only run the common test cases using second order accuracy because + it cannot differentiate the quartic accurately + """ + def __init__(self, *args, **kwargs): + self.num_nodes = 7 + self.integrator = 'bdf3' + super(BDFIntegratorEveryNode3PtTestCases, self).__init__(*args, **kwargs) \ No newline at end of file diff --git a/openconcept/utilities/math/tests/test_summation_comp.py b/openconcept/utilities/math/tests/test_summation_comp.py index 32098b7a..f6ca84ca 100644 --- a/openconcept/utilities/math/tests/test_summation_comp.py +++ b/openconcept/utilities/math/tests/test_summation_comp.py @@ -7,7 +7,7 @@ from openmdao.api import Problem, Group, IndepVarComp #from openmdao.components.sum_comp import SumComp from openconcept.utilities.math.sum_comp import SumComp -from openmdao.utils.assert_utils import assert_rel_error, assert_check_partials +from openmdao.utils.assert_utils import assert_near_equal, assert_check_partials class TestSummation1x1(unittest.TestCase): # this test case is nonsensical but should still pass @@ -37,7 +37,7 @@ def test_results(self): a = self.p['a'] out = self.p['sum_comp.sum_output'] expected = np.sum(a, axis=0) - assert_rel_error(self, out, expected,1e-16) + assert_near_equal(out, expected,1e-16) def test_partials(self): partials = self.p.check_partials(method='cs', out_stream=None) @@ -71,7 +71,7 @@ def test_results(self): a = self.p['a'] out = self.p['sum_comp.sum_output'] expected = np.sum(a, axis=None) - assert_rel_error(self, out, expected,1e-16) + assert_near_equal(out, expected,1e-16) def test_partials(self): partials = self.p.check_partials(method='cs', out_stream=None) @@ -104,7 +104,7 @@ def test_results(self): a = self.p['a'] out = self.p['sum_comp.sum_output'] expected = np.sum(a, axis=0) - assert_rel_error(self, out, expected,1e-16) + assert_near_equal(out, expected,1e-16) def test_partials(self): partials = self.p.check_partials(method='cs', out_stream=None) @@ -137,7 +137,7 @@ def test_results(self): a = self.p['a'] out = self.p['sum_comp.sum_output'] expected = np.sum(a, axis=None) - assert_rel_error(self, out, expected,1e-16) + assert_near_equal(out, expected,1e-16) def test_partials(self): partials = self.p.check_partials(method='cs', out_stream=None) @@ -173,7 +173,7 @@ def test_results(self): out = self.p['sum_comp.sum_output'] expected = self.sf*np.sum(a, axis=0) expected = expected.reshape((1,self.length)) - assert_rel_error(self, out, expected,1e-16) + assert_near_equal(out, expected,1e-16) def test_partials(self): partials = self.p.check_partials(method='cs', out_stream=None) @@ -209,7 +209,7 @@ def test_results(self): out = self.p['sum_comp.sum_output'] expected = self.sf*np.sum(a, axis=1) expected = expected.reshape((self.nn,)) - assert_rel_error(self, out, expected,1e-16) + assert_near_equal(out, expected,1e-16) def test_partials(self): partials = self.p.check_partials(method='cs', out_stream=None) @@ -245,7 +245,7 @@ def test_results(self): out = self.p['sum_comp.sum_output'] expected = self.sf*np.sum(a, axis=None) expected = expected.reshape((1,)) - assert_rel_error(self, out, expected,1e-16) + assert_near_equal(out, expected,1e-16) def test_partials(self): partials = self.p.check_partials(method='cs', out_stream=None) @@ -290,8 +290,8 @@ def test_results(self): expected_1 = expected_1.reshape((1,self.length)) expected_2 = np.sum(b, axis=0) * 1000 expected_2 = expected_2.reshape((1,self.length)) - assert_rel_error(self, out_1, expected_1,1e-15) - assert_rel_error(self, out_2, expected_2,1e-15) + assert_near_equal(out_1, expected_1,1e-15) + assert_near_equal(out_2, expected_2,1e-15) def test_partials(self): partials = self.p.check_partials(method='cs', out_stream=None) @@ -346,7 +346,7 @@ def test_results(self): out = self.p['sum_comp.sum_output'] expected = np.sum(a, axis=0) expected = expected.reshape((1,self.length)) - assert_rel_error(self, out, expected,1e-16) + assert_near_equal(out, expected,1e-16) def test_partials(self): partials = self.p.check_partials(method='cs', out_stream=None) @@ -361,7 +361,7 @@ def test(self): import numpy as np #from openmdao.api import Problem, Group, IndepVarComp from openconcept.utilities.math.sum_comp import SumComp - from openmdao.utils.assert_utils import assert_rel_error + from openmdao.utils.assert_utils import assert_near_equal n = 10 length = 1 @@ -393,7 +393,7 @@ def test(self): # Verify the results expected_i = np.sum(p['fuel_burn_by_seg'],axis=0) expected_i = expected_i.reshape((1,)) - assert_rel_error(self, p.get_val('totalfuelcomp.total_fuel', units='kg'), expected_i) + assert_near_equal(p.get_val('totalfuelcomp.total_fuel', units='kg'), expected_i) if __name__ == '__main__': diff --git a/openconcept/utilities/tests/test_dict_indepvarcomp.py b/openconcept/utilities/tests/test_dict_indepvarcomp.py index b17075f5..d3d5616d 100644 --- a/openconcept/utilities/tests/test_dict_indepvarcomp.py +++ b/openconcept/utilities/tests/test_dict_indepvarcomp.py @@ -7,7 +7,7 @@ from openmdao.api import Problem, Group #from openmdao.components.add_subtract_comp import AddSubtractComp from openconcept.utilities.dict_indepvarcomp import DictIndepVarComp -from openmdao.utils.assert_utils import assert_rel_error, assert_check_partials +from openmdao.utils.assert_utils import assert_near_equal, assert_check_partials # Test data to load # aero data ============================= @@ -56,20 +56,20 @@ def test_results(self): b = self.p['b'] out = self.p['geom|S_ref'] expected = 20 - assert_rel_error(self, out, expected,1e-16) + assert_near_equal(out, expected,1e-16) def test_units(self): out = self.p.get_val('geom|S_ref', units='m**2') expected = 20 * 0.3048**2 - assert_rel_error(self, out, expected,1e-4) + assert_near_equal(out, expected,1e-4) def test_twodeep(self): out = self.p.get_val('aero|CLmax|flaps30') expected = 1.7 - assert_rel_error(self, out, expected,1e-4) + assert_near_equal(out, expected,1e-4) out = self.p.get_val('aero|CLmax|flaps10') expected = 1.5 - assert_rel_error(self, out, expected,1e-4) + assert_near_equal(out, expected,1e-4) def test_partials(self): partials = self.p.check_partials(method='fd', out_stream=None) @@ -123,7 +123,7 @@ def test_no_value(self): # b = self.p['b'] # out = self.p['add_subtract_comp.adder_output'] # expected = a + b -# assert_rel_error(self, out, expected,1e-16) +# assert_near_equal(out, expected,1e-16) # def test_partials(self): # partials = self.p.check_partials(method='fd', out_stream=None) @@ -164,7 +164,7 @@ def test_no_value(self): # b = self.p['b'] # out = self.p['add_subtract_comp.adder_output'] # expected = a + b -# assert_rel_error(self, out, expected,1e-16) +# assert_near_equal(out, expected,1e-16) # def test_partials(self): # partials = self.p.check_partials(method='fd', out_stream=None) @@ -208,7 +208,7 @@ def test_no_value(self): # c = self.p['c'] # out = self.p['add_subtract_comp.adder_output'] # expected = a + b + c -# assert_rel_error(self, out, expected,1e-16) +# assert_near_equal(out, expected,1e-16) # def test_partials(self): # partials = self.p.check_partials(method='fd', out_stream=None) @@ -252,7 +252,7 @@ def test_no_value(self): # c = self.p['c'] # out = self.p['add_subtract_comp.adder_output'] # expected = 2*a + b - c -# assert_rel_error(self, out, expected,1e-16) +# assert_near_equal(out, expected,1e-16) # def test_partials(self): # partials = self.p.check_partials(method='fd', out_stream=None) @@ -297,7 +297,7 @@ def test_no_value(self): # out = self.p['add_subtract_comp.adder_output'] # m_to_ft = 3.280839895 # expected = a + b*m_to_ft + c*m_to_ft -# assert_rel_error(self, out, expected,1e-8) +# assert_near_equal(out, expected,1e-8) # def test_partials(self): # partials = self.p.check_partials(method='fd', out_stream=None) @@ -339,7 +339,7 @@ def test_no_value(self): # import numpy as np # #from openmdao.api import Problem, Group, IndepVarComp # from openconcept.utilities.math.add_subtract_comp import AddSubtractComp -# from openmdao.utils.assert_utils import assert_rel_error +# from openmdao.utils.assert_utils import assert_near_equal # n = 3 @@ -381,7 +381,7 @@ def test_no_value(self): # # Verify the results # expected_i = np.array([[100, 200, 300], [0, -1, -2]]).T -# assert_rel_error(self, p.get_val('totalforcecomp.total_force', units='kN'), expected_i) +# assert_near_equal(p.get_val('totalforcecomp.total_force', units='kN'), expected_i) if __name__ == '__main__': diff --git a/openconcept/utilities/tests/test_dvlabel.py b/openconcept/utilities/tests/test_dvlabel.py index 27c55ff2..7d426410 100644 --- a/openconcept/utilities/tests/test_dvlabel.py +++ b/openconcept/utilities/tests/test_dvlabel.py @@ -6,7 +6,7 @@ from openmdao.api import Problem, Group, IndepVarComp from openconcept.utilities.dvlabel import DVLabel -from openmdao.utils.assert_utils import assert_rel_error, assert_check_partials +from openmdao.utils.assert_utils import assert_near_equal, assert_check_partials class TestBasic(unittest.TestCase): @@ -38,8 +38,8 @@ def test_results(self): b_in = self.p['b_to_be_renamed'] a_out = self.p['a'] b_out = self.p['b'] - assert_rel_error(self, a_in, a_out,1e-16) - assert_rel_error(self, b_in, b_out,1e-16) + assert_near_equal(a_in, a_out,1e-16) + assert_near_equal(b_in, b_out,1e-16) def test_partials(self): partials = self.p.check_partials(method='fd', out_stream=None) @@ -75,8 +75,8 @@ def test_results(self): b_in = self.p['b_to_be_renamed'] a_out = self.p['a'] b_out = self.p['b'] - assert_rel_error(self, a_in, a_out,1e-16) - assert_rel_error(self, b_in*2.20462, b_out,1e-5) + assert_near_equal(a_in, a_out,1e-16) + assert_near_equal(b_in*2.20462, b_out,1e-5) def test_partials(self): partials = self.p.check_partials(method='fd', out_stream=None) @@ -111,8 +111,8 @@ def test_results(self): b_in = self.p['b_to_be_renamed'] a_out = self.p['a'] b_out = self.p['b'] - assert_rel_error(self, a_in, a_out,1e-16) - assert_rel_error(self, b_in*2.20462, b_out,1e-5) + assert_near_equal(a_in, a_out,1e-16) + assert_near_equal(b_in*2.20462, b_out,1e-5) def test_partials(self): partials = self.p.check_partials(method='fd', out_stream=None)