import flow360 as fl
from flow360.examples import TutorialBETDisk

TutorialBETDisk.get_files()

project = fl.Project.from_geometry(
    TutorialBETDisk.geometry, name="Tutorial BETDisk from Python"
)
geometry = project.geometry

geometry.group_faces_by_tag("faceName")
geometry.group_edges_by_tag("edgeName")

with fl.SI_unit_system:
    cylinder1 = fl.Cylinder(
        name="cylinder1",
        axis=[1, 0, 0],
        center=[-2.0, 5.0, 0],
        outer_radius=4.0,
        height=0.6,
    )

    bet_disk_refinement = fl.AxisymmetricRefinement(
        name="BET_Disk",
        spacing_axial=0.02,
        spacing_radial=0.03,
        spacing_circumferential=0.06,
        entities=cylinder1,
    )

with fl.SI_unit_system:
    cylinder2 = fl.Cylinder(
        name="cylinder2",
        axis=[1, 0, 0],
        center=[0, 5, 0],
        outer_radius=4.1,
        height=5,
    )
    farfield = fl.AutomatedFarfield()
    meshing_params = fl.MeshingParams(
        defaults=fl.MeshingDefaults(
            surface_max_edge_length=0.5,
            boundary_layer_growth_rate=1.15,
            boundary_layer_first_layer_thickness=1e-06,
        ),
        volume_zones=[farfield],
        refinements=[
            bet_disk_refinement,
            fl.UniformRefinement(
                name="cylinder_refinement", spacing=0.1, entities=[cylinder2]
            ),
        ],
    )

with fl.SI_unit_system:
    bet_cylinder = fl.Cylinder(
        name="bet_cylinder",
        axis=[-1, 0, 0],
        center=[-2.0, 5.0, 0.0],
        outer_radius=3.81,
        height=0.4572,
    )

    xrotor_file = fl.XROTORFile(file_path=TutorialBETDisk.extra["xrotor"])
    bet_disk_model = fl.BETDisk.from_xrotor(
        file=xrotor_file,
        rotation_direction_rule="rightHand",
        omega=460 * fl.u.rpm,
        chord_ref=0.3556,
        n_loading_nodes=20,
        entities=bet_cylinder,
        angle_unit=fl.u.deg,
        length_unit=fl.u.m,
    )

with fl.SI_unit_system:
    params = fl.SimulationParams(
        meshing=meshing_params,
        reference_geometry=fl.ReferenceGeometry(
            moment_center=[0.375, 0, 0],
            moment_length=[1.26666666, 1.26666666, 1.26666666],
            area=12.5,
        ),
        operating_condition=fl.AerospaceCondition.from_mach(
            mach=0.182,
            alpha=5 * fl.u.deg,
            reference_mach=0.54,
        ),
        time_stepping=fl.Steady(
            max_steps=10000, CFL=fl.RampCFL(initial=1, final=100, ramp_steps=2000)
        ),
        models=[
            fl.Wall(
                surfaces=[
                    geometry["wing"],
                    geometry["tip"],
                ],
            ),
            fl.Freestream(surfaces=farfield.farfield),
            fl.SlipWall(surfaces=farfield.symmetry_planes),
            fl.Fluid(
                navier_stokes_solver=fl.NavierStokesSolver(
                    absolute_tolerance=1e-12,
                ),
                turbulence_model_solver=fl.SpalartAllmaras(
                    absolute_tolerance=1e-10,
                    update_jacobian_frequency=1,
                    equation_evaluation_frequency=1,
                ),
            ),
            bet_disk_model,
        ],
        outputs=[
            fl.VolumeOutput(
                output_fields=[
                    "primitiveVars",
                    "betMetrics",
                    "qcriterion",
                ],
            ),
            fl.SurfaceOutput(
                surfaces=geometry["*"],
                output_fields=[
                    "primitiveVars",
                    "Cp",
                    "Cf",
                    "CfVec",
                ],
            ),
        ],
    )

case = project.run_case(
    params=params, name="Case of tutorial BETDisk from Python", use_beta_mesher=True
)

import matplotlib.pyplot as plt

case.wait()
results = case.results

total_forces = results.total_forces.as_dataframe()
total_forces.plot("pseudo_step", ["CL", "CD"])

non_linear = results.nonlinear_residuals.as_dataframe()
non_linear.plot(
    "pseudo_step",
    ["0_cont", "1_momx", "2_momy", "3_momz", "4_energ", "5_nuHat"],
    logy=True,
)
plt.show()