From 09f82d6d63208e3719f18febc885c2ed5720fa6f Mon Sep 17 00:00:00 2001 From: James Tocknell Date: Fri, 24 Jun 2022 12:21:46 +1000 Subject: [PATCH] Revert "Drop filter_files and jacobian_plot" This reverts commit a6abcdac8359b83a277a4611ed5728e611e14ea5. --- setup.py | 2 + src/disc_solver/analyse/jacobian_plot.py | 95 ++++++++++++++++++++++++ src/disc_solver/filter_files.py | 83 +++++++++++++++++++++ 3 files changed, 180 insertions(+) create mode 100644 src/disc_solver/analyse/jacobian_plot.py create mode 100644 src/disc_solver/filter_files.py diff --git a/setup.py b/setup.py index bc30c75f..55385ac1 100644 --- a/setup.py +++ b/setup.py @@ -59,11 +59,13 @@ "ds-combine-plot = disc_solver.analyse.combine_plot:combine_main", "ds-jacobian-eigenvalue-plot = disc_solver.analyse.compute_jacobian:jacobian_eigenvalues_main", "ds-jacobian-eigenvector-plot = disc_solver.analyse.compute_jacobian:jacobian_eigenvectors_main", + "ds-jacobian-plot = disc_solver.analyse.jacobian_plot:jacobian_main", "ds-validate-plot = disc_solver.analyse.validate_plot:validate_plot_main", "ds-hydro-check-plot = disc_solver.analyse.hydro_check_plot:hydro_check_plot_main", "ds-acc-plot = disc_solver.analyse.acc_plot:acc_main", "ds-diverge-plot = disc_solver.analyse.diverge_plot:diverge_main", "ds-conserve-plot = disc_solver.analyse.conserve_plot:conserve_main", + "ds-filter-files = disc_solver.filter_files:main", "ds-j-e-plot = disc_solver.analyse.j_e_plot:j_e_plot_main", "ds-plot-taylor-space = disc_solver.solve.taylor_space:main", "ds-component-plot = disc_solver.analyse.component_plot:plot_main", diff --git a/src/disc_solver/analyse/jacobian_plot.py b/src/disc_solver/analyse/jacobian_plot.py new file mode 100644 index 00000000..b18ef999 --- /dev/null +++ b/src/disc_solver/analyse/jacobian_plot.py @@ -0,0 +1,95 @@ +# -*- coding: utf-8 -*- +""" +Jacobian-plot command for DiscSolver +""" +import numpy as np +from numpy import degrees +import matplotlib.pyplot as plt +from scipy.linalg import eigvals + +from ..utils import ODEIndex +from .utils import ( + single_solution_plotter, common_plotting_options, analyse_main_wrapper, + get_common_plot_args, analysis_func_wrapper, plot_output_wrapper, +) + + +def plot_parser(parser): + """ + Add arguments for plot command to parser + """ + common_plotting_options(parser) + return parser + + +@analyse_main_wrapper( + "Plot jacobian for DiscSolver", + plot_parser, + cmd_parser_splitters={ + "common_plot_args": get_common_plot_args, + } +) +def jacobian_main(soln, *, soln_range, common_plot_args): + """ + Entry point for ds-jacobian-plot + """ + return jacobian_plot(soln, soln_range=soln_range, **common_plot_args) + + +@analysis_func_wrapper +def jacobian_plot( + soln, *, soln_range=None, plot_filename=None, show=False, stop=90, + figargs=None, linestyle='.', title=None, close=True, filename +): + """ + Show jacobian + """ + # pylint: disable=too-many-function-args,unexpected-keyword-arg + fig = generate_jacobian_plot( + soln, soln_range, linestyle=linestyle, stop=stop, figargs=figargs, + title=title, filename=filename, + ) + + return plot_output_wrapper( + fig, file=plot_filename, show=show, close=close + ) + + +@single_solution_plotter +def generate_jacobian_plot( + soln, *, linestyle='.', figargs=None, stop=90 +): + """ + Generate plot of jacobians + """ + if figargs is None: + figargs = {} + + jacobian_data = soln.internal_data.jacobian_data + angles = jacobian_data.angles + jacobians = jacobian_data.jacobians + npnot = np.logical_not + npand = np.logical_and + indexes = degrees(angles) <= stop + + data = np.array([eigvals(j) for j in jacobians]) + + fig, axes = plt.subplots( + nrows=2, ncols=6, constrained_layout=True, sharex=True, **figargs + ) + axes.shape = 12 + for param in ODEIndex: + ax = axes[param] + pos_data = data[:, param] >= 0 + ax.plot( + degrees(angles[npand(pos_data, indexes)]), + data[npand(pos_data, indexes), param], linestyle + "b", + ) + ax.plot( + degrees(angles[npand(npnot(pos_data), indexes)]), + - data[npand(npnot(pos_data), indexes), param], linestyle + "g", + ) + ax.set_xlabel("angle from plane (°)") + ax.set_yscale("log") + ax.set_title(param.name) + return fig diff --git a/src/disc_solver/filter_files.py b/src/disc_solver/filter_files.py new file mode 100644 index 00000000..0f4ba5e5 --- /dev/null +++ b/src/disc_solver/filter_files.py @@ -0,0 +1,83 @@ +# -*- coding: utf-8 -*- +""" +Solver component of Disc Solver +""" +import argparse +from pathlib import Path +from sys import stdin + +import logbook +from logbook.compat import redirected_warnings, redirected_logging + +from h5preserve import open as h5open + +from . import __version__ as ds_version +from .file_format import registries +from .logging import logging_options, log_handler +from .utils import ODEIndex + +log = logbook.Logger(__name__) + + +def get_level(file): + """ + Get value for level + """ + with h5open(file, registries) as f: + return max(f["run"].final_solution.solution[:, ODEIndex.v_θ]) + + +def level_wrapper(output_path, level_list): + """ + level generator based on input + """ + def output_wrapper(output_file): + """ + wrap output file for filter + """ + f = output_file.open("a") + + def filter_func(path): + """ + function for filter + """ + print(str(path), file=f) + return filter_func + + return { + float(item[0]): output_wrapper(Path(output_path, item[1])) + for item in level_list + } + + +def filter_files(*, files, levels): + """ + Filter files based on a property + """ + for file in files: + file_level = get_level(file) + for level, func in levels.items(): + if file_level >= level: + yield func(file) + + +def main(): + """ + Entry point for ds-filter-files + """ + parser = argparse.ArgumentParser( + description='Config Generator for DiscSolver' + ) + parser.add_argument( + '--version', action='version', version='%(prog)s ' + ds_version + ) + parser.add_argument("--output-path", default=".") + parser.add_argument("--level", action="append", nargs=2) + logging_options(parser) + args = vars(parser.parse_args()) + with log_handler(args), redirected_warnings(), redirected_logging(): + for output in filter_files( + files=(Path(f.strip()) for f in stdin), + levels=level_wrapper(args["output_path"], args["level"]), + ): + print(output)