FlowMol

(a) Total energy vs temperature from MD code compared with Rapaport (2004)
(b) Pressure vs temperature from MD code compared with Rapaport (2004)

Figure 12: (−) Rapaport Lennard-Jones (r_c = 2.5); (−−) Rapaport Weaks-Chandler-Anderson (WCA) (r_c = 2); (×) MD code Lennard-Jones (r_c = 2.5) (∘) MD code WCA; (■) ρ = 0.4; (■) ρ = 0.6; (■) ρ = 0.8; (■) ρ = 1.0

The velocity profile is compared with the analytical solution of the unsteady diffusion equation. Four layers of tethered molecules were used to model each wall, with the top wall given a sliding velocity of, U₀ = 1.0 at the start of the simulation, time t = 0. The temperature of both walls was controlled by applying the Nosé-Hoover (NH) thermostat to the wall atoms (Hoover, 1991). The two walls were thermostatted separately, and the equations of motion of the wall atoms were,

where n_x⁺ is a unit vector in the x− direction (only non-zero for the top wall), m_i ≡ m, and f_iext is the tethered atom force, obtained using the formula of Petravic & Harrowell (2006) (k₄ = 5 × 10³ and k₆ = 5 × 10⁶). The vector, r_i0 = r_i −r₀, is the displacement of the tethered atom, i, from its lattice site coordinate, r₀. The Nosé-Hoover thermostat dynamical variable is denoted by ξ, T₀ = 1.0 is the target temperature of the wall, and the effective time constant or damping coefficient, in Eq. (11d) was given the value, Q_ξ = NΔt. The simulation was carried out for a cubic domain of sidelength 27.40, of which the fluid region extent was 20.52 in the y−direction (see Figure 14a). Periodic boundaries were used in the streamwise (x) and spanwise (z) directions. The results presented in Figure 14b are the average of eight simulation trajectories starting with a different set of initial atom velocities. The lattice contained 16,384 molecules and was at a density of ρ = 0.8. The molecular simulation domain was sub-divided into 16 averaging slices of height 1.72 and the velocity was determined in each of them. The results show good agreement in both space and time, although some discrepancy is observed near the molecular walls. This is due to molecular layering and stick-slip behaviour near the walls in the molecular system which is not captured in the no-slip boundary assumed in the continuum solution.

The case of Couette flow is revisited repeatedly in this work, using the CFD code and in coupled simulations. This case is chosen because it is simple enough to afford an analytical solution and has been widely studied in the MD, CFD and coupling literature. The Couette flow case for a wide range of pressures has been studied fully, as part of a side project to the work presented herein (Heyes et al., 2012). The setup parameters used in this work are based on this published work Heyes et al. (2012), as are many of the insights into the behaviour of MD and coupled simulations.

3 Continuum code

3.1 Simple Finite Volume Solver

A two dimensional finite volume code which simulates the diffusive terms of the Navier Stokes equation has been developed to test the initial coupling. The unsteady diffusive equations are obtained by simplifying the 3 dimensional Navier Stokes Eq under the assumptions of developed laminar flow, a wide channel in the spanwise direction, negligible pressure gradient and no gravitational effects.

= .

(12)

uIJt+1 = u ijt + ,

(13)

In order to obtain the L₂ norm, the square of the discrepancy between analytical and numerical solutions at every point in the domain is obtained. The square root of the sum of these squared discrepencies is the L₂ norm. The L₂ norm error can be seen to decrease as the solution tends to the steady state and as the number of cells used (resolution) increases in figure 15b. The initial error is attributed to the use of an impulse started plate which the analytical solution would not be able to accurately re-create. The agreement between solutions appears to be good, verifying the CFD code outlined in this section and to be used in coupled Couette flow simulation.

4 Coupling of the DNS and MD algorithms

This section outlines the coupler library. The CPL library website is the most up to date source of information on the software. www.cpl-library.org This validation used version 1.0 of the CPL library code, and the text here is kept for historical records to highlight the validation cases performed during my thesis. Combining the capabilities of the DNS and MD codes requires interfacing and synchronising several elements of the two simulations. Each code models part of the physical problem, and is therefore responsible for a region of the computational domain. An overlap exists between those two domains; MD and continuum information must be “communicated” between the two codes in this region, in order to dynamically couple the two simulations.

Although this method has been demonstrated for a number of small cases, the most significant challenge is to carry out efficient and scalable (thousands of processors) computations of large continuum-molecular coupled systems. Due to its small size, the coupling case in Nie et al. (2004) has been reproduced using the minimal CFD code and the MD solveron a single processor. To move beyond this simple case, a framework must be developed to allow large scale coupling between the two parallel codes.

The proposed framework with separated coupler module and halo exchange is general enough to allow implementation of many of the existing coupling strategies (O’Connell & Thompson, 1995; Nie et al., 2004; Flekkoy et al., 2000; Hadjiconstantinou, 1998; Werder et al., 2005; Delgado-Buscalioni & Coveney, 2003). This generality ensures that the code development will continue to be useful as the field develops.

4.1 Coupler Outline

The coupler (CPL) consists of a series of function calls which are loosely based on the MPI framework in terms of both functionality and scope. They have been developed in FORTRAN 2008 with sufficient generality that they can be used as a language independent API through External functional interfaces. These routines are compiled into a library module which is linked to both the MD and CFD codes. Both codes are then run using the MPI multiple program multiple data paradigm (MPMD) with a call of the form,

mpiexec -n 128 ./md.exe : -n 8 ./cfd.exe.

The codes have entirely separate scopes and can only communicate through the coupler routines. The use of the MPMD framework enforces this separation and prevents conflicts between matching nomenclature in both the MD and CFD codes. It also pre-empts the use of dynamic processor allocation under MPI2, (MPI_spawn), providing a framework to support dynamic load balancing during a simulation.

The coupler consists of three classes of functions – setup, exchange and enquiry. Details of the structure and routines are included below:

• Setup The setup process is run once per simulation and entirely establishes the relationships between the two codes. The coupler setup routines and all coupler parameters are contained in a FORTRAN 2003 module with the protected attribute. This ensures only the three setup routines located inside this module can changes parameters internal to the coupler library. This prevents side-effects resulting from the interface with the CFD/MD codes as well as functions inside the coupler. The coupler setup functions have extensive error checking to ensure the inputs are of the correct form and the setup is completed correctly. A key part of the philosophy of the coupler is that both codes have copies of an identical set of read-only parameters including the complete knowledge of the processor topology of both regions. This ensures inter-communication between the codes is kept to a minimum.
  • CPL_create_comm – Setup MD/CFD INTRA_COMM and INTER_COMM between them. This should be called after both codes have initialised MPI. It defines a global communicator CPL_WORLD_COMM, local intra-communicators MD_COMM/CFD_COMM and an inter-communicator CPL_INTER_COMM between them to allow information exchange between the two domains Gropp et al. (1999b). The COUPLER.in input file is read to define relative processor topologies. The coupler input also includes parameters which are not specified from either of the coupled domains, such as the location and size of the overlap regions, the ratio of timestep between the two codes, as well as the facility to redefine variables in both codes in line with the coupled setup.
  • CPL_CFD/MD_init – Input CFD/MD parameters required to setup coupler parameters. Using a minimal set of parameters from both the CFD and MD code, together with the coupler input file, the entire coupling relationship between both domains is setup. The parameters are a minimal set of required and an extended range of optional arguments passed to the CPL_CFD/MD_init function. A hierarchy of communicators are established as part of the mapping process, including communicators containing only topologically overlapping processors CPL_OLAP_COMM and an optimised communicator based on the network topology of the overlapping processors CPL_GRAPH_COMM derived from the MPI_graph function. The details of the communicators are largely hidden from the user as described in the next section.
• Exchange The exchange routines are the core part of the coupler library. They include the following subroutines:
  • CPL_gather – CFD processors gather required data from overlapping MD processors.
  • CPL_scatter – CFD processors scatter required data to overlapping MD processors.
  • CPL_send – Send from current realm to overlapping processors on other realm.
  • CPL_recv – Receive on current realm to overlapping processors on other realm.

  All routines require two inputs; an array containing data for all the cells on the calling processor and the global cell numbers specifying the location and range (global extents) of the subset of these cells that the user wants to exchange. The use of global extents prevents the user from requiring any knowledge of the processor’s local numbering system, communicators or the mapping between domains. The global cell number conventions used in exchanges are identical on all processors, so the same input is used on the corresponding exchange routine in the other realm. Only processors with cells in the specified region send and receive the data between the domains. The data array can be of arbitrary size, allowing any number of properties to be exchanged (e.g. density (1), velocity (3) or stresses (9)). As an example, consider sending the MD velocity to be used as the CFD boundary condition for the setup discussed in the background chapter of Smith (2014) thesis. Each MD processor averages the velocity and passes its local array of 3 averaged velocity components per cell to CPL_send. The extents are specified to include all global CFD cells in xz plane at the CFD minimum cell location in y. Every CFD processor posts a corresponding CPL_recv with the same global extents. Every CFD processor on the bottom of the domain then receives an array containing the required boundary cells per processor.

• Enquiry A library of functions to retrieve parameters defined in the other realm or internal to the coupler. The use of intent out on all enquiry routines forces the user to make local copies of the returned variables. This is a further safeguard against corruption of data in the coupler and the associated side effects.
  • CPL_get – Returns any parameter considered protected public in the coupler library.
  • CPL_Cart_coords – Returns processor cartesian topology for any processor on either realm.
  • CPL_COMM_rank – Returns processor rank in any of the communicators
  • CPL_extents – Various extent routines to return cells on the current processor, in the overlap region or for arbitrary INPUT=limits subject to the limits of the current processor and overlap region.

    As the coupled topology is set-up at the beginning, all processors have a copy of all processor information in both regions. As a results, the enquiry routines require no communication.

As data from the MD and CFD codes must be manipulated before being passed to the CPL, it is recommended that coupler calls be done in socket routines which contain code written specifically for a given CFD or MD implementation. Calls to the coupler routines in either code can be removed by the use of pre-processing flags or dummy socket routines to allow each code to continue to function separately.

The following flowchart is available in Smith (2014) thesis shows the conceptual operation of a typical coupled setup.

4.2 Benchmarking and scalability of the coupled algorithm

As the coupling module is expected to load balance between both simulations and interface them efficiently, it is reasonable to expect that the coupled code will scale as well as the individual codes. The coupling region is local and can be decomposed spatially for parallel computations in the same manner as both of the codes it couples. As there is no greater requirement for communication in the coupling region, the observed scaling should be similar.

For the case of laminar Couette flow, the computational requirements of the continuum solver are almost negligible, and the continuum solver is run on a single processor with 64 cells in x and 7 in y. The speedup of the coupled code therefore depends almost entirely on the scaling of the MD solver and the coupler. If the coupler is performing efficiently, this combined speedup can be expected to be similar to the scaling of the MD solver alone. The scaling of the coupler is compared to the MD solver in Figure 17a, up to 1024 processes. The system size (N = 3,317,760) is the same as used in the previous all MD scaling tests, section 2.4. The coupled timestep ratio is 50.

Figure 17a demonstrates that the couplers performance only slightly deteriorates speedup. It is therefore not a severe bottleneck of performance of the coupled codes.

4.3 Verification of the Coupler

A part of the development of the coupler library, simple functions (e.g. sinusoidal profiles) are used to verify the passing is implemented correctly. A testing suite was written to ‘soak test’ the various coupler routines for a wide range of processor topologies in both the CFD and MD domains.

In this subsection, the coupler will be verified using three different types of coupling scheme from the literature. These include the state coupling of O’Connell & Thompson (1995) and Nie et al. (2004) both discussed in the state coupling section of Smith (2014) thesis as well as the flux coupling scheme of Flekkoy et al. (2000) discussed in the flux coupling section (thesis). A general schematic for the cases in this section is shown in Figure 18.

4.3.1 O’Connell & Thompson (1995)

The paper by O’Connell & Thompson (1995) presented an investigation of a sliding bottom wall in the molecular dynamics region with a fixed top wall in the continuum region. The bottom wall accelerates from stationary with a decaying exponential velocity of the form u_w(t) = U(1 − e^−t∕t₀) where U = 1 and t₀ = 160 in MD units. The results of O’Connell & Thompson (1995) are recreated in Figure 19 as part of the coupler verification. These results are compared to the time evolving analytical solution provided in the thesis of O’Connell (1995). The essential parameters of the setup are identical and the reader is referred to O’Connell & Thompson (1995) for full details. Only the main parameters and differences are summarised here. The molecular domain had dimensions 11.9 × 18.7 × 8.5 at ρ = 0.81 resulting in N = 1330 molecules. The cutoff length was r_c = 2.2 and the temperature T = 1.1. A bottom wall of thickness 2.34 was kept in place using the Petravic & Harrowell (2006) tethering potential introduced in section 2.5.4. The same density was used for the solid walls and the wall/fluid interaction potential ϵ = 1.0 was chosen instead of the 0.6 used by O’Connell & Thompson (1995). The domain was terminated by the O’Connell & Thompson (1995) boundary force applied at the top to prevent molecules escaping. Application of this boundary force results in a compression of the MD system causing the propagation of a shock wave. To prevent this effect, a number of molecules are removed to ensure the density remains at the intended ρ = 0.81 in the presence of the boundary force. A Nose-Hoover thermostat was applied instead of the Langevin thermostat used in O’Connell & Thompson (1995). A number of thermostatting strategies were tested including z-component only, profile unbiased thermostat (PUT) and wall only thermostatting. No discernible difference was found in the behaviour of the profile when compared to the analytical solution.

The setup of the coupled case was specified by the various parameters in Figure 18, including molecular domain top, y_MD⁺ = 18.7, CFD domain bottom y_CFD⁻ = 9.34, total coupled domain height H = 25.68 and cell size Δy = 2.34. The overlap region consists of four overlap cells. The arbitrary coupling strength coefficient is set to ξ = 0.01, as in O’Connell & Thompson (1995).

Molecules in the top MD cell experience the boundary force to prevent them escaping the domain. In addition, the top cell, together with the cell below, experiences a constrained force to force the average molecular velocity inside to match the continuum. The halo cell for the CFD was obtained by averaging the momentum in the pure MD region at the location covered by the bottom grey cell in Figure 18. Please see Smith (2014) thesis) for details of all forces and averaging methods.

The coupling was verified first using the 2D simple finite volume CFD solver outlined in section 3.1 with fixed top boundary condition and bottom boundary passed from the MD region. The CFD solver used a viscosity of μ = 2.14 (or equivalent Reynolds number of Re = 0.378) which was matched to the MD value obtained from Green-Kubo (Green, 1954; Kubo, 1957) calibration calculations.

Good agreement was observed to the analytical solution for an accelerating plate. The setup of this problem applies the moving boundary condition in the molecular region, which in turn drives the continuum due to the coupling between them.

This model implements no mass flux and the scaling parameters ξ is said to cause lagging behind the continuum solution in time evolving flows by Nie et al. (2004). The coupling method presented in the work of Nie et al. (2004) is claimed to ameliorate this problem. The implementation of this method using the CPL library is presented in the next section.

4.3.2 Nie et al. (2004)

This section implements the coupled simulation of Nie et al. (2004). The setup of the coupled simulation in this section aims to be similar to that of the previous section while matching the essential parameters of Nie et al. (2004). As before, the reader is referred to the original work (Nie et al., 2004) for full details of the setup used. The simulation is impulse started Couette flow, where the top wall in the CFD domain was moved at a velocity of u_w(t) = 1 and the bottom wall in the molecular region was stationary.

The MD domain size was 34.06 × 34.06 × 34.06 with a density of ρ = 0.81 giving N = 27,200 molecules. The molecular domain was simulated on 4 processes with topology 2 × 2 × 1. The bigger domain size provides better statistics and each cell was twice the size (in each direction) of those used by O’Connell & Thompson (1995). Nie et al. (2004) use the average of 10 MD ensembles, each of width 4.82, to provide boundary conditions for a single CFD instance. The implementation here instead uses a wider domain in z and a single MD simulation. A bottom wall of thickness 2.62 was kept in place using the Petravic & Harrowell (2006) tethering potential introduced in the NEMD section 2.5.4. The same density was used for the solid walls and fluid regions and for simplicity the wall/fluid interaction potential is set to ϵ = 1.0, instead of the 0.6 used by Nie et al. (2004). The cutoff length was r_c = 2.2 and the temperature T = 1.1. A Nosè-Hoover thermostat was applied instead of the Langevin. Again application of the thermostat in z only, in a profile unbiased manner or simply thermostatting the wall results in similar agreement to the Couette flow time-evolving analytical solution shown in Figure 20.

The setup of the coupled case was specified by molecular domain top, y_MD⁺ = 31.3, CFD domain bottom y_CFD⁻ = 15.6, total coupled domain height H = 52.1 and cell size Δy = 5.12. The overlap was of size 15.6, split into four overlap cells. Molecules in the top MD cell are subjected to the boundary forces to prevent them from escaping. In addition, the top cell, together with the cell below, experiences a constrained forc to ensure average molecular velocties are matched to the continuum. The halo for the CFD was obtained by averaging the pure MD region below the overlap region. The ratio of timesteps between the MD and CFD simulations was 50. Please see Smith (2014) thesis) for details of all forces and averaging methods.

The coupling was verified using the 2D simple finite volume CFD solver and full 3D DNS code ranslow (run on a single processor). Both CFD solvers used a viscosity of μ = 2.14 (or equivalent Reynolds number).

The velocity against time shows good agreement to the analytical solution in Figure 21a.

The stress in both the CFD and MD simulation can also be evaluated for the coupled case and the agreement compared to the analytical solution for stress (see appendix thesis). In the MD region, the stress was obtained from the volume average (VA) localisation of the virial stress (Cormier et al., 2001) discussed in the VA section of Smith (2014) thesis. In the CFD region, the stress can be calculated using the hydrodynamic pressure, the gradient of velocity and the shear viscosity. The stress in both simulations and across the coupled region can be seen to match the analytical solution as shown in Figure 21b. The constraint on velocity can be seen to ensure the stress is also matched between both domain. In the next section, the stress itself is constrained, using the flux coupling schemes of Flekkoy et al. (2000).

4.3.3 Flekkoy et al. (2000)

The final verification case is based on the flux coupling work of Flekkoy et al. (2000). They model Couette flow using an MD domain between two CFD regions. The results they present are only for steady state Couette flow. This section extends these results to verify the time evolution of a coupled Couette flow with flux coupling. A much larger molecular domain was used to test the coupler for large system sizes. This also addresses the higher averaging requirement of flux coupling outlined in statistical error section of Smith (2014) thesis. The MD domain size was 164.2 × 51.3 × 41.0 at density ρ = 0.8, resulting in N = 276,480 molecules. The molecular domain was simulated on 32 processes with topology 8 × 2 × 2. The bottom wall was tethered with a thickness of 5.13 and a Nosé-Hoover thermostat was applied to only the wall molecules with a temperature setpoint of T = 1.0. The WCA potential was used (r_c = 2) for computational efficiency.

The paper by Flekkoy et al. (2000) is unclear of exact domain sizes but states that cells of size 7.6 are employed with 10 cells shown in the MD domain on their Fig. 2 (i.e. an MD domain of 76.0). This MD region is sandwiched between two CFD domains in Flekkoy et al. (2000) with a coupling region on each side of width 22.8. Only one sided coupling was used here, again based on the setup shown in Figure 18. A domain of similar height without the bottom coupling region was used – y_MD⁺ = 46.17. The cells were chosen to be slightly smaller than Flekkoy et al. (2000) and consistent with Nie et al. (2004) at Δy = 5.13. The continuum domain bottom was located at y_CFD⁻ = 25.65 and the total domain was twice the size of Nie et al. (2004) with a height of H = 107.8

The constraint force of Flekkoy et al. (2000) with weighting function given in Smith (2014) thesis) applied both the appropriate direct pressure and shear pressure to the MD domain. As the pressure applied is not designed to prevent molecules from escaping, only to match the continuum pressure, specular walls are also employed. These reflect molecules back with identical x and z velocities and exactly opposite y components of velocity. Specular walls are chosen as the continuum solver was incompressible and does not expect a mass flow into the bottom of the continuum domain, therefore restricting mass flux is consistent. Moving these walls with a mean velocity, as in the work of (Werder et al., 2005), was unnecessary as the mean flow rate in the wall normal direction is, on average, zero for this case.

The viscosities are also matched in this simulation with continuum viscosity set to μ = 1.6 (again from Green-Kubo autocorrelation functions). The coupling was between the MD and the full 3D DNS code ranslow with Reynolds number Re = 0.5. The parallel nature of the DNS code was employed with the 32 × 16 × 8 CFD cells split between 4 × 2 × 1 processes.

The coupled solution for velocity is compared to the analtical solution in Figure 22. Good agreement is observed at equivalent times to those shown for Nie et al. (2004) in the previous section. As the domain is over 3 times the size of the previous Nie et al. (2004) case, the time scales are proportional to the square of the domain height and are scaled accordingly. Note the greater number of symbolds due to the much greater numbers of cells in both domains.

As the stress in both regions was matched by the coupler, this is seen to match in Figure 23b. The velocity can also be seen to agrees in Figure 23a as a result of this stress matching scheme. The near wall agreement in Figure 23a appears to be superior to the Nie et al. (2004) coupling shown in Figure 21a. This is a consequence of the larger system size which reduces the effect of near wall molecular behaviour and improves velocity averages.

The halo for the CFD (bottom grey cells in Figure 18) was obtained by averaging the momentum in the pure MD domain below the overlap region in line. This is a departure from the work of Flekkoy et al. (2000) who average the stresses from the MD region. The manner in which to obtain the stress tensor from the MD is not clear (as discussed in section flux pressure in Smith E (2014) thesis. This mixed coupling is based on the work of Ren (2007), where flux-state coupling is found to be stable, while flux-flux coupling was found to be unstable. There are also concerns about the sample requirements, which are said to be prohibitive for flux schemes in Hadjiconstantinou et al. (2003). However, a parameter study of error in appendix of (thesis) which demonstrate these averages requirements for flux shemes, while often large, are not prohibitive for all cases. Despite this, the mixed flux-state exchange was deemed to be a simpler and more robust way to test the coupling software.

Having verified the coupler gives identical results using both flux and state coupling, the impact of specifying the wrong viscosity is evaluated in the next section. The agreement of state based coupling schemes is shown to be contingent on a consistent viscosity being specified in both domains. The impact of this finding is then discussed.

4.4 Mismatched Viscosity in Couette Flow

The case of coupled Couette flow with a mismatch of viscosity between the continuum and molecular systems is considered in this subsection.

When coupling the continuum and molecular systems, the scaling between the two regions must be consistent. As discussed in class="cmbx-10x-x-109">(thesis) the continuum and molecular domains are scaled based on different parameters. The molecular region is based on length scale ℓ, mass m and energy ϵ, while the continuum is dependent on the Reynolds number Re = ρUL∕μ. The two regions are matched by ensuring the length scales ℓ vs. L and densities ρ_CFD vs. ρ_MD = ρ_MD(m,ℓ) are consistent. The velocity or stress matching is then ensured by the coupling mechanism. Viscosity is a tuneable parameter in the continuum but it is a measureable quantity in a molecular simulation, not an input parameter. Nie et al. (2004) use a viscosity of μ = 2.14 based on previous simulations. As shown in the previous subsection, using this pre-matching of viscosity, good agreement is observed for both the velocity and stress between both systems, Figure 21. For the case where the viscosities are not the same in both regions, the velocity profiles can be seen to match but the stresses do not, Figure 24.

Consider next the flux couping of Flekkoy et al. (2000). The use of the same viscosity in both domain (μ = 1.6) obtained by prior simulation, Figure 23 in the previous subsection showed good agreement for both velocity and stresses. The use of different viscosity in both regions results in disagreement in the velocity but good agreement in the stresses, Figure 25.

The results for a mismatched stress in figure 25 are consistent with the simulation of two phase fluids with different viscosities. Molecular fluid will typically experience different viscosity measurements in the vicinity of the walls or at high shear rates (shear thinning). By using state coupling, this behaviour would not be correctly transmitted to the CFD realm as a constant viscosity is implicitly assumed throughout the MD. In addition, it is possible to unknowingly simulate cases with viscosity difference and still obtain apparently correct velocity profiles. The conclusion is that the flux coupling scheme are more general as they require no assumption of the viscosity in the MD region, as they simply apply a force consistent with the current state of stress in the CFD.

The results from this study of mismatched viscosity suggest that a more stringent set of tests are required in order to establish the consistency between the two realms of a coupled simulation, a problem that motivates the next page with the development of a consistent control volume framework for both systems. In addition, flux coupling scheme are more general and should be preferred in the computational implementation of a general purpose coupling simulation tool, which was the aim of this section. As they are not currently derived from a constrained dynamics approach, some theoretical development is required and this motivates the work on control volumes (thesis).

This section has introduced a coupling library for simulation of large scale case on many processors. The intention was to couple simulations which approach scales of engineering interest, including the modelling of turbulent phenomenon. However, there are a number of problems to overcome before a generalised coupling scheme can be implemented.

5 Overview

This section includes the details of the software developments and verifications in the development of a coupled continuum-molecular solver. In order to couple molecular dynamics (MD) to continuum computational fluid mechanics (CFD), a new MD code was developed, a CFD code was adapted and a coupling library designed to manage the data exchange between them.

In section 2, the details of the development and testing of a new MD code were presented. The MD code was designed to use cell/neighbour lists, Hilbert curve re-ordering and other optimisations to obtain similar serial speeds to LAMMPS. The code is then parallelised using MPI to allow running on high performance computing platforms. The parallel performance is tested and shows 94% efficiency when scaling from 8 processors to 1024. The MD code was verified by checking the energy conservation, trajectory agreement (and divergence), radial distribution function, phase diagrams and simulation of Couette flow. The optimised MD code is therefore ready for large-scale simulation of molecular dynamics and coupling problems.

Section, 3, outlined the development and verification of a simple CFD code. This was used as part of the development of the coupling library and for simple CFD test cases.

The final part of this page, section 4, outlined the development of a robust and minimal coupling (CPL) library. This was designed as a language independent APIs using FORTRAN 2008, which implements the data exchange between the CFD and MD codes. Implementation of coupling should be as simple as initialising the CPL module by passing the key simulation parameters from both MD/CFD codes. The data can be exchanged by CPL_send/CPL_recv calls which ensure all data is passed correctly. The library was soak tested and verified for range of cases and topologies from the literature (O’Connell & Thompson, 1995; Nie et al., 2004; Flekkoy et al., 2000). In order to test the scaling of the coupler, performance of a CPL simulation was compared to an uncoupled case. The difference in speedup on 1024 processor was found to be 575 for the full MD and 449 for a coupled case, compared to the ideal 1024. The coupler is therefore seen to have minimal impact on the scaling of the two codes.

References