Skip to content
Snippets Groups Projects
Commit b744085c authored by Laura Christine Kühle's avatar Laura Christine Kühle
Browse files

Changed maximal line length in code to 100.

parent 55c61934
No related branches found
No related tags found
No related merge requests found
......@@ -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
......
......@@ -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)
......
......@@ -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()
......
......@@ -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)
......
......@@ -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,9 +129,10 @@ 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_' \
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
......@@ -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)
......@@ -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):
......
......@@ -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()
......
......@@ -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)
......
......@@ -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,
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())
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,
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())
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,
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'])
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):
......
......@@ -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)
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment