Merge pull request #28 from PCoulomPoMtl/main

Integration of proper bubble interactions structure

.gitignore

@@ -78,4 +78,18 @@ dkms.conf
 *.cmake
 CMakeCache.txt
 CMakeFiles/
-Makefile
+Makefile
+
+# Results
+*.txt
+!CMakeLists.txt
+!LICENSE.txt
+*.pdf
+!APECSS_Documentation.pdf
+!JOSS-Paper.pdf
+!apecssbubble.pdf
+!emissions_results.pdf
+!LagrangianWaveTracking_withShock.pdf
+!LagrangianWaveTracking.pdf
+!ultrasound_lipdcoated_simple.pdf
+*.png

README.md

@@ -20,6 +20,7 @@ APECSS is a software toolbox to compute pressure-driven bubble dynamics and the
 Key features of APECSS are:
 - Bubble dynamics using widely-used models (Rayleigh-Plesset, Keller-Miksis, Gilmore), solved using an in-built 5th-order Runge-Kutta scheme with adaptive time stepping.
 - Acoustic emissions of the bubble under different assumptions (incompressible, quasi-acoustic, fully compressible).
+- Interbubble interactions through their acoustic emissions
 - Prediction of the formation and attenuation of shock fronts emitted by the bubble.
 - Viscoelastic media (Kelvin-Voigt, Zener, Oldroyd-B).
 - Lipid monolayer coating of the bubble as used for ultrasound contrast agents.

examples/binaryinteraction/run.apecss

@@ -6,7 +6,7 @@
 
 BUBBLE
 RPModel KM
-Emissions IC 250.0e-6
+Emissions IC 300.0e-6
 PressureAmbient 1.0e5
 END
 

examples/binaryinteraction/src/binaryinteraction_apecss.c

@@ -21,6 +21,7 @@
 
 // Declaration of additional case-dependent functions
 APECSS_FLOAT interaction_bubble_pressure_infinity(APECSS_FLOAT t, struct APECSS_Bubble *Bubble);
+APECSS_FLOAT interaction_bubble_pressurederivative_infinity(APECSS_FLOAT t, struct APECSS_Bubble *Bubble);
 
 // Declaration of the structure holding the interaction variables of each bubble
 struct Interaction
@@ -95,18 +96,6 @@ int main(int argc, char **args)
   for (register int i = 0; i < nBubbles; i++) apecss_bubble_setdefaultoptions(Bubbles[i]);
   for (register int i = 0; i < nBubbles; i++) apecss_bubble_readoptions(Bubbles[i], OptionsDir);
 
-  // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-  // Allocate and set bubble-associated variables for the interaction
-  for (register int i = 0; i < nBubbles; i++)
-  {
-    struct Interaction *interaction_data = (struct Interaction *) malloc(sizeof(struct Interaction));
-    interaction_data->dp_neighbor = 0.0;  // Pressure excerted by the neighbor bubbles
-
-    // Hook interaction-structure to the void data pointer
-    Bubbles[i]->user_data = interaction_data;
-  }
-  // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-
   /* Allocate the structures for the fluid properties and ODE solver parameters */
   struct APECSS_Gas *Gas = (struct APECSS_Gas *) malloc(sizeof(struct APECSS_Gas));
   struct APECSS_Liquid *Liquid = (struct APECSS_Liquid *) malloc(sizeof(struct APECSS_Liquid));
@@ -164,6 +153,21 @@ int main(int argc, char **args)
   // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
   // Use the revised pressure at infinity, including neighbor contributions
   for (register int i = 0; i < nBubbles; i++) Bubbles[i]->get_pressure_infinity = interaction_bubble_pressure_infinity;
+  for (register int i = 0; i < nBubbles; i++) Bubbles[i]->get_pressurederivative_infinity = interaction_bubble_pressurederivative_infinity;
+
+  // Initialize interaction structure
+  for (register int i = 0; i < nBubbles; i++) Bubbles[i]->Interaction = (struct APECSS_Interaction *) malloc(sizeof(struct APECSS_Interaction));
+
+  // Update interaction structure
+  for (register int i = 0; i < nBubbles; i++)
+  {
+    Bubbles[i]->Interaction->nBubbles = nBubbles;
+    Bubbles[i]->Interaction->dp_neighbor = 0.0;
+    Bubbles[i]->Interaction->last_t_1 = 0.0;
+    Bubbles[i]->Interaction->last_t_2 = 0.0;
+    Bubbles[i]->Interaction->last_p_1 = 0.0;
+    Bubbles[i]->Interaction->last_p_2 = 0.0;
+  }
 
   // Define the size of each bubble
   for (register int i = 0; i < nBubbles; i++)
@@ -175,6 +179,23 @@ int main(int argc, char **args)
 
     Bubbles[i]->r_hc = Bubbles[i]->R0 / 8.54;
   }
+
+  // Define center location for each bubble
+  for (register int i = 0; i < nBubbles; i++)
+  {
+    if (0 == i)
+    {
+      Bubbles[i]->Interaction->location[0] = 0.0;
+      Bubbles[i]->Interaction->location[1] = 0.0;
+      Bubbles[i]->Interaction->location[2] = 0.0;
+    }
+    else if (1 == i)
+    {
+      Bubbles[i]->Interaction->location[0] = bubble_bubble_dist;
+      Bubbles[i]->Interaction->location[1] = 0.0;
+      Bubbles[i]->Interaction->location[2] = 0.0;
+    }
+  }
   // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 
   clock_t starttimebubble = clock();
@@ -203,17 +224,7 @@ int main(int argc, char **args)
 
     // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
     // Update the contribution of the neighbor bubble
-    struct Interaction *data_0 = Bubbles[0]->user_data;
-    data_0->dp_neighbor =
-        Liquid->rhoref *
-        (APECSS_POW2(Bubbles[1]->R) * Bubbles[1]->ode[0](Bubbles[1]->ODEsSol, Bubbles[1]->t, Bubbles[1]) + 2.0 * Bubbles[1]->R * APECSS_POW2(Bubbles[1]->U)) /
-        bubble_bubble_dist;
-
-    struct Interaction *data_1 = Bubbles[1]->user_data;
-    data_1->dp_neighbor =
-        Liquid->rhoref *
-        (APECSS_POW2(Bubbles[0]->R) * Bubbles[0]->ode[0](Bubbles[0]->ODEsSol, Bubbles[0]->t, Bubbles[0]) + 2.0 * Bubbles[0]->R * APECSS_POW2(Bubbles[0]->U)) /
-        bubble_bubble_dist;
+    apecss_interactions_instantaneous(Bubbles);
     // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
   }
 
@@ -246,6 +257,12 @@ int main(int argc, char **args)
 
 APECSS_FLOAT interaction_bubble_pressure_infinity(APECSS_FLOAT t, struct APECSS_Bubble *Bubble)
 {
-  struct Interaction *temp_struct = Bubble->user_data;
-  return (Bubble->p0 - Bubble->Excitation->dp * APECSS_SIN(2.0 * APECSS_PI * Bubble->Excitation->f * t) + temp_struct->dp_neighbor);
+  return (Bubble->p0 - Bubble->Excitation->dp * APECSS_SIN(2.0 * APECSS_PI * Bubble->Excitation->f * t) + Bubble->Interaction->dp_neighbor);
+}
+
+APECSS_FLOAT interaction_bubble_pressurederivative_infinity(APECSS_FLOAT t, struct APECSS_Bubble *Bubble)
+{
+  // Approximate numerical computation of p_infinity derivative
+  APECSS_FLOAT derivative = -Bubble->Excitation->dp * 2.0 * APECSS_PI * Bubble->Excitation->f * APECSS_COS(2.0 * APECSS_PI * Bubble->Excitation->f * t);
+  return (derivative);
 }

examples/cavitationonset/IC/build/CMakeLists.txt

@@ -0,0 +1,41 @@
+cmake_minimum_required(VERSION 3.12)
+
+project (cavitationonset_apecss)
+set(CMAKE_C_COMPILER /usr/bin/gcc)
+
+include_directories($ENV{APECSS_DIR}/include)
+FILE (GLOB_RECURSE myfiles ABSOLUTE ../src/*.c)
+
+set (mylibs m apecss)
+link_directories($ENV{USRLIB_DIR} $ENV{APECSS_DIR}/lib)
+
+foreach(arg ${myincludes})
+  IF (arg MATCHES "-I")
+    STRING(REGEX REPLACE "-I" "" myinc ${arg})
+    message("Additional include: ${myinc}")
+    include_directories(${myinc})
+  ENDIF(arg MATCHES "-I")
+endforeach(arg ${myincludes})
+
+foreach(arg ${mylibs})
+  STRING(REGEX REPLACE "lib" "" myl1 ${arg})
+	STRING(REGEX REPLACE ".a$" "" myl2 ${myl1})
+	set(mylibs ${myl2} ${mylibs} ${myl2})
+endforeach(arg ${mylibs})
+
+if(NOT CMAKE_BUILD_TYPE)
+   message(STATUS "Optimization: No optimization specified.")
+   set(CMAKE_BUILD_TYPE Release)
+endif()
+
+if(CMAKE_BUILD_TYPE STREQUAL "Debug")
+   message(STATUS "Optimization: Debug")
+   set(CMAKE_C_FLAGS_DEBUG "-Wall -g")
+elseif(CMAKE_BUILD_TYPE STREQUAL "Release")
+   message(STATUS "Optimization: Release")
+   add_definitions(-DNDEBUG)
+   set(CMAKE_C_FLAGS_RELEASE "-Wall -Werror -O3")
+endif()
+
+add_executable(cavitationonset_apecss ${myfiles})
+target_link_libraries(cavitationonset_apecss ${mylibs})

examples/cavitationonset/IC/build/compile.sh

@@ -0,0 +1,5 @@
+#!/bin/bash
+rm cavitationonset_apecss
+cmake CMakeLists.txt -DCMAKE_BUILD_TYPE=Release
+make
+rm -r CMakeCache.txt CMakeFiles Makefile cmake_install.cmake

examples/cavitationonset/IC/gather_results.py

@@ -0,0 +1,24 @@
+import os
+
+# File designed to recover data for plotting results in cavitation onset case
+
+file = open("Ida2009_results.txt", "r")
+lines = file.readlines()
+file.close()
+
+count = int(lines[0].split(" ")[0])
+png = float(lines[0].split(" ")[5])
+cl_size = float(lines[0].split(" ")[7])
+
+file_name = "{}_{:.4E}_{:.2f}.txt".format(count, png, cl_size)
+
+working_path = os.getcwd()
+results_path = os.path.join(working_path, "results")
+if not os.path.exists(results_path) :
+    os.makedirs(results_path)
+file_results = open(os.path.join(results_path, file_name), "w")
+
+for line in lines :
+    file_results.write(line)
+
+file_results.close()

examples/cavitationonset/IC/run.apecss

@@ -0,0 +1,36 @@
+#########################################################
+#                                                       #
+#  APECSS Options File                                  #
+#                                                       #
+#########################################################
+
+BUBBLE
+RPModel KM
+Emissions IC 900.0e-06
+PressureAmbient 0.1013e6
+END
+
+GAS
+EoS IG
+ReferenceDensity 1.2
+PolytropicExponent 1
+END
+
+LIQUID
+ReferencePressure 0.1013e6
+ReferenceDensity 1000.0
+ReferenceSoundSpeed 1500.0
+Viscosity 1.002e-3
+END
+
+INTERFACE
+SurfaceTensionCoeff 0.0728
+END
+
+RESULTS
+Bubble
+END
+
+ODESOLVER
+tolerance 1.0e-10
+END