Page MenuHomec4science
Files F62928354

ddm_wave_prop.py
No OneTemporary
Actions

File Metadata

Created
Thu, May 16, 14:21

ddm_wave_prop.py
View Options

import numpy as np
import matplotlib.pyplot as plt
import random
import imageio # to create gif
import pandas as pd
# MAIN COMPUTATION
# Fixed node at the top
def wave_propagation_old(dataset, Nnodes, E, A, L, displacement, Ttime, dt):
"""
INPUT
dataset = array of points [epsilon,sigma]
Nnodes = number of nodes in the 1D mesh
E = Young modulus
A = beam section
L = beam length
displacement = imposed displacement at bottom node 0 (function of time)
Ttime = total time of simulation
dt = time step for simumation
OUTPUT
coords = nodes coordinates
all_displacements = nodes displacements at each timestep
KE_list = kinetic energy at each timestep
SE_list = strain energy at each timestep
TE_list = total energy at each timestep
"""
# Dataset properties
Ndata = dataset.shape[0]
# Distance function
Ce = 1
We = lambda x : 0.5*Ce*x**2
We_s = lambda x : 0.5/Ce*x**2
def F_cost(epsis, sigmas, weights):
return np.multiply(weights, We(epsis)+We_s(sigmas))
# Geometry
coords = L*np.arange(Nnodes)
Nel = Nnodes - 1
# Time parameters
Nt = int(Ttime/dt)
beta = 0.25
gamma = 0.5
# Generate matrices
M, Ktot = compute_M_K(Nnodes, E, A, L)
B, weights = compute_B_weights_old(Nnodes, coords)
Keq = compute_Keq(Ce, B, weights)
# Apply boundary conditions
Nfree = Nnodes - 2
F_free = np.zeros(Nfree)
Keq_free = Keq[1:-1,1:-1]
M_free = M[1:-1,1:-1]
B_free = B[:,1:-1]
BigMat = compute_BigMat(Keq_free, M_free, beta, dt)
virtual_force = np.hstack((Keq[1:-1,0],M[1:-1,0]/(beta*dt**2)))
# Pick initial guess
random.seed(42)
guess_list = [random.randint(0,Ndata-1) for _ in range(Nel)]
initial_guess = dataset[guess_list,:]
# Initialize displacement, velocity and acceleration vectors
u, v, a = np.zeros(Nnodes), np.zeros(Nnodes), np.zeros(Nnodes)
eta = np.zeros(Nnodes)
converged = False
all_displacements = np.zeros(Nnodes)
guess = initial_guess
# Energy (for plotting only)
KE_list = []
SE_list = []
TE_list = []
for k in range(Nt): # Loop for time integration
u_previous, v_previous , a_previous = u.copy(), v.copy(), a.copy()
# Predictor step
u_pred = u_previous + dt*v_previous + (0.5-beta)*dt**2*a_previous
v_pred = v_previous + (1-gamma)*dt*a_previous
# Solving corrector step with fixed point:
for i in range(100): # Loop for fixed point
# Imposed displacement
disp = displacement(k*dt)
rhs_up = Ce*np.einsum('e,e,ei->i',weights,guess[:,0],B_free)
rhs_down = F_free.ravel() - np.einsum('e,e,ei->i',weights,guess[:,1],B_free) + M_free.dot(u_pred[1:-1])/(beta*dt**2)
rhs = np.vstack((rhs_up.reshape((-1,1)), rhs_down.reshape((-1,1)))).ravel() - disp*virtual_force
u_eta = np.linalg.solve(BigMat, rhs)
u[1:-1]=u_eta[:Nfree]
u[0] = disp
eta[1:-1] = u_eta[Nfree:]
epsi_comp = B.dot(u)
sigma_comp = guess[:,1] + Ce*B.dot(eta)
closest = find_closest(epsi_comp, sigma_comp, weights, dataset, F_cost)
if closest == guess_list:
converged = True
break
guess_list=closest
guess=dataset[guess_list,:]
if converged:
all_displacements = np.vstack([all_displacements,u])
a = (u-u_pred)/(beta*dt**2)
v = v_pred+gamma*dt*a
KE = v.T @ M @ v
SE = u.T @ Ktot @ u
KE_list.append(KE)
SE_list.append(SE)
TE_list.append(KE+SE)
else: break
print('DD computation finished')
return coords, all_displacements, KE_list, SE_list, TE_list
# Free node at the top
def wave_propagation(dataset, Nnodes, E, A, L, displacement, Ttime, dt):
"""
INPUT
dataset = array of points [epsilon,sigma]
Nnodes = number of nodes in the 1D mesh
E = Young modulus
A = beam section
L = beam length
displacement = imposed displacement at bottom node 0 (function of time)
Ttime = total time of simulation
dt = time step for simumation
OUTPUT
coords = nodes coordinates
all_displacements = nodes displacements at each timestep
KE_list = kinetic energy at each timestep
SE_list = strain energy at each timestep
TE_list = total energy at each timestep
"""
# Dataset properties
Ndata = dataset.shape[0]
# Distance function
Ce = 1
We = lambda x : 0.5*Ce*x**2
We_s = lambda x : 0.5/Ce*x**2
def F_cost(epsis, sigmas, weights):
return np.multiply(weights, We(epsis)+We_s(sigmas))
# Geometry
coords = L*np.arange(Nnodes)
Nel = Nnodes - 1
# Time parameters
Nt = int(Ttime/dt)
beta = 0.25
gamma = 0.5
# Generate matrices
M, Ktot = compute_M_K(Nnodes, E, A, L) # Ktot only used for energy verification
B, weights = compute_B_weights_old(Nnodes, coords)
Keq = compute_Keq(Ce, B, weights)
# Apply boundary conditions
Nfree = Nnodes - 1
F_free = np.zeros(Nfree)
Keq_free = Keq[1:,1:]
M_free = M[1:,1:]
B_free = B[:,1:]
BigMat = compute_BigMat(Keq_free, M_free, beta, dt)
virtual_force = np.hstack((Keq[1:,0],M[1:,0]/(beta*dt**2)))
# Pick initial guess
random.seed(42)
guess_list = [random.randint(0,Ndata-1) for _ in range(Nel)]
initial_guess = dataset[guess_list,:]
# Initialize displacement, velocity and acceleration vectors
u, v, a = np.zeros(Nnodes), np.zeros(Nnodes), np.zeros(Nnodes)
eta = np.zeros(Nnodes)
converged = False
all_displacements = np.zeros(Nnodes)
guess = initial_guess
# Energy (for plotting only)
KE_list = []
SE_list = []
TE_list = []
for k in range(Nt): # Loop for time integration
u_previous, v_previous , a_previous = u.copy(), v.copy(), a.copy()
# Predictor step
u_pred = u_previous + dt*v_previous + (0.5-beta)*dt**2*a_previous
v_pred = v_previous + (1-gamma)*dt*a_previous
# Solving corrector step with fixed point:
for i in range(100): # Loop for fixed point
# Imposed displacement
disp = displacement(k*dt)
rhs_up = Ce*np.einsum('e,e,ei->i',weights,guess[:,0],B_free)
rhs_down = F_free.ravel() - np.einsum('e,e,ei->i',weights,guess[:,1],B_free) + M_free.dot(u_pred[1:])/(beta*dt**2)
rhs = np.vstack((rhs_up.reshape((-1,1)), rhs_down.reshape((-1,1)))).ravel() - disp*virtual_force
u_eta = np.linalg.solve(BigMat, rhs)
u[1:]=u_eta[:Nfree]
u[0] = disp
eta[1:] = u_eta[Nfree:]
epsi_comp = B.dot(u)
sigma_comp = guess[:,1] + Ce*B.dot(eta)
closest = find_closest(epsi_comp, sigma_comp, weights, dataset, F_cost)
if closest == guess_list:
converged = True
break
guess_list=closest
guess=dataset[guess_list,:]
if converged:
all_displacements = np.vstack([all_displacements,u])
a = (u-u_pred)/(beta*dt**2)
v = v_pred+gamma*dt*a
KE = v.T @ M @ v
SE = u.T @ Ktot @ u
KE_list.append(KE)
SE_list.append(SE)
TE_list.append(KE+SE)
else: break
print('DD computation finished')
return coords, all_displacements, KE_list, SE_list, TE_list
# Free node at the top, confinement pressure
def wave_propagation_pressure(dataset, E, A, L, displacement, Ttime, dt, rho=1):
"""
INPUT
dataset = array of points [epsilon,sigma]
Nnodes = number of nodes in the 1D mesh
E = Young modulus
A = beam section
L = beam length
displacement = imposed displacement at bottom node 0 (function of time)
Ttime = total time of simulation
dt = time step for simumation
OUTPUT
coords = nodes coordinates
all_displacements = nodes displacements at each timestep
KE_list = kinetic energy at each timestep
SE_list = strain energy at each timestep
TE_list = total energy at each timestep
"""
# Dataset properties
pressures = [p for p in dataset.keys()]
pressures.sort(reverse=True) # all pressures from largest to smallest
Nel = len(pressures)
Nnodes = Nel + 1
# Distance function
Ce = E
We = lambda x : 0.5*Ce*x**2
We_s = lambda x : 0.5/Ce*x**2
def F_cost(epsis, sigmas, weights):
return np.multiply(weights, We(epsis)+We_s(sigmas))
# Geometry
coords = L*np.arange(Nnodes)
# Time parameters
Nt = int(Ttime/dt)
beta = 0.25
gamma = 0.5
# Generate matrices
M, Ktot = compute_M_K(Nnodes, E, A, L, rho) # Ktot only used for energy verification
B, weights = compute_B_weights(Nnodes, A, L)
Keq = compute_Keq(Ce, B, weights)
# Apply boundary conditions
Nfree = Nnodes - 1
F_free = np.zeros(Nfree)
Keq_free = Keq[1:,1:]
M_free = M[1:,1:]
B_free = B[:,1:]
BigMat = compute_BigMat(Keq_free, M_free, beta, dt)
virtual_force = np.hstack((Keq[1:,0],M[1:,0]/(beta*dt**2)))
# Pick initial guess
initial_guess = np.zeros((Nel,2))
guess_indices = np.zeros(Nel,dtype=int)
for i, p in enumerate(pressures):
Ndata = dataset[p].shape[0]
guess_index = random.randint(0,Ndata-1)
guess_indices[i] = guess_index
initial_guess[i,:] = dataset[p][guess_index,:]
# Initialize displacement, velocity and acceleration vectors
u, v, a = np.zeros(Nnodes), np.zeros(Nnodes), np.zeros(Nnodes)
eta = np.zeros(Nnodes)
converged = False
all_displacements = np.zeros(Nnodes)
guess = initial_guess
guess_indices_old = np.copy(guess_indices)
# Energy (for plotting only)
KE_list = []
SE_list = []
TE_list = []
for k in range(Nt): # Loop for time integration
converged = False # NEW BUT NECESSARY
u_previous, v_previous , a_previous = u.copy(), v.copy(), a.copy()
# Predictor step
u_pred = u_previous + dt*v_previous + (0.5-beta)*dt**2*a_previous
v_pred = v_previous + (1-gamma)*dt*a_previous
# Solving corrector step with fixed point:
for i in range(100): # Loop for fixed point
# Imposed displacement
disp = displacement(k*dt)
rhs_up = Ce*np.einsum('e,e,ei->i',weights,guess[:,0],B_free)
rhs_down = F_free.ravel() - np.einsum('e,e,ei->i',weights,guess[:,1],B_free) + M_free.dot(u_pred[1:])/(beta*dt**2)
rhs = np.vstack((rhs_up.reshape((-1,1)), rhs_down.reshape((-1,1)))).ravel() - disp*virtual_force
u_eta = np.linalg.solve(BigMat, rhs)
u[1:]=u_eta[:Nfree]
u[0] = disp
eta[1:] = u_eta[Nfree:]
epsi_comp = B.dot(u)
sigma_comp = guess[:,1] + Ce*B.dot(eta)
for j, p in enumerate(pressures):
closest_index = find_closest(epsi_comp[j], sigma_comp[j], weights[j], dataset[p], F_cost)
guess_indices[j] = closest_index
if (guess_indices_old == guess_indices).all():
converged = True
break
guess_indices_old = np.copy(guess_indices)
for j, p in enumerate(pressures):
guess[j] = dataset[p][guess_indices[j],:]
if converged:
all_displacements = np.vstack([all_displacements,u])
a = (u-u_pred)/(beta*dt**2)
v = v_pred+gamma*dt*a
KE = v.T @ M @ v
SE = u.T @ Ktot @ u
KE_list.append(KE)
SE_list.append(SE)
TE_list.append(KE+SE)
else:
print('Convergence failed!')
break
print('DD computation finished')
return coords, all_displacements, KE_list, SE_list, TE_list
# LOAD DATASET
def load_dataset(path):
"""
From pickle file, load dataset in the form of a dictionnary (key = confinement pressure)
"""
df = pd.read_pickle(path + 'dataset')
col = df.columns
dataset = {}
for i in range(int(len(col)/2)):
Pressure = int(col[2*i].split()[1][2:-2])
Strain = df[col[2*i]].dropna()
Stress = df[col[2*i+1]].dropna()
dataset[Pressure] = np.array([[e,s] for e,s in zip(Strain,Stress)])
return dataset
# PLOTS
def plot_energy(KE, SE, TE, N0=0):
plt.figure()
plt.plot(KE, label='Kinetic energy')
plt.plot(SE, label='Strain energy')
plt.plot(TE, label='Total energy')
plt.vlines(N0,np.amin(TE),np.amax(TE),colors='r', linestyles='dashed', label='Pulse end')
plt.legend()
plt.show()
def plot_disp(absolute_path,coords,disp,xmin,xmax,num=0):
Nnodes = coords.shape[0]
Nel = Nnodes-1
fig = plt.figure(figsize=(5, 10))
plt.xlim(xmin,xmax)
for bar in range(Nel):
plt.plot([disp[bar],disp[bar+1]],[coords[bar],coords[bar+1]],'ko-')
plt.savefig(absolute_path+'{0}.png'.format(str(num).zfill(3)))
plt.close(fig)
def save_gif(absolute_path, coords, all_displacements):
disp_min = np.amin(all_displacements)
disp_max = np.amax(all_displacements)
nbtot = all_displacements.shape[0]
with imageio.get_writer('displacement.gif', mode='I') as writer:
for k,disp in enumerate(all_displacements):
if k%20 == 0: print('Plotting displacement {nb} out of {tot}'.format(nb=k,tot=nbtot))
plot_disp(absolute_path,coords,disp,disp_min,disp_max,k)
filename = absolute_path+'{0}.png'.format(str(k).zfill(3))
image = imageio.imread(filename)
writer.append_data(image)
def plot_disp_T(absolute_path,coords,disp,xmin,xmax,num=0):
"""
disp is a dictionnary where the keys are the sinus periods
"""
Nnodes = coords.shape[0]
Nel = Nnodes-1
periods = [i for i in disp.keys()]
fig = plt.figure(figsize=(5, 10))
plt.xlim(xmin,xmax)
color = plt.cm.rainbow(np.linspace(0, 1, len(periods)))
for i, T in enumerate(periods):
for bar in range(Nel-1):
plt.plot([disp[T][bar],disp[T][bar+1]],[coords[bar],coords[bar+1]],'o-',c=color[i])
plt.plot([disp[T][Nel-1],disp[T][Nel]],[coords[Nel-1],coords[Nel]],'o-',c=color[i],label=f'T={T}s')
plt.legend(loc='lower right')
plt.savefig(absolute_path+'{0}.png'.format(str(num).zfill(3)))
plt.close(fig)
def save_gif_T(absolute_path, coords, all_displacements):
"""
all_displacements is a dictionnary where the keys are the sinus periods
"""
periods = [i for i in all_displacements.keys()]
disp_min = np.amin(list(all_displacements.values()))
disp_max = np.amax(list(all_displacements.values()))
nbtot = all_displacements[periods[0]].shape[0]
with imageio.get_writer('displacement.gif', mode='I') as writer:
for k in range(nbtot):
disp = {T:all_displacements[T][k] for T in periods}
if k%10 == 0: print('Plotting displacement {nb} out of {tot}'.format(nb=k,tot=nbtot))
plot_disp_T(absolute_path,coords,disp,disp_min,disp_max,k)
filename = absolute_path+'{0}.png'.format(str(k).zfill(3))
image = imageio.imread(filename)
writer.append_data(image)
# TOOLS
def compute_M_K(Nnodes, E, A, L, rho=1):
"""
Compute mass and rigidity matrices
"""
M = np.identity(Nnodes)
M[0,0] = 1/2
M[-1,-1] = 1/2
M = M*(L*A*rho)
Ktot = 2*np.identity(Nnodes) - np.diag(np.ones(Nnodes-1),1) - np.diag(np.ones(Nnodes-1),-1)
Ktot[0,0] = 1
Ktot[-1,-1] = 1
Ktot = Ktot*(E*A/L)
return M, Ktot
def compute_B_weights_old(Nnodes, coords):
"""
Compute strain-displacement matrix (B) and weights of each element
"""
Nel = Nnodes - 1
B = np.zeros((Nel, Nnodes)) # link element deformation to node displacement
weights = np.zeros((Nel,1)) # weights of each element
for i in range(Nel):
dX = np.abs(coords[i+1] - coords[i])
B[i,i] = -1/dX
B[i,i+1] = 1/dX
weights[i] = dX
return B, weights.ravel()
def compute_B_weights(Nnodes, A, L):
"""
Compute strain-displacement matrix (B) and weights of each element
"""
Nel = Nnodes - 1
B = np.zeros((Nel, Nnodes)) # link element deformation to node displacement
weights = np.zeros((Nel,1)) # weights of each element
for i in range(Nel):
B[i,i] = -1/L
B[i,i+1] = 1/L
weights[i] = A*L
return B, weights.ravel()
def compute_Keq(Ce, B, weights):
"""
Compute equivalent stiffness matrix (see DDM algorithm)
"""
return Ce*np.einsum('e,ei,ej->ij',weights,B,B)
def compute_BigMat(Keq, M, beta, dt):
"""
Compute big matrix used in DDM algorithm
"""
Mmod = M/(beta*dt**2)
return np.hstack((np.vstack((Keq,Mmod)),np.vstack((-Mmod,Keq))))
def compute_distance(epsi, sigma, weight, dataset, F_cost):
"""
Computes distances between (epsi,sigma) and the whole dataset (can be intensive)
"""
depsi = epsi-dataset[:,0]
dsigma = sigma-dataset[:,1]
return F_cost(depsi, dsigma, weight)
def find_closest(epsis, sigmas, weights, dataset, F_cost):
"""
Returns indices of closest dataset points to each point of (epsis,sigmas)
"""
if isinstance(weights, float):
e,s,w = epsis, sigmas, weights
distances = compute_distance(e,s,w,dataset,F_cost)
closest_indices = np.argmin(distances)
else:
Nel = weights.shape[0]
closest_indices = []
for i in range(Nel):
e,s,w = epsis[i], sigmas[i], weights[i]
distances = compute_distance(e,s,w,dataset,F_cost)
closest_index = np.argmin(distances)
closest_indices.append(closest_index)
return closest_indices

Event Timeline

c4science · Help