diff --git a/CADRE/attitude.py b/CADRE/attitude.py index 74cf176..c79fd63 100644 --- a/CADRE/attitude.py +++ b/CADRE/attitude.py @@ -51,7 +51,7 @@ def compute(self, inputs, outputs): w_B[1, i] = np.dot(Odot_BI[0, :, i], O_BI[2, :, i]) w_B[2, i] = np.dot(Odot_BI[1, :, i], O_BI[0, :, i]) - def compute_partials(self, inputs, partials): + def _compute_partials(self, inputs): """ Calculate and save derivatives. (i.e., Jacobian) """ @@ -72,6 +72,7 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): """ Matrix-vector product with the Jacobian. """ + self._compute_partials(inputs) dw_B = d_outputs['w_B'] if mode == 'fwd': @@ -235,7 +236,7 @@ def compute(self, inputs, outputs): O_RI[1, :, i] = jB O_RI[2, :, i] = -v - def compute_partials(self, inputs, partials): + def _compute_partials(self, inputs): """ Calculate and save derivatives. (i.e., Jacobian) """ @@ -302,6 +303,7 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): """ Matrix-vector product with the Jacobian. """ + self._compute_partials(inputs) dO_RI = d_outputs['O_RI'] if mode == 'fwd': @@ -357,7 +359,7 @@ def compute(self, inputs, outputs): O_BR[1, 1, :] = O_BR[0, 0, :] O_BR[2, 2, :] = np.ones(self.n) - def compute_partials(self, inputs, partials): + def _compute_partials(self, inputs): """ Calculate and save derivatives. (i.e., Jacobian) """ @@ -373,6 +375,7 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): """ Matrix-vector product with the Jacobian. """ + self._compute_partials(inputs) dO_BR = d_outputs['O_BR'] if mode == 'fwd': @@ -563,7 +566,7 @@ def compute(self, inputs, outputs): v_e2b_B[:] = computepositionrotd(self.n, r_e2b_I[3:, :], O_BI) - def compute_partials(self, inputs, partials): + def _compute_partials(self, inputs): """ Calculate and save derivatives. (i.e., Jacobian) """ @@ -578,6 +581,7 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): """ Matrix-vector product with the Jacobian. """ + self._compute_partials(inputs) dv_e2b_B = d_outputs['v_e2b_B'] if mode == 'fwd': @@ -665,7 +669,7 @@ def compute(self, inputs, outputs): T_tot[:, i] = np.dot(self.J, wdot_B[:, i]) + \ np.dot(wx, np.dot(self.J, w_B[:, i])) - def compute_partials(self, inputs, partials): + def _compute_partials(self, inputs): """ Calculate and save derivatives. (i.e., Jacobian) """ @@ -690,6 +694,7 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): """ Matrix-vector product with the Jacobian. """ + self._compute_partials(inputs) dT_tot = d_outputs['T_tot'] if mode == 'fwd': diff --git a/CADRE/battery.py b/CADRE/battery.py index 44e480f..d3fe375 100644 --- a/CADRE/battery.py +++ b/CADRE/battery.py @@ -161,7 +161,7 @@ def compute(self, inputs, outputs): outputs['I_bat'] = P_bat / self.V - def compute_partials(self, inputs, partials): + def _compute_partials(self, inputs): """ Calculate and save derivatives. (i.e., Jacobian) """ @@ -183,7 +183,7 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): """ Matrix-vector product with the Jacobian. """ - + self._compute_partials(inputs) dI_bat = d_outputs['I_bat'] if mode == 'fwd': @@ -257,7 +257,7 @@ def compute(self, inputs, outputs): outputs['ConS0'] = KSfunction.compute(self.SOC0 - SOC, self.rho).flatten() outputs['ConS1'] = KSfunction.compute(SOC - self.SOC1, self.rho).flatten() - def compute_partials(self, inputs, partials): + def _compute_partials(self, inputs): """ Calculate and save derivatives. (i.e., Jacobian) """ @@ -280,8 +280,9 @@ def compute_partials(self, inputs, partials): def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): """ - Matrix-vector product with the Jacobian. - """ + Matrix-vector product with the Jacobian. + """ + self._compute_partials(inputs) if mode == 'fwd': if 'I_bat' in d_inputs: if 'ConCh' in d_outputs: diff --git a/CADRE/comm.py b/CADRE/comm.py index 2ddd70d..fe37881 100644 --- a/CADRE/comm.py +++ b/CADRE/comm.py @@ -97,7 +97,7 @@ def compute(self, inputs, outputs): q_A[2, :] = - np.sin(antAngle/2.) / rt2 q_A[3, :] = 0.0 - def compute_partials(self, inputs, partials): + def _compute_partials(self, inputs): """ Calculate and save derivatives. (i.e., Jacobian) """ @@ -114,6 +114,7 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): """ Matrix-vector product with the Jacobian. """ + self._compute_partials(inputs) if mode == 'fwd': if 'antAngle' in d_inputs: for k in range(4): @@ -169,7 +170,7 @@ def compute(self, inputs, outputs): O_AB[:, :, i] = np.dot(A.T, B) - def compute_partials(self, inputs, partials): + def _compute_partials(self, inputs): """ Calculate and save derivatives. (i.e., Jacobian) """ @@ -241,6 +242,7 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): """ Matrix-vector product with the Jacobian. """ + self._compute_partials(inputs) if 'q_A' in d_inputs: dq_A = d_inputs['q_A'] dO_AB = d_outputs['O_AB'] @@ -315,7 +317,7 @@ def compute(self, inputs, outputs): Dr[i] = self.alpha * P_comm[i] * gain[i] * \ CommLOS[i] / S2 ** 2 - def compute_partials(self, inputs, partials): + def _compute_partials(self, inputs): """ Calculate and save derivatives. (i.e., Jacobian) """ @@ -350,6 +352,7 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): """ Matrix-vector product with the Jacobian. """ + self._compute_partials(inputs) dDr = d_outputs['Dr'] if mode == 'fwd': @@ -402,7 +405,7 @@ def compute(self, inputs, outputs): for i in range(0, self.n): GSdist[i] = np.dot(r_b2g_A[:, i], r_b2g_A[:, i])**0.5 - def compute_partials(self, inputs, partials): + def _compute_partials(self, inputs): """ Calculate and save derivatives. (i.e., Jacobian) """ @@ -420,6 +423,7 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): """ Matrix-vector product with the Jacobian. """ + self._compute_partials(inputs) if mode == 'fwd': if 'r_b2g_A' in d_inputs: for k in range(3): @@ -462,7 +466,7 @@ def compute(self, inputs, outputs): q_E[0, :] = np.cos(theta) q_E[3, :] = -np.sin(theta) - def compute_partials(self, inputs, partials): + def _compute_partials(self, inputs): """ Calculate and save derivatives. (i.e., Jacobian) """ @@ -481,6 +485,7 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): """ Matrix-vector product with the Jacobian. """ + self._compute_partials(inputs) if mode == 'fwd': if 't' in d_inputs: for k in range(4): @@ -534,7 +539,7 @@ def compute(self, inputs, outputs): O_IE[:, :, i] = np.dot(A.T, B) - def compute_partials(self, inputs, partials): + def _compute_partials(self, inputs): """ Calculate and save derivatives. (i.e., Jacobian) """ @@ -609,6 +614,7 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): """ Matrix-vector product with the Jacobian. """ + self._compute_partials(inputs) if 'q_E' in d_inputs: dO_IE = d_outputs['O_IE'] dq_E = d_inputs['q_E'] @@ -676,7 +682,7 @@ def compute(self, inputs, outputs): self.x[:, 1] = result[1] outputs['gain'] = self.MBI.evaluate(self.x)[:, 0] - def compute_partials(self, inputs, partials): + def _compute_partials(self, inputs): """ Calculate and save derivatives. (i.e., Jacobian) """ @@ -687,6 +693,7 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): """ Matrix-vector product with the Jacobian. """ + self._compute_partials(inputs) dgain = d_outputs['gain'] if mode == 'fwd': @@ -745,7 +752,7 @@ def compute(self, inputs, outputs): r_e2g_E[1, :] = r_GS * cos_lat * np.sin(self.d2r*lon) r_e2g_E[2, :] = r_GS * np.sin(self.d2r*lat) - def compute_partials(self, inputs, partials): + def _compute_partials(self, inputs): """ Calculate and save derivatives. (i.e., Jacobian) """ @@ -780,6 +787,7 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): """ Matrix-vector product with the Jacobian. """ + self._compute_partials(inputs) dr_e2g_E = d_outputs['r_e2g_E'] if mode == 'fwd': @@ -838,7 +846,7 @@ def compute(self, inputs, outputs): for i in range(0, self.n): r_e2g_I[:, i] = np.dot(O_IE[:, :, i], r_e2g_E[:, i]) - def compute_partials(self, inputs, partials): + def _compute_partials(self, inputs): """ Calculate and save derivatives. (i.e., Jacobian) """ @@ -857,6 +865,7 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): """ Matrix-vector product with the Jacobian. """ + self._compute_partials(inputs) dr_e2g_I = d_outputs['r_e2g_I'] if mode == 'fwd': @@ -932,7 +941,7 @@ def compute(self, inputs, outputs): x = (proj - 0) / (-Rb - 0) CommLOS[i] = 3 * x ** 2 - 2 * x ** 3 - def compute_partials(self, inputs, partials): + def _compute_partials(self, inputs): """ Calculate and save derivatives. (i.e., Jacobian) """ @@ -967,6 +976,7 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): """ Matrix-vector product with the Jacobian. """ + self._compute_partials(inputs) dCommLOS = d_outputs['CommLOS'] if mode == 'fwd': @@ -1016,7 +1026,7 @@ def compute(self, inputs, outputs): outputs['r_b2g_A'] = computepositionrotd(self.n, inputs['r_b2g_B'], inputs['O_AB']) - def compute_partials(self, inputs, partials): + def _compute_partials(self, inputs): """ Calculate and save derivatives. (i.e., Jacobian) """ @@ -1027,6 +1037,7 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): """ Matrix-vector product with the Jacobian. """ + self._compute_partials(inputs) dr_b2g_A = d_outputs['r_b2g_A'] if mode == 'fwd': @@ -1090,7 +1101,7 @@ def compute(self, inputs, outputs): for i in range(0, self.n): r_b2g_B[:, i] = np.dot(O_BI[:, :, i], r_b2g_I[:, i]) - def compute_partials(self, inputs, partials): + def _compute_partials(self, inputs): """ Calculate and save derivatives. (i.e., Jacobian) """ @@ -1109,6 +1120,7 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): """ Matrix-vector product with the Jacobian. """ + self._compute_partials(inputs) dr_b2g_B = d_outputs['r_b2g_B'] if mode == 'fwd': @@ -1222,7 +1234,7 @@ def compute(self, inputs, outputs): outputs['azimuthGS'] = azimuthGS outputs['elevationGS'] = elevationGS - def compute_partials(self, inputs, partials): + def _compute_partials(self, inputs): """ Calculate and save derivatives. (i.e., Jacobian) """ @@ -1240,6 +1252,7 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): """ Matrix-vector product with the Jacobian. """ + self._compute_partials(inputs) if mode == 'fwd': if 'r_b2g_A' in d_inputs: r_b2g_A = d_inputs['r_b2g_A'].reshape((3 * self.n), order='F') diff --git a/CADRE/power.py b/CADRE/power.py index dba79ac..493b79e 100644 --- a/CADRE/power.py +++ b/CADRE/power.py @@ -90,7 +90,7 @@ def compute(self, inputs, outputs): for c in range(7): outputs['V_sol'] += self.raw[:, c, :].T - def compute_partials(self, inputs, partials): + def _compute_partials(self, inputs): """ Calculate and save derivatives. (i.e., Jacobian) """ @@ -120,6 +120,7 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): """ Matrix-vector product with the Jacobian. """ + self._compute_partials(inputs) dV_sol = d_outputs['V_sol'] if mode == 'fwd': diff --git a/CADRE/reactionwheel.py b/CADRE/reactionwheel.py index 6edbdd3..8c9f621 100644 --- a/CADRE/reactionwheel.py +++ b/CADRE/reactionwheel.py @@ -57,7 +57,7 @@ def compute(self, inputs, outputs): T_m[:, i] = -T_RW[:, i] - np.dot(w_Bx, h_RW[:, i]) - def compute_partials(self, inputs, partials): + def _compute_partials(self, inputs): """ Calculate and save derivatives. (i.e., Jacobian) """ @@ -98,6 +98,7 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): """ Matrix-vector product with the Jacobian. """ + self._compute_partials(inputs) dT_m = d_outputs['T_m'] if mode == 'fwd': @@ -165,7 +166,7 @@ def compute(self, inputs, outputs): self.b * T_RW[k, i])**2 + self.V * self.I0) - def compute_partials(self, inputs, partials): + def _compute_partials(self, inputs): """ Calculate and save derivatives. (i.e., Jacobian) """ @@ -184,6 +185,7 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): """ Matrix-vector product with the Jacobian. """ + self._compute_partials(inputs) dP_RW = d_outputs['P_RW'] if mode == 'fwd': diff --git a/CADRE/rk4.py b/CADRE/rk4.py index e9de1f5..e98454e 100644 --- a/CADRE/rk4.py +++ b/CADRE/rk4.py @@ -194,7 +194,7 @@ def compute(self, inputs, outputs): state_var_name = self.name_map['y'] outputs[state_var_name][:] = self.y.T.reshape((n_time, n_state)).T - def compute_partials(self, inputs, partials): + def _compute_partials(self): """ Calculate and save derivatives. (i.e., Jacobian) """ @@ -255,6 +255,7 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): """ Matrix-vector product with the Jacobian. """ + self._compute_partials() if mode == 'fwd': result_ext = self._applyJext(d_inputs, d_outputs) diff --git a/CADRE/solar.py b/CADRE/solar.py index 6a2ccce..ee212fe 100644 --- a/CADRE/solar.py +++ b/CADRE/solar.py @@ -121,7 +121,7 @@ def setx(self, inputs): self.x[:, 1] = result[0] self.x[:, 2] = result[1] - def compute_partials(self, inputs, partials): + def _compute_partials(self): """ Calculate and save derivatives. (i.e., Jacobian) """ @@ -133,6 +133,7 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): """ Matrix-vector product with the Jacobian. """ + self._compute_partials() deA = d_outputs['exposedArea'] if mode == 'fwd': diff --git a/CADRE/sun.py b/CADRE/sun.py index e83b91b..4a8af86 100644 --- a/CADRE/sun.py +++ b/CADRE/sun.py @@ -66,7 +66,7 @@ def compute(self, inputs, outputs): x = (dist - self.r1) / (self.r2 - self.r1) LOS[i] = 3*x**2 - 2*x**3 - def compute_partials(self, inputs, partials): + def _compute_partials(self, inputs): """ Calculate and save derivatives. (i.e., Jacobian) """ @@ -143,6 +143,7 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): """ Matrix-vector product with the Jacobian. """ + self._compute_partials(inputs) dLOS = d_outputs['LOS'] if mode == 'fwd': @@ -203,7 +204,7 @@ def compute(self, inputs, outputs): outputs['r_e2s_B'] = computepositionrotd(self.n, inputs['r_e2s_I'], inputs['O_BI']) - def compute_partials(self, inputs, partials): + def _compute_partials(self, inputs): """ Calculate and save derivatives. (i.e., Jacobian) """ @@ -214,6 +215,7 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): """ Matrix-vector product with the Jacobian. """ + self._compute_partials(inputs) dr_e2s_B = d_outputs['r_e2s_B'] if mode == 'fwd': @@ -285,7 +287,7 @@ def compute(self, inputs, outputs): r_e2s_I[1, i] = np.sin(Lambda)*np.cos(eps) r_e2s_I[2, i] = np.sin(Lambda)*np.sin(eps) - def compute_partials(self, inputs, partials): + def _compute_partials(self, inputs): """ Calculate and save derivatives. (i.e., Jacobian) """ @@ -321,6 +323,7 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): """ Matrix-vector product with the Jacobian. """ + self._compute_partials(inputs) dr_e2s_I = d_outputs['r_e2s_I'] if mode == 'fwd': @@ -372,7 +375,7 @@ def compute(self, inputs, outputs): outputs['azimuth'] = azimuth outputs['elevation'] = elevation - def compute_partials(self, inputs, partials): + def _compute_partials(self, inputs): """ Calculate and save derivatives. (i.e., Jacobian) """ @@ -389,6 +392,7 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): """ Matrix-vector product with the Jacobian. """ + self._compute_partials(inputs) if mode == 'fwd': if 'r_e2s_B' in d_inputs: r_e2s_B = d_inputs['r_e2s_B'].reshape((3*self.n), order='F') diff --git a/CADRE/test/test_mdp.py b/CADRE/test/test_mdp.py index e5bd96f..aaf8f9d 100644 --- a/CADRE/test/test_mdp.py +++ b/CADRE/test/test_mdp.py @@ -16,7 +16,7 @@ from openmdao import __version__ as om_version from openmdao.api import Problem, PETScKrylov -from openmdao.utils.mpi import MPI +from openmdao.utils.mpi import MPI, multi_proc_exception_check from CADRE.CADRE_mdp import CADRE_MDP_Group @@ -87,6 +87,7 @@ def test_CADRE_MDP(self): abs_names = model._var_allprocs_prom2abs_list['output'] checked = 0 + bad = {} for var in data: # we changed constraint names @@ -103,7 +104,7 @@ def test_CADRE_MDP(self): # make sure var is local before we try to look it up compname = abs_names[xvar][0].rsplit('.', 1)[0] comp = model._get_subsystem(compname) - if comp and not (comp.comm is None or comp.comm == MPI.COMM_NULL): + if comp and not (comp.comm is None or (MPI is not None and comp.comm == MPI.COMM_NULL)): computed = prob[xvar] actual = data[var] if isinstance(computed, np.ndarray): @@ -117,13 +118,19 @@ def test_CADRE_MDP(self): print(actual) if np.mean(actual) > 1e-3 or np.mean(computed) > 1e-3: - assert rel <= 1e-3 + if not rel <= 1e-3: + bad[xvar] = rel checked += 1 + with multi_proc_exception_check(prob.comm): + if bad: + self.fail(f"Bad norms: {bad}") + # make sure we checked everything if MPI: # objective + con5 for both points + con1-4 on this proc - self.assertEqual(checked, 7) + with multi_proc_exception_check(prob.comm): + self.assertEqual(checked, 7) else: # objective + 5 constraints for each point self.assertEqual(checked, 11) @@ -184,7 +191,12 @@ def test_CADRE_MDP(self): print(actual) if np.mean(actual) > 1e-3 or np.mean(computed) > 1e-3: - assert rel <= 1e-3 + if not rel <= 1e-3: + bad[bkey1, bkey2] = rel + + with multi_proc_exception_check(prob.comm): + if bad: + self.fail(f'Bad derivatives: {bad}') if __name__ == '__main__':