From 1f5ead3bf1a0b211aa25afb35ec781016316a23d Mon Sep 17 00:00:00 2001 From: Aadi B <122440845+aadi-bh@users.noreply.github.com> Date: Fri, 24 Nov 2023 14:43:22 +0530 Subject: [PATCH] Restructure --- common.py | 19 +++++++- filters.py | 18 ++++++++ ic.py | 16 +++++++ main.py | 124 +++++++++++++++++++++++++++++++++++++++++++++++++++++ op.py | 59 +++++++++++++++++++++++++ 5 files changed, 235 insertions(+), 1 deletion(-) create mode 100644 filters.py create mode 100644 ic.py create mode 100644 main.py create mode 100644 op.py diff --git a/common.py b/common.py index 854c282..fe3e779 100644 --- a/common.py +++ b/common.py @@ -1,4 +1,5 @@ import numpy as np +from numpy.fft import * def elrk4(SemiGroup,Nonlinear,y0,tinterval,dt,args): y = y0 @@ -54,8 +55,24 @@ def elrk4(SemiGroup,Nonlinear,y0,tinterval,dt,args): flag = False elif np.any(np.isnan(y)): flag = False -# print("t, dt = ", t, dt) + print("t, dt = ", t, dt) times = np.array(time) solution.append(y) return times, solution + +def create_filter(k, sigma, args): + K = np.max(np.abs(k)) + filter = sigma(k / K, **args) + return filter + +def mask(c, tol=1e-14): + return c * np.where(np.abs(c) np.pi/2, 1, 0) + return y + +def bump(x): + return np.exp(-6 * (x-pi)**2) diff --git a/main.py b/main.py new file mode 100644 index 0000000..fcf868f --- /dev/null +++ b/main.py @@ -0,0 +1,124 @@ +#!/usr/bin/env python +# coding: utf-8 + +import numpy as np +from numpy.fft import * +from scipy.integrate import solve_ivp +from scipy import sparse +from numpy import pi +import matplotlib.pyplot as plt +import argparse + +from common import * +import ic +import op +import filters + + +xmin, xmax = 0, 2 * pi + +parser = argparse.ArgumentParser() +# parser.add_argument('-N', type=int, help='Number of cells', default=100) +parser.add_argument('--cfl', type=float, help='CFL number', default=0.1) +# parser.add_argument('-scheme', +# choices=('C','LF','GLF','LLF','LW','ROE','EROE','GOD'), +# help='Scheme', default='LF') +parser.add_argument('--ic', + choices=('saw_tooth','sin','step', 'bump'), + help='Initial condition', default='sin') +parser.add_argument('--Tf', type=float, help='Final time', default=1.0) +parser.add_argument('--pde', choices=('linadv', 'burgers'), default='burgers') +parser.add_argument('--add_visc', type=bool, default=False) +parser.add_argument('--use_filter', type=bool, default=False) +parser.add_argument('--max_lgN', type=int, default=7) +parser.add_argument('--integrator', choices=('solve_ivp', 'elrk4'), default='elrk4') +args = parser.parse_args() + +if __name__ == '__main__': + rhs = op.linadv + initial_condition = ic.sin + visc = op.semigroup_none + USE_FILTER = False + tf = 2 + cfl = 0.5 + if args.pde == 'burgers': + rhs = op.burgers + if args.ic == 'saw_tooth': + initial_condition = ic.saw_tooth + elif args.ic == 'sin': + initial_condition = ic.sin + elif args.ic == 'step': + initial_condition = ic.step + elif args.ic == 'bump': + initial_condition = ic.bump + if args.add_visc == True: + visc = op.semigroup_heat + tf = max(0, args.Tf) + cfl = min(1, args.cfl) + USE_FILTER = args.use_filter + + print("PDE :", args.pde) + print("TF :", tf) + print("CFL :", cfl) + print("IC :", args.ic) + print("FILTER:", args.use_filter) + print("VISC :", args.add_visc) + + sols = [] + for i in range(0, args.max_lgN - 4 + 1): + N = np.power(2, i + 4); + print(f"N={N}") + M = 3 * N // 2; + m = M // 2; + NN = (2 * m) + N + dx = (xmax-xmin) / N + dt = (cfl * dx) + x = np.arange(xmin, xmax, dx) + k = freqs(N) + kk = freqs(NN) + filter = np.ones(len(kk)) + p = i + if USE_FILTER == True: + filter = create_filter(kk, sigma, {'p':p}) + args = (N, M, filter) + u_hat_init = fft(initial_condition(x)) + S_half, S = visc(dt, k, eps = 1e-2) + ''' + output = solve_ivp(fun = rhs, + t_span = [0, tf], + t_eval = [0, tf], + y0 = u_hat_init, + args = args) + times, u = output.t, output.y.transpose() + ''' + times, u = elrk4([S_half, S], rhs, u_hat_init, (0, tf), dt, args) + sols.append((times, u)) + + filename = f"{N}.txt" + np.savetxt(filename, np.vstack((k, u[-1].real))) + print("Saved solution to " + filename) + + # PLOT + fig, ax = plt.subplots(nrows=1, ncols=2, figsize = (14, 8), width_ratios=[3,1]) + if rhs == op.burgers and initial_condition == ic.sin: + if tf >= 1 or True: + god = np.loadtxt("burg3_GOD_5.txt").transpose() + elif tf == 0.6: + god = np.loadtxt("burg3_GOD_3.txt").transpose() + elif tf == 0.2: + god = np.loadtxt("burg3_GOD_1.txt").transpose() + ax[0].plot(god[0] * 2 * np.pi, god[1], 'ko', markersize=0.1, label="Godunov flux") + + fig.tight_layout() + for (times, u) in sols: + ax[0].plot(x, ifft(u[-1]).real, label=str(N)+f", t={np.round(times[-1], 3)}") + plot_resolution(u[-1], ax[1], {'linewidth':0.5, 'markersize':0.5}) + fig, ax = plt.subplots(nrows=1, ncols=2, figsize = (14, 8), width_ratios=[3,1]) + x = np.linspace(xmin, xmax, 1000) + ax[0].plot(x, initial_condition(x), linewidth=0.1, color='k', label="init") + ax[0].legend() + fig.savefig(f"{rhs}-{str(USE_FILTER)}.png") + + plt.close() + print("Done.") + diff --git a/op.py b/op.py new file mode 100644 index 0000000..ab043d9 --- /dev/null +++ b/op.py @@ -0,0 +1,59 @@ +# Operators + +import numpy as np +from numpy.fft import * +from scipy import sparse +from numpy import pi +from common import * + +def pad(c, m): + # ASSUME c is fft(u) / len(u) + N = len(c) + newN = 2*m + N + r = fftshift(c) + r = np.r_[np.zeros(m), r, np.zeros(m)] + r *= newN / N + r = ifftshift(r) + return r + +def unpad(c, m): + N = len(c) + newN = N - 2 * m + r = fftshift(c) + r = r[m:m + newN] + r *= newN / N + r = ifftshift(r) + return r + +def linadv(t, u_hat, N, M, filter, a=1): + NN = (2 * M//2) + N + u_hat = pad(u_hat, M//2) + kk = freqs(NN) + nonlinear = -1j * kk * u_hat + nonlinear *= a + nonlinear *= filter + return unpad(mask(nonlinear), M//2) + +def burgers(t, u_hat, N, M, filter, a=1): + NN = (2 * M//2) + N + u_hat = pad(u_hat, M//2) + u = ifft(u_hat) + kk = freqs(NN) + nonlinear = -1 * fft(u**2/2) * 1j * kk + nonlinear *= a + nonlinear *= filter + return unpad(mask(nonlinear), M//2) + +def semigroup_heat(dt, k, eps): + S_half = sparse.diags(np.exp(-1 * eps * k**2 * dt / 2)) + S = sparse.diags(np.exp(-1 * eps * k**2 * dt)) + return S_half, S + +def semigroup_none(dt, k, eps): + S = sparse.diags(np.ones(len(k))) + S_half = S + return S_half, S + +def zero(t, u_hat, N, M, eps, filter, a=1): + return np.zeros(len(u_hat)); +