import time import numpy as np from vertical_discr import * from firedrake import * """ ********************************************* * Definition of the space/time domain * ********************************************* """ #______________________ Beach ______________________# xb = 7.5 # Start of the beach sb = 0.2 # Slope of the beach #______________________ Basin ______________________# g = 9.81 # Gravitational constant H0 = 1.0 # Depth at rest (flat bottom) Hend = 0.5 # Depth at the end of the beach Lx = xb +(H0-Hend)/sb # Length in x Ly = 2.0 # Length in y Lw = 4.0 # End of the x-transform #____________________ Wavemaker ____________________# lamb = 2.0 # Wavelength k = 2*pi/lamb # Wave number w = sqrt(g*k*tanh(k*H0)) # Wave frequency Tw = 2*pi/w # Wave period gamma = 0.03 # Wave amplitude t_stop = 5.0*Tw # When to stop the wavemaker #______________________ Time _______________________# T0 = 0.0 # Initial time t = T0 # Temporal variable Tend = 15*Tw # Final time dt = 0.002 # Time-step scheme ="SE" # Time scheme """ **************************************************** * Definition of the mesh * **************************************************** """ #_________________ Vertical discretization ________________# n_z = 4 # Order of the expansion Nz = n_z+1 # Number of point in one element #________________ Horizontal discretization _______________# res = 0.05 # resolution Nx = int(Lx/res) + 1 # Number of elements in x Ny = int(Ly/res) + 1 # Number of elements in y #___________________________ Mesh _________________________# if Ly==0: #(x,z)-waves hor_mesh = IntervalMesh(Nx,Lx) else: #(x,y,z)-waves hor_mesh = RectangleMesh(Nx,Ny,Lx,Ly,quadrilateral=True) """ ************************************************ * Definition of the function spaces * ************************************************ """ #___________________ For h and psi_1 __________________# V = FunctionSpace(hor_mesh, "CG", 1) #_____________________ For hat_psi ____________________# Vec = VectorFunctionSpace(hor_mesh, "CG", 1,dim=n_z) #________________ Mixed function space ________________# V_mixed = V*Vec # to update h and hat_psi simultaneously """ ****************************************************** * Definition of the functions * ****************************************************** """ if scheme=="SE": #_________ Symplectic-Euler scheme _________# #______________________ At time t^n _____________________# h_n0 = Function(V) # h^n psi_1_n0 = Function(V) # psi_1^n hat_psi_n0 = Function(Vec) # hat_psi^n #________________ At time t^{n+1} and t^* _______________# psi_1_n1 = Function(V) # psi_1^{n+1} w_n1 = Function(V_mixed) h_n1, hat_psi_star = split(w_n1) # h^{n+1}, hat_psi^* hat_psi_n1 = Function(Vec) # hat_psi^{n+1} else: #________________ Stormer-Verlet scheme _______________# #______________________ At time t^n _____________________# h_n0 = Function(V) # h^n psi_1_n0 = Function(V) # psi_1^n hat_psi_n0 = Function(Vec) # hat_psi^n #_______________ At time t^{n+1/2} and t^* ______________# w_half = Function(V_mixed) h_half, hat_psi_star = split(w_half)# h^{n+1/2}, hat_psi^* #____________________ At time t^{n+1} ___________________# psi_1_n1 = Function(V) # psi_1^{n+1} w_n1 = Function(V_mixed) h_n1, hat_psi_n1 = split(w_n1) # h^{n+1}, hat_psi^{n+1} x_coord = Function(V).interpolate(Expression('x[0]')) #___________________________ Beach __________________________# beach = Function(V) # b(x) #_______________________ Depth at rest ______________________# H = Function(V) # H(x) #_________________________ Wavemaker ________________________# WM = Function(V) # R(x,y;t^n) dWM_dt = Function(V) # (dR/dt)^n dWM_dy = Function(V) # (dR/dy)^n WM_n1 = Function(V) # NEED TO EVALUATE WM DIFFERENTLY SELON LE TIME SCHEME #______________________ Trial functions ______________________# psi_1 = TrialFunction(V) # psi_1^{n+1} for linear solvers hat_psi = TrialFunction(Vec) # hat_psi^{n+1} for linear solvers #_______________________ Test functions ______________________# delta_h = TestFunction(V) # from dH/dh delta_hat_psi = TestFunction(Vec) # from dH/dhat_psi w_t = TestFunction(V_mixed) # from dH/dpsi_1... delta_psi, delta_hat_star = split(w_t) # ...and dH/dhat_psi """ ************************************************************************* * Initialisation of the Functions * *************************************************************************""" #______________________________ Beach topography ______________________________# beach_expr = \ Expression("0.5*(1+copysign(1.0,x[0]-xb))*slope*(x[0]-xb)",xb=xb, slope=sb) beach.interpolate(beach_expr) # b(x) #________________________________ Depth at rest _______________________________# #H.assign(H0-beach) # H(x) H.interpolate(Expression("H0-0.5*(1+copysign(1.0,x[0]-xb))*slope*(x[0]-xb)",xb=xb, slope=sb)) dim = 3 if dim==2: #______________________ Wavemaker motion and derivatives ______________________# WM_expr = \ Expression("-0.5*(1+copysign(1.0,Lw-x[0]))*A*cos(w*t)",\ A=gamma, Lw = Lw, w=w, t=t) WM.interpolate(WM_expr) # \tilde{R}(x,y;t) WM_n1_expr = \ Expression("-0.5*(1+copysign(1.0,Lw-x[0]))*A*cos(w*t)",\ A=gamma, Lw = Lw, w=w, t=t+dt) WM_n1.interpolate(WM_n1_expr) dWM_dt_expr = \ Expression("0.5*(1+copysign(1.0,Lw-x[0]))*A*w*sin(w*t)",\ A=gamma,Lw=Lw, w=w, t=t) dWM_dt.interpolate(dWM_dt_expr) # d\tilde{R}/dt dWM_dy_expr = Expression("0.0") dWM_dy.interpolate(dWM_dy_expr) else: #______________________ Wavemaker motion and derivatives ______________________# WM_expr = \ Expression("-0.5*(1+copysign(1.0,Lw-x[0]))*A*(x[1]-0.5*Ly)/(0.5*Ly)*cos(w*t)",\ A=gamma, Ly=Ly, Lw = Lw, w=w, t=t) WM.interpolate(WM_expr) # \tilde{R}(x,y;t) dWM_dt_expr = \ Expression("0.5*(1+copysign(1.0,Lw-x[0]))*A*w*(x[1]-0.5*Ly)/(0.5*Ly)*sin(w*t)",\ A=gamma, Ly=Ly, Lw=Lw, w=w, t=t) dWM_dt.interpolate(dWM_dt_expr) # d\tilde{R}/dt dWM_dy_expr = Expression("-0.5*(1+copysign(1.0,Lw-x[0]))*A*cos(w*t)/(0.5*Ly)",\ A=gamma, Ly=Ly, Lw=Lw, w=w, t=t) dWM_dy.interpolate(dWM_dy_expr) # d\tilde{R}/dy #____________________________________ Depth ___________________________________# h_n0.assign(H) # h(x,y;t=0) = H(x) #________________ Velocity pot. at the surface: phi(x,y,z=h;t) ________________# psi_1_n0.assign(0.0) # \psi_1(x,y;t=0) = 0 #________________ Velocity pot. in depth: phi(x,y,z