How to use the openmdao.api.DirectSolver function in openmdao

#
        # Set the options for our variables.
        #
        phase.set_time_options(fix_initial=True, fix_duration=True, units='s')
        phase.add_state('x', fix_initial=True, rate_source='v', units='m')
        phase.add_state('v', fix_initial=True, fix_final=True, rate_source='u', units='m/s')

        phase.add_control('u', units='m/s**2', scaler=0.01, continuity=False, rate_continuity=False,
                          rate2_continuity=False, lower=-1.0, upper=1.0)

        #
        # Maximize distance travelled.
        #
        phase.add_objective('x', loc='final', scaler=-1)

        p.model.linear_solver = om.DirectSolver()

        #
        # Setup the problem and set our initial values.
        #
        p.setup(check=True)

        p['traj.phase0.t_initial'] = 0.0
        p['traj.phase0.t_duration'] = 1.0

        p['traj.phase0.states:x'] = phase.interpolate(ys=[0, 0.25], nodes='state_input')
        p['traj.phase0.states:v'] = phase.interpolate(ys=[0, 0], nodes='state_input')
        p['traj.phase0.controls:u'] = phase.interpolate(ys=[1, -1], nodes='control_input')

        #
        # Solve the problem.
        #

#
        # Connect the two phases
        #
        p.model.connect('phase0.t_duration', 'phase1.t_duration')

        p.model.connect('phase0.timeseries2.controls:theta', 'phase1.controls:theta')
        p.model.connect('phase0.timeseries2.states:v', 'phase1.controls:v')

        # Minimize time at the end of the phase
        # phase1.add_objective('time', loc='final', scaler=1)
        # phase1.add_boundary_constraint('S', loc='final', upper=12)

        phase1.add_objective('S', loc='final', ref=1)

        p.model.linear_solver = om.DirectSolver()
        p.setup(check=True)

        p['phase0.t_initial'] = 0.0
        p['phase0.t_duration'] = 2.0

        p['phase0.states:x'] = phase0.interpolate(ys=[0, 10], nodes='state_input')
        p['phase0.states:y'] = phase0.interpolate(ys=[10, 5], nodes='state_input')
        p['phase0.states:v'] = phase0.interpolate(ys=[0, 9.9], nodes='state_input')
        p['phase0.controls:theta'] = phase0.interpolate(ys=[5, 100], nodes='control_input')
        p['phase0.input_parameters:g'] = 9.80665

        p['phase1.states:S'] = 0.0

        dm.run_problem(p)

        expected = np.sqrt((10-0)**2 + (10 - 5)**2)

# In this case, by omitting targets, we're connecting these parameters to parameters
        # with the same name in each phase.
        traj.add_input_parameter('S', units='m**2', val=0.005)

        # Link Phases (link time and all state variables)
        traj.link_phases(phases=['ascent', 'descent'], vars=['*'])

        # Issue Connections
        p.model.connect('external_params.radius', 'size_comp.radius')
        p.model.connect('external_params.dens', 'size_comp.dens')

        p.model.connect('size_comp.mass', 'traj.input_parameters:m')
        p.model.connect('size_comp.S', 'traj.input_parameters:S')

        # Finish Problem Setup
        p.model.linear_solver = om.DirectSolver()

        p.driver.add_recorder(om.SqliteRecorder('ex_two_phase_cannonball.db'))

        p.setup()

        # Set Initial Guesses
        p.set_val('external_params.radius', 0.05, units='m')
        p.set_val('external_params.dens', 7.87, units='g/cm**3')

        p.set_val('traj.design_parameters:CD', 0.5)
        p.set_val('traj.design_parameters:CL', 0.0)
        p.set_val('traj.design_parameters:T', 0.0)

        p.set_val('traj.ascent.t_initial', 0.0)
        p.set_val('traj.ascent.t_duration', 10.0)

newton = prob.model.nonlinear_solver = om.NewtonSolver()
newton.options['atol'] = 1e-6
newton.options['rtol'] = 1e-6
newton.options['iprint'] = 2
newton.options['maxiter'] = 20
newton.options['solve_subsystems'] = True
newton.options['max_sub_solves'] = 10
newton.options['err_on_non_converge'] = True
newton.options['reraise_child_analysiserror'] = False
newton.linesearch =  om.BoundsEnforceLS()
newton.linesearch.options['bound_enforcement'] = 'scalar'
newton.linesearch.options['iprint'] = -1

prob.model.linear_solver = om.DirectSolver(assemble_jac=True)

# setup the optimization
prob.driver = om.ScipyOptimizeDriver()
prob.driver.options['optimizer'] = 'SLSQP'
prob.driver.options['debug_print'] = ['desvars', 'nl_cons', 'objs']
prob.driver.opt_settings={'Major step limit': 0.05}

prob.model.add_design_var('fan:PRdes', lower=1.20, upper=1.4)
prob.model.add_design_var('lpc:PRdes', lower=2.0, upper=4.0)
prob.model.add_design_var('OPR', lower=40.0, upper=70.0, ref0=40.0, ref=70.0)
prob.model.add_design_var('RTO:T4max', lower=3000.0, upper=3600.0, ref0=3000.0, ref=3600.0)
prob.model.add_design_var('CRZ:VjetRatio', lower=1.35, upper=1.45, ref0=1.35, ref=1.45)
prob.model.add_design_var('TR', lower=0.5, upper=0.95, ref0=0.5, ref=0.95)

prob.model.add_objective('TOC.perf.TSFC')

self.connect('Fl_O:stat:area', 'perf_calcs.A_actual')
            self.connect('Fl_O:stat:P', 'perf_calcs.Ps_actual')

        if self.options['internal_solver']:
            newton = self.nonlinear_solver = om.NewtonSolver()
            newton.options['atol'] = 1e-10
            newton.options['rtol'] = 1e-10
            newton.options['maxiter'] = 20
            newton.options['iprint'] = 2
            newton.options['solve_subsystems'] = True
            newton.options['reraise_child_analysiserror'] = False
            newton.linesearch = om.BoundsEnforceLS()
            newton.linesearch.options['bound_enforcement'] = 'scalar'

            newton.linesearch.options['iprint'] = -1
            self.linear_solver = om.DirectSolver(assemble_jac=True)

phase.add_design_parameter('S', val=49.2386, units='m**2', opt=False)
    phase.add_design_parameter('Isp', val=1600.0, units='s', opt=False)
    phase.add_design_parameter('throttle', val=1.0, opt=False)

    phase.add_boundary_constraint('h', loc='final', equals=20000, scaler=1.0E-3, units='m')
    phase.add_boundary_constraint('aero.mach', loc='final', equals=1.0)
    phase.add_boundary_constraint('gam', loc='final', equals=0.0, units='rad')

    phase.add_path_constraint(name='h', lower=100.0, upper=20000, ref=20000)
    phase.add_path_constraint(name='aero.mach', lower=0.1, upper=1.8)
    phase.add_path_constraint(name='alpha', lower=-8, upper=8)

    # Minimize time at the end of the phase
    phase.add_objective('time', loc='final', ref=1.0)

    p.model.linear_solver = DirectSolver()

    p.setup(check=True, force_alloc_complex=force_alloc_complex)

    p['phase0.t_initial'] = 0.0
    p['phase0.t_duration'] = 300.0

    p['phase0.states:r'] = phase.interpolate(ys=[0.0, 111319.54], nodes='state_input')
    p['phase0.states:h'] = phase.interpolate(ys=[100.0, 20000.0], nodes='state_input')
    p['phase0.states:v'] = phase.interpolate(ys=[135.964, 283.159], nodes='state_input')
    p['phase0.states:gam'] = phase.interpolate(ys=[0.0, 0.0], nodes='state_input')
    p['phase0.states:m'] = phase.interpolate(ys=[19030.468, 16841.431], nodes='state_input')
    p['phase0.controls:alpha'] = phase.interpolate(ys=[0.0, 0.0], nodes='control_input')

    p.run_driver()

    if SHOW_PLOTS:

coupled.ln_solver = LinearGaussSeidel()
coupled.ln_solver.options['maxiter'] = 100

## Krylov Solver - No preconditioning
# coupled.ln_solver = ScipyGMRES()
# coupled.ln_solver.options['iprint'] = 1

## Krylov Solver - LNGS preconditioning
# coupled.ln_solver = ScipyGMRES()
# coupled.ln_solver.options['iprint'] = 1
# coupled.ln_solver.preconditioner = LinearGaussSeidel()
# coupled.vlmstates.ln_solver = LinearGaussSeidel()
# coupled.spatialbeamstates.ln_solver = LinearGaussSeidel()

##Direct Solver
coupled.ln_solver = DirectSolver()

# adds the MDA to root (do not remove!)
root.add('coupled',
         coupled,
         promotes=['*'])

# Add functional components here
root.add('vlmfuncs',
        vlmfuncs_comp,
        promotes=['*'])
root.add('spatialbeamfuncs',
        spatialbeamfuncs_comp,
        promotes=['*'])
root.add('fuelburn',
        fuelburn_comp,
        promotes=['*'])

# group to converge for the impulse balance
        conv = self.add_subsystem('impulse_converge', om.Group(), promotes=['*'])

        if self.options['internal_solver']:
            newton = conv.nonlinear_solver = om.NewtonSolver()
            newton.options['maxiter'] = 30
            newton.options['atol'] = 1e-2
            newton.options['solve_subsystems'] = True
            newton.options['max_sub_solves'] = 20
            newton.options['reraise_child_analysiserror'] = False
            newton.linesearch = om.BoundsEnforceLS()
            newton.linesearch.options['bound_enforcement'] = 'scalar'
            newton.linesearch.options['iprint'] = -1
            conv.linear_solver = om.DirectSolver(assemble_jac=True)

        out_tot = SetTotal(thermo_data=thermo_data, mode='h', init_reacts=self.options['Fl_I1_elements'],
                        fl_name="Fl_O:tot")
        conv.add_subsystem('out_tot', out_tot, promotes_outputs=['Fl_O:tot:*'])
        self.connect('mix_flow.n_mix', 'out_tot.init_prod_amounts')
        self.connect('mix_flow.ht_mix', 'out_tot.h')
        # note: gets Pt from the balance comp

        out_stat = SetStatic(mode="area", thermo_data=thermo_data,
                             init_reacts=self.options['Fl_I1_elements'],
                             fl_name="Fl_O:stat")
        conv.add_subsystem('out_stat', out_stat, promotes_outputs=['Fl_O:stat:*'], promotes_inputs=['area', ])
        self.connect('mix_flow.n_mix', 'out_stat.init_prod_amounts')
        self.connect('mix_flow.W_mix','out_stat.W')
        conv.connect('Fl_O:tot:S', 'out_stat.S')
        self.connect('mix_flow.ht_mix', 'out_stat.ht')

#Specify solver settings:
        newton = self.nonlinear_solver = om.NewtonSolver()
        newton.options['atol'] = 1e-8
        newton.options['rtol'] = 1e-8
        newton.options['iprint'] = 2
        newton.options['maxiter'] = 50
        newton.options['solve_subsystems'] = True
        newton.options['max_sub_solves'] = 100
        newton.options['reraise_child_analysiserror'] = False
        # ls = newton.linesearch = BoundsEnforceLS()
        ls = newton.linesearch = om.ArmijoGoldsteinLS()
        ls.options['maxiter'] = 3
        ls.options['bound_enforcement'] = 'scalar'
        # ls.options['print_bound_enforce'] = True

        self.linear_solver = om.DirectSolver(assemble_jac=True)