core

import openmdao.api as om
import numpy as np
import matplotlib.pyplot as plt

import Inviscid as inviscid

if __name__ == "__main__":

    def plot_pressure(soln_filename=None):
        plt.plot(p.get_val('aerofoil_ccw_coordinates')[0, :], p.get_val('aerofoil_ccw_coordinates')[1, :])
        plt.plot(p.get_val('InviscidSubsystem.collocation_point_coordinates')[0, :], p.get_val('InviscidSubsystem.collocation_point_coordinates')[1, :], marker=".", linestyle='None')
        plt.plot(p.get_val('InviscidSubsystem.extra_collocation_point_coordinates')[0, :], p.get_val('InviscidSubsystem.extra_collocation_point_coordinates')[1, :], marker=".", linestyle='None')
        plt.plot(p.get_val('InviscidSubsystem.collocation_point_coordinates')[0, :], -p.get_val('pressure_coefficients'), label="TASWiG")
        if soln_filename:
            javafoil = np.loadtxt(soln_filename)
            plt.plot(javafoil[:, 0], -javafoil[:, 2], label="Javafoil")
        plt.xlim(-0.1, 1.1)
        plt.ylim(-1, 1.8)
        plt.xlabel("x/c")
        plt.ylabel("-Cp")
        plt.legend()
        plt.show()

    def nice_check_partials():
        with np.printoptions(linewidth=1024, edgeitems=20, formatter={'float_kind': '{:5.5f}'.format}):
            # formatter={'float_kind': '{:5.2f}'.format})
            p.check_partials(show_only_incorrect=True, compact_print=False)

    aerofoil_path = r"DOWNLOAD THE AEROFOIL COORDINATES AND PUT THE FILEPATH HERE"
    aerofoil = np.loadtxt(aerofoil_path).T
    num_points = aerofoil.shape[1]
    # num_points = 61
    # num_upper_weights = 3
    # num_lower_weights = 3

    p = om.Problem(model=om.Group())
    p.model.add_subsystem('InviscidSubsystem',
                          inviscid.InviscidGroup(num_points=num_points),
                          promotes_inputs=['aerofoil_ccw_coordinates', 'alpha'],
                          promotes_outputs=['external_tangential_velocities',
                                            'pressure_coefficients',
                                            'stagnation_point_position'])

    p.model.set_input_defaults('alpha', 0.5, units='deg')
    p.model.set_input_defaults('aerofoil_ccw_coordinates', aerofoil)

    p.model.add_subsystem('obj_cmp', om.ExecComp('obj = (stagnation_point_position-1.0)**2'), promotes=['stagnation_point_position', 'obj'])

    p.model.add_objective('obj')
    p.model.add_design_var('alpha', lower=0, upper=10, units='deg')

    p.driver = om.ScipyOptimizeDriver()
    p.driver.options['optimizer'] = 'SLSQP'
    #p.driver.options['maxiter'] = 100
    p.driver.options['tol'] = 1e-8
    p.driver.options['debug_print'] = ['desvars', 'ln_cons', 'nl_cons', 'objs']

    p.setup()
    print("Setup")
    p.run_model()
    print("Ran model")
    p.run_driver()
    print("Finished optimisation")