# MFS based MD Simulation

### A Sound scattering simulator based on the method of fundamental solutions


#### Contributors:
*Nicholas St. Clair, Physics (Kleckner Lab), University of California, Merced*

*Arnold Kim, Applied Mathematics, University of California, Merced*

*Dustin Kleckner, Physics, University of California, Merced*

#### Description:
The following notebook makes use of the method of fundamental solutions to compute the sound scattered from multiple sound hard scattering bodies. It determines the forces resultant from said scattering events, some of which are non-linear interactions brought about by interference of acoustic fields. The force data from the MFS scattering code is then fed into a Runge-Kutta based step integrator to produce a molecular dynamics type simulation of particle dynamics.

# 1. Relevant Libraries and Functions


```python
#Import relevant python libraries, as well as MFS master class
import numpy as np
import scipy.special as sp
import time
import matplotlib.pyplot as plt
from mfs import mfs as MFS
import csv
from scipy.integrate import solve_ivp
```


```python
'''A function to assign random initial positions to a group of particles,
the x & y coordinates are both drawn from a Gaussian distribution, whose spread relative to 
the particle size can be specified. A few if/else statements ensure that none of the particles
are initially overlapping.'''

def init_config(a,Np,σ=3):
    #σ: spread in distribution size
    #a: particle radius
    #Np: Number of particles for which to definie initial locations
    
    #X,Y are arrays that contain all of the coordinates of the particles in the config
    #Initialize the coordinates of the first particle of config in X,Y
    X = np.random.normal(loc = 0.0, scale = σ * a, size = 1 )
    Y = np.random.normal(loc = 0.0, scale = σ * a, size = 1 )
    #X_t,Y_t are temporary coordinates of the next particle to be added to the config
    X_t = np.random.normal(loc = 0.0, scale = σ * a, size = 1 )
    Y_t = np.random.normal(loc = 0.0, scale = σ * a, size = 1 )
    #ch is a checker used to indicate overlapping particles
    ch = 0

    #hold this loop open until coordinate arrays contain positions for Np particles
    while len(X) < Np:
        for i in range(len(X)):
            #Check if the next particle to be added is overlapping with any of the current particles
            if np.sqrt((X[i] - X_t)**2 + (Y[i] - Y_t)**2) < (2.0 * a):
                #if so add 1 to the checker variable
                ch += 1
            else:
                #otherwise do nothing
                ch += 0
        #if the checker variable is zero, it means the temporary particle position does not
        #overlap with any of the current particles 
        if ch == 0:
            #so we add the temporary position to the existing array of particle locations
            X = np.append(X,X_t)
            Y = np.append(Y,Y_t) 
        #Otherwise if the checker variable is nonzero, we have an overlap
        elif ch > 0:
            #so we define a new set of temporary coordinates to be tested
            X_t = np.random.normal(loc = 0.0, scale = σ * a, size = 1 )
            Y_t = np.random.normal(loc = 0.0, scale = σ * a, size = 1 )
            #and reset the checker variable
            ch = 0
    #Once we have a satisfying array of particle coordinates (size of array = number of desired particles,
    #and none of the particles overlap, return that array)
    return(X,Y)
```

## 2. Base differential equation
The following cell defines a function which takes a 1 X (6*Np) array (called M below) containing current position and velocity info and returns an array of the same shape with info updated according to the current velocity and force data (Np: Particle Number). The function also takes particle radius (a), particle mass (m), and dynamic viscosity of host fluid ($\mu$), in order to determine the Stokes Drag force on a moving particle.

*Unfortunately the integrator only accepts 1D arrays...so the differential equation to be solved must take the following form:*

\begin{equation}
     \begin{bmatrix}
         x_{i,t+\Delta t}\\
         y_{i,t+\Delta t}\\
         z_{i,t+\Delta t}\\
         v_{i,t+\Delta t}^{x}\\
         v_{i,t+\Delta t}^{y}\\
         v_{i,t+\Delta t}^{z}\\
         \vdots
     \end{bmatrix} \  = \
     \begin{bmatrix}
         v_{i,t}^{x}\\
         v_{i,t}^{y}\\
         v_{i,t}^{z}\\
         a_{i,t}^{x}\\
         a_{i,t}^{y}\\
         a_{i,t}^{z}\\
         \vdots
     \end{bmatrix}
     \Delta t, \ \ \ i=1,2,...Np
\end{equation}

The MFS class is responsible for determing the acoustic forces resultant on all particles, as well as any contact forces between pairs of particles that come too close. This in invoked by calling the force() and contact() methods from the MFS class. An instance of the scatter class is passed to the solver, this instance only need be called once before starting the solver.

The method of determining an acoustic radiation force on one or multiple particles is detailed in the Jupyter notebooks titled: "Acoustic Radiation Force on a single sound hard sphere", and "Acoustic forces between two and three spheres". The user is referred to those notebooks for a comprehensive understanding of the Method of Fundamental Solutions for determining the acoustically induced forces on sound-hard, spherical, scattering bodies.

The contact forces computed between pairs of particles are determined according to a power law scheme, whereby the force between two particles goes as the separation distance between those particles raised to some power:

\begin{equation}
        \vec{F}_{contact} = C*F_{0} \left[ \frac{\alpha - (\delta /d)}{\alpha - 1} \right]^{\beta} \hat{\delta}
\end{equation}

In the above $\alpha$ determines the "cutoff distance", i.e. the location within which to begin calculating the contact force. $\delta$ is the interparticle separation distance, d is the particle diameter, $\beta$ controls the sharpness of the force as it goes from zero to it's max value, and $C*F_{0}$ is the max value of the force, with $C$ ultimately controlling that value. $F_{0}$ is the reference acoustic force experienced by a perfect sound absorbing sphere within an external pressure field. The reference acoustic force is simply the acoustic energy density contained in the external field, multiplied by the particle volume, and is defined as:

\begin{equation}
    F_{0} = \frac{P_{0}^{2}}{4\rho_{0}c_{0}^2}\frac{4}{3}\pi a^{3}
\end{equation}

where $P_{0}^{2}$ is the pressure amplitude of the external field, $\rho_{0}$ is the mass density of the host medium, and $c_{0}$ is the sound speed in the host medium.


```python
def Binding(t, M, Np, a, m, μ, scat):
    #store current position and velocity data to separate arrays
    X = M.reshape(Np,6)[:,:3]
    V = M.reshape(Np,6)[:,3:6]
    #feed the class the current position data and use to compute forces
    scat.solve(X)
    Ft = scat.force() + scat.contact()
    #Determine current acceleration data (including viscous contributions)
    A = (Ft - 6*np.pi*μ*a*V)/m 
    #Add acceleration due to gravity in the z-direction
    A[:,2] -= 9.8
    #zip together current velocity and acceleration data to update particle kinematics
    Out = np.array(list(zip(V,A))).flatten()
    return Out;
```

## 3. Runge-Kutta based solver

The following cell defines initial conditions, as well as any global parameters necessary to set up the simulation.
It then feeds the initial conditions to the solver to be run for a predetermined amount of time.


```python
#time range of simulation
t_0 = 0
t_f = 10

#parameters for MFS class, as well as damping parameters used in differential equation
#MFS number (number of source points)
N = 92
#quadrature number (actual number of integration pointns is 2*Nq^2)
Nq = 8
# acoustic wavenumber (1/m)
k = 733
#wave speed (m/s)
c = 343
#angular frequency (rad/s)
ω = k*c
#size parameter
ka = 1.466
#particle radius (m)
a = ka / k
#mass density of air (kg/m^3)
ρ = 1.225
#incident field (velocity potential) amplitude (m^2 / s)
phi_a = 4e2 / (ρ * ω)
#particle mass density (kg/m^3)
ρp = 10
#particle mass (kg)
m = ((4 * np.pi / 3) * a ** 3) * ρp
#dynamic viscosity of air (kg/(m*s))
μ = 1.81e-5
#z level of group (m), generally set to be a stable levitation plane
zl = -0.0012542461099942367  
#number of particles
Np = 5
#variance in initial particle configuration spacing (in terms of 'a', the particle radius)
σ = 3

#Loop through number of desired trials
num_trials = 100
for n in range(num_trials):
    
    #set particle locations in initial config
    pos = np.empty((Np,3))
    pos[:,0], pos[:,1] = init_config(a,Np,σ)
    pos[:,2] = np.full(Np,zl)
    
    #set particle velocities in initial config (generally nil)
    vel = np.zeros((Np,3))
    ICs = np.empty((Np,6))
    
    #load particle locations and velocities into an initial condition array, fed to the integrator
    ICs[:, :3] = pos
    ICs[:, 3:] = vel

    #call MFS scatter class with relevant parameters
    
    scat = MFS.Scatter(N=N, k=k, a=a, phi_a = phi_a, rho = ρ, source_depth=0.5)
    scat.incoming_planewaves([1, 1j], [(0, 0, 1), (0, 0, -1)])
    #run the simulation
    nsol = solve_ivp(lambda t,M: Binding(t, M, Np, a, m, μ, scat), [t_0,t_f], ICs.flatten(), max_step = 0.05 )
    
    #uncomment the lines below in order to save the simulation data from each of the trials.
    
    #specifiy path of folder within which to save simulation data for each trial
    path = 'MD_Data_raw/'
    
    #specify a unique file name for each trial run, and include size parameter in the file name
    filename = f'MD_Data_ka={ka}_({n}).csv'
    
    #write trial data to the unique file using the python csv class
    with open(path+filename, mode = 'w') as csvfile:
        writer = csv.writer(csvfile, delimiter = ',')
        #write the time data as first row
        writer.writerow(nsol.t)
        
        #loop over number of time steps in simulation
        for i in range(len(nsol.t)):
            #write the kinematic data from each time step
            writer.writerow(nsol.y[:,i])
    
    #A separate notebook is used to analyze the data, said notebook is also included in the github repo.
```


```python

```

e681e2f6274979088da20f7e6fdefaed4e6fc7f8
ipynb
Jupyter Notebook
Jupyter Notebooks/MFS -- Molecular Dynamics Simulation -- Data creation.ipynb
klecknerlab/mfs
42fe3a5b1029c740a8139db92c710e69761cf224
Apache-2.0
Jupyter Notebooks/MFS -- Molecular Dynamics Simulation -- Data creation.ipynb
klecknerlab/mfs
42fe3a5b1029c740a8139db92c710e69761cf224
Apache-2.0
Jupyter Notebooks/MFS -- Molecular Dynamics Simulation -- Data creation.ipynb
klecknerlab/mfs
42fe3a5b1029c740a8139db92c710e69761cf224
Apache-2.0
Qwen/Qwen-72B
1. YES
2. YES