From b744085ca7d196b67f6d59f139c98f97f5d993d6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?K=C3=BChle=2C=20Laura=20Christine=20=28lakue103=29?= <laura.kuehle@uni-duesseldorf.de> Date: Tue, 31 Aug 2021 16:55:53 +0200 Subject: [PATCH] Changed maximal line length in code to 100. --- ANN_Data_Generator.py | 65 ++++++++++++++++------------ ANN_Model.py | 7 +-- ANN_Training.py | 34 +++++++++------ Basis_Function.py | 20 ++++++--- DG_Approximation.py | 35 +++++++++------ Initial_Condition.py | 3 +- Main.py | 11 ++--- Plotting.py | 6 ++- Troubled_Cell_Detector.py | 89 +++++++++++++++++++++++---------------- Update_Scheme.py | 6 ++- 10 files changed, 167 insertions(+), 109 deletions(-) diff --git a/ANN_Data_Generator.py b/ANN_Data_Generator.py index 860ef55..629b08c 100644 --- a/ANN_Data_Generator.py +++ b/ANN_Data_Generator.py @@ -14,8 +14,8 @@ import DG_Approximation class TrainingDataGenerator(object): - def __init__(self, initial_conditions, left_bound=-1, right_bound=1, balance=0.5, stencil_length=3, - distribution=None, directory=None): + def __init__(self, initial_conditions, left_bound=-1, right_bound=1, balance=0.5, + stencil_length=3, distribution=None, directory=None): self._balance = balance self._left_bound = left_bound self._right_bound = right_bound @@ -52,8 +52,8 @@ class TrainingDataGenerator(object): data = {} for set_name in self._distribution: print('Calculating ' + set_name + ' data...') - input_data, output_data = self._calculate_data_set(round(self._distribution[set_name]*num_samples), - normalize) + input_data, output_data = self._calculate_data_set( + round(self._distribution[set_name]*num_samples), normalize) data[set_name] = [input_data, output_data] print('Finished calculating ' + set_name + ' data!') @@ -74,11 +74,12 @@ class TrainingDataGenerator(object): def _calculate_data_set(self, num_samples, normalize): num_smooth_samples = round(num_samples * self._balance) - smooth_input, smooth_output = self._generate_cell_data(num_smooth_samples, self._smooth_functions, True) + smooth_input, smooth_output = self._generate_cell_data(num_smooth_samples, + self._smooth_functions, True) num_troubled_samples = num_samples - num_smooth_samples - troubled_input, troubled_output = self._generate_cell_data(num_troubled_samples, self._troubled_functions, - False) + troubled_input, troubled_output = self._generate_cell_data(num_troubled_samples, + self._troubled_functions, False) # Normalize data if normalize: @@ -119,19 +120,21 @@ class TrainingDataGenerator(object): initial_condition.induce_adjustment(-h[0]/3) left_bound, right_bound = interval - dg_scheme = DG_Approximation.DGScheme('NoDetection', polynomial_degree=polynomial_degree, - num_grid_cells=self._stencil_length, left_bound=left_bound, - right_bound=right_bound, quadrature='Gauss', - quadrature_config={'num_eval_points': polynomial_degree+1}) + dg_scheme = DG_Approximation.DGScheme( + 'NoDetection', polynomial_degree=polynomial_degree, + num_grid_cells=self._stencil_length, left_bound=left_bound, right_bound=right_bound, + quadrature='Gauss', quadrature_config={'num_eval_points': polynomial_degree+1}) if initial_condition.is_smooth(): - input_data[i] = dg_scheme.build_training_data(0, self._stencil_length, initial_condition) + input_data[i] = dg_scheme.build_training_data( + 0, self._stencil_length, initial_condition) else: - input_data[i] = dg_scheme.build_training_data(centers[self._stencil_length//2], self._stencil_length, - initial_condition) + input_data[i] = dg_scheme.build_training_data( + centers[self._stencil_length//2], self._stencil_length, initial_condition) # Update Function ID - if (i % num_function_samples == num_function_samples - 1) and (function_id != len(initial_conditions)-1): + if (i % num_function_samples == num_function_samples - 1) \ + and (function_id != len(initial_conditions)-1): function_id = function_id + 1 count += 1 @@ -166,7 +169,8 @@ class TrainingDataGenerator(object): interval = np.array([point - self._stencil_length/2 * grid_spacing, point + self._stencil_length/2 * grid_spacing]) stencil = np.array([point + factor * grid_spacing - for factor in range(-(self._stencil_length//2), self._stencil_length//2 + 1)]) + for factor in range(-(self._stencil_length//2), + self._stencil_length//2 + 1)]) return interval, stencil, grid_spacing @staticmethod @@ -182,17 +186,24 @@ class TrainingDataGenerator(object): np.random.seed(1234) boundary = [-1, 1] -functions = [{'function': Initial_Condition.Sine(boundary[0], boundary[1], {}), 'config': {'factor': 2}}, - {'function': Initial_Condition.Linear(boundary[0], boundary[1], {}), 'config': {}}, - {'function': Initial_Condition.Polynomial(boundary[0], boundary[1], {}), 'config': {}}, - {'function': Initial_Condition.Continuous(boundary[0], boundary[1], {}), 'config': {}}, - {'function': Initial_Condition.LinearAbsolut(boundary[0], boundary[1], {}), 'config': {}}, - {'function': Initial_Condition.HeavisideOneSided(boundary[0], boundary[1], {}), 'config': {}}, - {'function': Initial_Condition.HeavisideTwoSided(boundary[0], boundary[1], {}), 'config': { - 'adjustment': 0}}] - -generator = TrainingDataGenerator(functions, distribution={'train': 0.727, 'valid': 0.243, 'test': 0.03}, - left_bound=boundary[0], right_bound=boundary[1]) +functions = [{'function': Initial_Condition.Sine(boundary[0], boundary[1], {}), + 'config': {'factor': 2}}, + {'function': Initial_Condition.Linear(boundary[0], boundary[1], {}), + 'config': {}}, + {'function': Initial_Condition.Polynomial(boundary[0], boundary[1], {}), + 'config': {}}, + {'function': Initial_Condition.Continuous(boundary[0], boundary[1], {}), + 'config': {}}, + {'function': Initial_Condition.LinearAbsolut(boundary[0], boundary[1], {}), + 'config': {}}, + {'function': Initial_Condition.HeavisideOneSided(boundary[0], boundary[1], {}), + 'config': {}}, + {'function': Initial_Condition.HeavisideTwoSided(boundary[0], boundary[1], {}), + 'config': {'adjustment': 0}}] + +generator = TrainingDataGenerator( + functions, distribution={'train': 0.727, 'valid': 0.243, 'test': 0.03}, left_bound=boundary[0], + right_bound=boundary[1]) # generator = TrainingDataGenerator(functions, left_bound=boundary[0], right_bound=boundary[1]) sample_number = 66000 diff --git a/ANN_Model.py b/ANN_Model.py index 15685ba..7347911 100644 --- a/ANN_Model.py +++ b/ANN_Model.py @@ -25,13 +25,14 @@ class ThreeLayerReLu(torch.nn.Module): if not hasattr(torch.nn.modules.activation, activation_function): raise ValueError('Invalid activation function: "%s"' % activation_function) - self._name = self.__class__.__name__ + '_' + str(first_hidden_size) + '_' + str(second_hidden_size) + '_'\ - + activation_function + self._name = self.__class__.__name__ + '_' + str(first_hidden_size) + '_' \ + + str(second_hidden_size) + '_' + activation_function self.input_linear = torch.nn.Linear(input_size, first_hidden_size) self.middle_linear = torch.nn.Linear(first_hidden_size, second_hidden_size) self.output_linear = torch.nn.Linear(second_hidden_size, output_size) - self._output_layer = getattr(torch.nn.modules.activation, activation_function)(**activation_config) + self._output_layer = getattr(torch.nn.modules.activation, activation_function)( + **activation_config) def forward(self, input_data): prediction = self.input_linear(input_data).clamp(min=0) diff --git a/ANN_Training.py b/ANN_Training.py index 3171fb3..e9c32c0 100644 --- a/ANN_Training.py +++ b/ANN_Training.py @@ -32,10 +32,13 @@ class ModelTrainer(object): data_dir = config.pop('data_dir', 'data') self._model_dir = config.pop('model_dir', 'model') self._plot_dir = config.pop('plot_dir', 'fig') - self._training_file = config.pop('training_set', 'smooth_0.01k__troubled_0.01k__normalized.npy') - self._validation_file = config.pop('validation_set', 'smooth_4.0095k__troubled_4.0095k__normalized.npy') + self._training_file = config.pop('training_set', + 'smooth_0.01k__troubled_0.01k__normalized.npy') + self._validation_file = config.pop('validation_set', + 'smooth_4.0095k__troubled_4.0095k__normalized.npy') self._test_file = config.pop('test_set', 'smooth_0.495k__troubled_0.495k__normalized.npy') - file_names = {'train': self._training_file, 'validation': self._validation_file, 'test': self._test_file} + file_names = {'train': self._training_file, 'validation': self._validation_file, + 'test': self._test_file} self._read_training_data(data_dir, file_names) self._batch_size = config.pop('batch_size', 500) @@ -60,8 +63,10 @@ class ModelTrainer(object): raise ValueError('Invalid optimizer: "%s"' % self._optimizer) self._model = getattr(ANN_Model, self._model)(self._model_config) - self._loss_function = getattr(torch.nn.modules.loss, self._loss_function)(**self._loss_config) - self._optimizer = getattr(torch.optim, self._optimizer)(self._model.parameters(), **self._optimizer_config) + self._loss_function = getattr(torch.nn.modules.loss, self._loss_function)( + **self._loss_config) + self._optimizer = getattr(torch.optim, self._optimizer)( + self._model.parameters(), **self._optimizer_config) self._validation_loss = torch.zeros(self._num_epochs//100) def _read_training_data(self, directory, file_names): @@ -76,7 +81,8 @@ class ModelTrainer(object): return map(torch.tensor, (np.load(input_file), np.load(output_file))) def epoch_training(self): - # Get Training/validation Datasets from Saved Files and Map to Torch Tensor and batched-datasets + # Get Training/validation Datasets from Saved Files + # and Map to Torch Tensor and batched-datasets x_train, y_train = self._training_data['train'] train_ds = TensorDataset(x_train, y_train) train_dl = DataLoader(train_ds, batch_size=self._batch_size, shuffle=True) @@ -133,8 +139,9 @@ class ModelTrainer(object): # Set file name test_name = self._validation_file.split('.npy')[0] - name = self._model.get_name() + '__' + self._optimizer.__class__.__name__ + '_' + str(self._learning_rate)\ - + '__' + self._loss_function.__class__.__name__ + '__' + test_name + name = self._model.get_name() + '__' + self._optimizer.__class__.__name__ + '_' \ + + str(self._learning_rate) + '__' + self._loss_function.__class__.__name__ + '__' \ + + test_name # Set paths for plot files if not existing already if not os.path.exists(self._plot_dir): @@ -181,7 +188,8 @@ class ModelTrainer(object): precision = true_positive / (true_positive+false_positive) recall = true_positive / (true_positive+false_negative) - accuracy = (true_positive+true_negative) / (true_positive+true_negative+false_positive+false_negative) + accuracy = (true_positive+true_negative) / (true_positive+true_negative + + false_positive+false_negative) # print(true_positive+true_negative+false_positive+false_negative) return precision, recall, accuracy @@ -189,8 +197,9 @@ class ModelTrainer(object): # Saving Model train_name = self._training_file.split('.npy')[0] valid_name = self._validation_file.split('.npy')[0] - path = self._model.get_name() + '__' + self._optimizer.__class__.__name__ + '_' + str(self._learning_rate)\ - + '__' + self._loss_function.__class__.__name__ + '__' + train_name + '__' + valid_name + '.pt' + path = self._model.get_name() + '__' + self._optimizer.__class__.__name__ + '_' \ + + str(self._learning_rate) + '__' + self._loss_function.__class__.__name__ + '__' \ + + train_name + '__' + valid_name + '.pt' # Set paths for plot files if not existing already if not os.path.exists(self._model_dir): @@ -203,7 +212,8 @@ class ModelTrainer(object): pass -# Loss Functions: BCELoss, BCEWithLogitsLoss, CrossEntropyLoss (not working), MSELoss (with reduction='sum') +# Loss Functions: BCELoss, BCEWithLogitsLoss, +# CrossEntropyLoss (not working), MSELoss (with reduction='sum') # Optimizer: Adam, SGD trainer = ModelTrainer({'num_epochs': 100}) # trainer.epoch_training() diff --git a/Basis_Function.py b/Basis_Function.py index f0f8046..e5b726e 100644 --- a/Basis_Function.py +++ b/Basis_Function.py @@ -2,6 +2,8 @@ """ @author: Laura C. Kühle +TODO: Contemplate whether calculating projections during initialization can save time + """ import numpy as np from sympy import Symbol, integrate @@ -48,7 +50,8 @@ class Legendre(Vector): if degree == 1: vector.append(eval_point) else: - poly = (2.0*degree - 1)/degree * eval_point * vector[-1] - (degree-1)/degree * vector[-2] + poly = (2.0*degree - 1)/degree * eval_point * vector[-1] \ + - (degree-1)/degree * vector[-2] vector.append(poly) return vector @@ -56,7 +59,8 @@ class Legendre(Vector): class OrthonormalLegendre(Legendre): def _build_basis_vector(self, eval_point): leg_vector = self._calculate_legendre_vector(eval_point) - return [leg_vector[degree] * np.sqrt(degree+0.5) for degree in range(self._polynomial_degree+1)] + return [leg_vector[degree] * np.sqrt(degree+0.5) + for degree in range(self._polynomial_degree+1)] def _build_wavelet_vector(self, eval_point): degree = self._polynomial_degree @@ -71,9 +75,12 @@ class OrthonormalLegendre(Legendre): 1/3 * np.sqrt(2.5) * (4 - 15*eval_point + 12*(eval_point**2))] if degree == 3: return [np.sqrt(15/34) * (1 + 4*eval_point - 30*(eval_point**2) + 28*(eval_point**3)), - np.sqrt(1/42) * (-4 + 105*eval_point - 300*(eval_point**2) + 210*(eval_point**3)), - 1/2 * np.sqrt(35/34) * (-5 + 48*eval_point - 105*(eval_point**2) + 64*(eval_point**3)), - 1/2 * np.sqrt(5/42) * (-16 + 105*eval_point - 192*(eval_point**2) + 105*(eval_point**3))] + np.sqrt(1/42) * (-4 + 105*eval_point - 300*(eval_point**2) + + 210*(eval_point**3)), + 1/2 * np.sqrt(35/34) * (-5 + 48*eval_point - 105*(eval_point**2) + + 64*(eval_point**3)), + 1/2 * np.sqrt(5/42) * (-16 + 105*eval_point - 192*(eval_point**2) + + 105*(eval_point**3))] if degree == 4: return [np.sqrt(1/186) * (1 + 30*eval_point + 210*(eval_point**2) - 840*(eval_point**3) + 630*(eval_point**4)), @@ -116,7 +123,8 @@ class OrthonormalLegendre(Legendre): for i in range(self._polynomial_degree+1): row = [] for j in range(self._polynomial_degree+1): - entry = integrate(self._basis[i].subs(x, first_param) * self._wavelet[j].subs(z, second_param), + entry = integrate(self._basis[i].subs(x, first_param) + * self._wavelet[j].subs(z, second_param), (z, -1, 1)) if is_left_matrix: entry = entry * (-1)**(j + self._polynomial_degree + 1) diff --git a/DG_Approximation.py b/DG_Approximation.py index d0d1bd5..925aa6a 100644 --- a/DG_Approximation.py +++ b/DG_Approximation.py @@ -10,6 +10,7 @@ TODO: Check whether 'projection' is always a np.array() TODO: Check whether all instance variables sensible TODO: Use cfl_number for updating, not just time TODO: Adjust code to allow classes for all equations (Burger, linear advection, 1D Euler) +TODO: Change maximal line length to 100 -> Done """ import os @@ -78,15 +79,16 @@ class DGScheme(object): # Replace the string names with the actual class instances # (and add the instance variables for the quadrature) - self._init_cond = getattr(Initial_Condition, self._init_cond)(self._left_bound, self._right_bound, - self._init_config) + self._init_cond = getattr(Initial_Condition, self._init_cond)( + self._left_bound, self._right_bound, self._init_config) self._limiter = getattr(Limiter, self._limiter)(self._limiter_config) self._quadrature = getattr(Quadrature, self._quadrature)(self._quadrature_config) self._detector = getattr(Troubled_Cell_Detector, self._detector)( - self._detector_config, self._mesh, self._wave_speed, self._polynomial_degree, self._num_grid_cells, - self._final_time, self._left_bound, self._right_bound, self._basis, self._init_cond, self._quadrature) - self._update_scheme = getattr(Update_Scheme, self._update_scheme)(self._polynomial_degree, self._num_grid_cells, - self._detector, self._limiter) + self._detector_config, self._mesh, self._wave_speed, self._polynomial_degree, + self._num_grid_cells, self._final_time, self._left_bound, self._right_bound, + self._basis, self._init_cond, self._quadrature) + self._update_scheme = getattr(Update_Scheme, self._update_scheme)( + self._polynomial_degree, self._num_grid_cells, self._detector, self._limiter) def approximate(self): """ @@ -127,10 +129,11 @@ class DGScheme(object): plt.show() def save_plots(self): - name = self._init_cond.get_name() + '__' + self._detector.get_name() + '__' + self._limiter.get_name() \ - + '__' + self._update_scheme.get_name() + '__' + self._quadrature.get_name() + '__final_time_' \ - + str(self._final_time) + '__wave_speed_' + str(self._wave_speed) + '__number_of_cells_' \ - + str(self._num_grid_cells) + '__polynomial_degree_' + str(self._polynomial_degree) + name = self._init_cond.get_name() + '__' + self._detector.get_name() + '__' \ + + self._limiter.get_name() + '__' + self._update_scheme.get_name() + '__' \ + + self._quadrature.get_name() + '__final_time_' + str(self._final_time) \ + + '__wave_speed_' + str(self._wave_speed) + '__number_of_cells_' \ + + str(self._num_grid_cells) + '__polynomial_degree_' + str(self._polynomial_degree) # Set paths for plot files if not existing already if not os.path.exists(self._plot_dir): @@ -155,7 +158,8 @@ class DGScheme(object): self._limiter_config['cell_len'] = self._cell_len # Set mesh with one ghost point on each side - self._mesh = np.arange(self._left_bound - (3/2*self._cell_len), self._right_bound + (5/2*self._cell_len), + self._mesh = np.arange(self._left_bound - (3/2*self._cell_len), + self._right_bound + (5/2*self._cell_len), self._cell_len) # +3/2 # Set inverse mass matrix @@ -181,8 +185,10 @@ class DGScheme(object): for degree in range(self._polynomial_degree + 1): new_entry = sum( - initial_condition.calculate( - eval_point + self._cell_len/2 * self._quadrature.get_eval_points()[point] - adjustment) + initial_condition.calculate(eval_point + + self._cell_len/2 + * self._quadrature.get_eval_points()[point] + - adjustment) * basis_vector[degree].subs(x, self._quadrature.get_eval_points()[point]) * self._quadrature.get_weights()[point] for point in range(self._quadrature.get_num_points())) @@ -203,4 +209,5 @@ class DGScheme(object): initial_condition = self._init_cond projection = self._do_initial_projection(initial_condition, adjustment) - return self._detector.calculate_cell_average_and_reconstructions(projection[:, 1:-1], stencil_length) + return self._detector.calculate_cell_average_and_reconstructions(projection[:, 1:-1], + stencil_length) diff --git a/Initial_Condition.py b/Initial_Condition.py index 33ec6ab..28ef065 100644 --- a/Initial_Condition.py +++ b/Initial_Condition.py @@ -226,7 +226,8 @@ class HeavisideTwoSided(InitialCondition): left_factor = config.pop('left_factor', np.random.choice([-1, 1])) right_factor = config.pop('right_factor', np.random.choice([-1, 1])) adjustment = config.pop('adjustment', np.random.uniform(low=-1, high=1)) - config = {'left_factor': left_factor, 'right_factor': right_factor, 'adjustment': adjustment} + config = {'left_factor': left_factor, 'right_factor': right_factor, + 'adjustment': adjustment} self._reset(config) def _get_point(self, x): diff --git a/Main.py b/Main.py index 5ae4083..3e96072 100644 --- a/Main.py +++ b/Main.py @@ -14,7 +14,7 @@ tic = timeit.default_timer() alpha = 1 p = 2 cfl = 0.2 -N = 4 # 40 elements work well for Condition 3 +N = 32 # 40 elements work well for Condition 3 finalTime = 1 xL = -1 xR = 1 @@ -44,10 +44,11 @@ update_scheme = 'SSPRK3' init_cond = 'Sine' limiter = 'ModifiedMinMod' quadrature = 'Gauss' -dg_scheme = DGScheme(detector, detector_config=detector_config, init_cond=init_cond, init_config=init_config, - limiter=limiter, limiter_config=limiter_config, quadrature=quadrature, - quadrature_config=quadrature_config, update_scheme=update_scheme, wave_speed=alpha, - polynomial_degree=p, cfl_number=cfl, num_grid_cells=N, final_time=finalTime, +dg_scheme = DGScheme(detector, detector_config=detector_config, init_cond=init_cond, + init_config=init_config, limiter=limiter, limiter_config=limiter_config, + quadrature=quadrature, quadrature_config=quadrature_config, + update_scheme=update_scheme, wave_speed=alpha, polynomial_degree=p, + cfl_number=cfl, num_grid_cells=N, final_time=finalTime, verbose=False, left_bound=xL, right_bound=xR) dg_scheme.approximate() diff --git a/Plotting.py b/Plotting.py index 43a1663..ec4c888 100644 --- a/Plotting.py +++ b/Plotting.py @@ -94,7 +94,8 @@ def calculate_approximate_solution(projection, points, polynomial_degree, basis) return np.reshape(np.array(approx), (1, len(approx) * num_points)) -def calculate_exact_solution(mesh, cell_len, wave_speed, final_time, interval_len, quadrature, init_cond): +def calculate_exact_solution(mesh, cell_len, wave_speed, final_time, interval_len, quadrature, + init_cond): grid = [] exact = [] num_periods = np.floor(wave_speed * final_time / interval_len) @@ -104,7 +105,8 @@ def calculate_exact_solution(mesh, cell_len, wave_speed, final_time, interval_le eval_values = [] for point in range(len(eval_points)): - new_entry = init_cond.calculate(eval_points[point] - wave_speed * final_time + num_periods * interval_len) + new_entry = init_cond.calculate(eval_points[point] - wave_speed * final_time + + num_periods * interval_len) eval_values.append(new_entry) grid.append(eval_points) diff --git a/Troubled_Cell_Detector.py b/Troubled_Cell_Detector.py index b73e1f9..9bd0dcd 100644 --- a/Troubled_Cell_Detector.py +++ b/Troubled_Cell_Detector.py @@ -13,16 +13,16 @@ import torch from sympy import Symbol import ANN_Model -from Plotting import plot_solution_and_approx, plot_semilog_error, plot_error, plot_shock_tube, plot_details, \ - calculate_approximate_solution, calculate_exact_solution +from Plotting import plot_solution_and_approx, plot_semilog_error, plot_error, plot_shock_tube, \ + plot_details, calculate_approximate_solution, calculate_exact_solution x = Symbol('x') z = Symbol('z') class TroubledCellDetector(object): - def __init__(self, config, mesh, wave_speed, polynomial_degree, num_grid_cells, final_time, left_bound, right_bound, - basis, init_cond, quadrature): + def __init__(self, config, mesh, wave_speed, polynomial_degree, num_grid_cells, final_time, + left_bound, right_bound, basis, init_cond, quadrature): self._mesh = mesh self._wave_speed = wave_speed self._polynomial_degree = polynomial_degree @@ -57,16 +57,17 @@ class TroubledCellDetector(object): def calculate_cell_average_and_reconstructions(self, projection, stencil_length): """ - Calculate the cell averages of all cells in a projection. Reconstructions are only calculated for the middle - cell and added left and right to it, respectively. + Calculate the cell averages of all cells in a projection. Reconstructions are only + calculated for the middle cell and added left and right to it, respectively. Here come some parameter. """ - cell_averages = calculate_approximate_solution(projection, [0], 0, self._basis.get_basis_vector()) - left_reconstructions = calculate_approximate_solution(projection, [-1], self._polynomial_degree, - self._basis.get_basis_vector()) - right_reconstructions = calculate_approximate_solution(projection, [1], self._polynomial_degree, - self._basis.get_basis_vector()) + cell_averages = calculate_approximate_solution( + projection, [0], 0, self._basis.get_basis_vector()) + left_reconstructions = calculate_approximate_solution( + projection, [-1], self._polynomial_degree, self._basis.get_basis_vector()) + right_reconstructions = calculate_approximate_solution( + projection, [1], self._polynomial_degree, self._basis.get_basis_vector()) middle_idx = stencil_length//2 return np.array(list(map(np.float64, zip(cell_averages[:, :middle_idx], left_reconstructions[:, middle_idx], cell_averages[:, middle_idx], @@ -81,10 +82,12 @@ class TroubledCellDetector(object): print('maximum error =', max_error) def _plot_mesh(self, projection): - grid, exact = calculate_exact_solution(self._mesh[2:-2], self._cell_len, self._wave_speed, self._final_time, - self._interval_len, self._quadrature, self._init_cond) - approx = calculate_approximate_solution(projection[:, 1:-1], self._quadrature.get_eval_points(), - self._polynomial_degree, self._basis.get_basis_vector()) + grid, exact = calculate_exact_solution( + self._mesh[2:-2], self._cell_len, self._wave_speed, self._final_time, + self._interval_len, self._quadrature, self._init_cond) + approx = calculate_approximate_solution( + projection[:, 1:-1], self._quadrature.get_eval_points(), self._polynomial_degree, + self._basis.get_basis_vector()) pointwise_error = np.abs(exact-approx) max_error = np.max(pointwise_error) @@ -108,11 +111,11 @@ class ArtificialNeuralNetwork(TroubledCellDetector): self._stencil_len = config.pop('stencil_len', 3) self._model = config.pop('model', 'ThreeLayerReLu') - self._model_config = config.pop('model_config', {'input_size': self._stencil_len+2, 'first_hidden_size': 8, - 'second_hidden_size': 4, 'output_size': 2, - 'activation_function': 'Softmax', - 'activation_config': {'dim': 1}}) - self._model_state = config.pop('model_state', 'Train24k24k_Valid8k8k_Norm12ReLU8+4nodesSM1Adamlr1e-2MSE.pt') + self._model_config = config.pop('model_config', { + 'input_size': self._stencil_len+2, 'first_hidden_size': 8, 'second_hidden_size': 4, + 'output_size': 2, 'activation_function': 'Softmax', 'activation_config': {'dim': 1}}) + self._model_state = config.pop( + 'model_state', 'Train24k24k_Valid8k8k_Norm12ReLU8+4nodesSM1Adamlr1e-2MSE.pt') if not hasattr(ANN_Model, self._model): raise ValueError('Invalid model: "%s"' % self._model) @@ -122,7 +125,8 @@ class ArtificialNeuralNetwork(TroubledCellDetector): # Reset ghost cells to adjust for stencil length num_ghost_cells = self._stencil_len//2 projection = projection[:, 1:-1] - projection = np.concatenate((projection[:, -num_ghost_cells:], projection, projection[:, :num_ghost_cells]), + projection = np.concatenate((projection[:, -num_ghost_cells:], projection, + projection[:, :num_ghost_cells]), axis=1) # Calculate input data depending on stencil length @@ -136,7 +140,8 @@ class ArtificialNeuralNetwork(TroubledCellDetector): # Return troubled cells model_output = torch.round(self._model(input_data.float())) - return [cell for cell in range(len(model_output)) if model_output[cell, 0] == torch.tensor([1])] + return [cell for cell in range(len(model_output)) + if model_output[cell, 0] == torch.tensor([1])] class WaveletDetector(TroubledCellDetector): @@ -151,7 +156,8 @@ class WaveletDetector(TroubledCellDetector): # Set additional necessary parameter self._num_coarse_grid_cells = self._num_grid_cells//2 - self._wavelet_projection_left, self._wavelet_projection_right = self._basis.get_multiwavelet_projections() + self._wavelet_projection_left, self._wavelet_projection_right \ + = self._basis.get_multiwavelet_projections() def get_cells(self, projection): multiwavelet_coeffs = self._calculate_wavelet_coeffs(projection[:, 1: -1]) @@ -171,8 +177,9 @@ class WaveletDetector(TroubledCellDetector): def plot_results(self, projection, troubled_cell_history, time_history): multiwavelet_coeffs = self._calculate_wavelet_coeffs(projection) coarse_projection = self._calculate_coarse_projection(projection) - plot_details(projection[:, 1:-1], self._mesh[2:-2], coarse_projection, self._basis.get_basis_vector(), - self._basis.get_wavelet_vector(), multiwavelet_coeffs, self._num_coarse_grid_cells, + plot_details(projection[:, 1:-1], self._mesh[2:-2], coarse_projection, + self._basis.get_basis_vector(), self._basis.get_wavelet_vector(), + multiwavelet_coeffs, self._num_coarse_grid_cells, self._polynomial_degree) super().plot_results(projection, troubled_cell_history, time_history) @@ -193,16 +200,19 @@ class WaveletDetector(TroubledCellDetector): return coarse_projection def _plot_mesh(self, projection): - grid, exact = calculate_exact_solution(self._mesh[2:-2], self._cell_len, self._wave_speed, self._final_time, - self._interval_len, self._quadrature, self._init_cond) - approx = calculate_approximate_solution(projection[:, 1:-1], self._quadrature.get_eval_points(), - self._polynomial_degree, self._basis.get_basis_vector()) + grid, exact = calculate_exact_solution( + self._mesh[2:-2], self._cell_len, self._wave_speed, self._final_time, + self._interval_len, self._quadrature, self._init_cond) + approx = calculate_approximate_solution( + projection[:, 1:-1], self._quadrature.get_eval_points(), self._polynomial_degree, + self._basis.get_basis_vector()) pointwise_error = np.abs(exact-approx) max_error = np.max(pointwise_error) self._plot_coarse_mesh(projection) - plot_solution_and_approx(grid, exact, approx, self._colors['fine_exact'], self._colors['fine_approx']) + plot_solution_and_approx(grid, exact, approx, self._colors['fine_exact'], + self._colors['fine_approx']) plt.legend(['Exact (Coarse)', 'Approx (Coarse)', 'Exact (Fine)', 'Approx (Fine)']) plot_semilog_error(grid, pointwise_error) plot_error(grid, exact, approx) @@ -211,17 +221,21 @@ class WaveletDetector(TroubledCellDetector): def _plot_coarse_mesh(self, projection): coarse_cell_len = 2*self._cell_len - coarse_mesh = np.arange(self._left_bound - (0.5*coarse_cell_len), self._right_bound + (1.5*coarse_cell_len), + coarse_mesh = np.arange(self._left_bound - (0.5*coarse_cell_len), + self._right_bound + (1.5*coarse_cell_len), coarse_cell_len) coarse_projection = self._calculate_coarse_projection(projection) # Plot exact and approximate solutions for coarse mesh - grid, exact = calculate_exact_solution(coarse_mesh[1:-1], coarse_cell_len, self._wave_speed, self._final_time, - self._interval_len, self._quadrature, self._init_cond) - approx = calculate_approximate_solution(coarse_projection, self._quadrature.get_eval_points(), - self._polynomial_degree, self._basis.get_basis_vector()) - plot_solution_and_approx(grid, exact, approx, self._colors['coarse_exact'], self._colors['coarse_approx']) + grid, exact = calculate_exact_solution( + coarse_mesh[1:-1], coarse_cell_len, self._wave_speed, self._final_time, + self._interval_len, self._quadrature, self._init_cond) + approx = calculate_approximate_solution( + coarse_projection, self._quadrature.get_eval_points(), self._polynomial_degree, + self._basis.get_basis_vector()) + plot_solution_and_approx( + grid, exact, approx, self._colors['coarse_exact'], self._colors['coarse_approx']) class Boxplot(WaveletDetector): @@ -282,7 +296,8 @@ class Theoretical(WaveletDetector): def _get_cells(self, multiwavelet_coeffs, projection): troubled_cells = [] - max_avg = np.sqrt(0.5) * max(1, max(abs(projection[0][cell+1]) for cell in range(self._num_coarse_grid_cells))) + max_avg = np.sqrt(0.5) * max(1, max(abs(projection[0][cell+1]) + for cell in range(self._num_coarse_grid_cells))) for cell in range(self._num_coarse_grid_cells): if self._is_troubled_cell(multiwavelet_coeffs, cell, max_avg): diff --git a/Update_Scheme.py b/Update_Scheme.py index d21c09f..4df2451 100644 --- a/Update_Scheme.py +++ b/Update_Scheme.py @@ -75,11 +75,13 @@ class SSPRK3(UpdateScheme): current_projection, __ = self._apply_limiter(current_projection) current_projection = self._enforce_boundary_condition(current_projection) - current_projection = self._apply_second_step(original_projection, current_projection, cfl_number) + current_projection = self._apply_second_step(original_projection, current_projection, + cfl_number) current_projection, __ = self._apply_limiter(current_projection) current_projection = self._enforce_boundary_condition(current_projection) - current_projection = self._apply_third_step(original_projection, current_projection, cfl_number) + current_projection = self._apply_third_step(original_projection, current_projection, + cfl_number) current_projection, troubled_cells = self._apply_limiter(current_projection) current_projection = self._enforce_boundary_condition(current_projection) -- GitLab