diff --git a/.gitignore b/.gitignore index de32c1e..1b55076 100644 --- a/.gitignore +++ b/.gitignore @@ -78,4 +78,18 @@ dkms.conf *.cmake CMakeCache.txt CMakeFiles/ -Makefile \ No newline at end of file +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 \ No newline at end of file diff --git a/README.md b/README.md index c24b551..7bb5a8b 100644 --- a/README.md +++ b/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. diff --git a/examples/binaryinteraction/run.apecss b/examples/binaryinteraction/run.apecss index 97fd792..779d9bc 100644 --- a/examples/binaryinteraction/run.apecss +++ b/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 diff --git a/examples/binaryinteraction/src/binaryinteraction_apecss.c b/examples/binaryinteraction/src/binaryinteraction_apecss.c index 6fc2e9e..c4c91a3 100644 --- a/examples/binaryinteraction/src/binaryinteraction_apecss.c +++ b/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); } \ No newline at end of file diff --git a/examples/cavitationonset/IC/build/CMakeLists.txt b/examples/cavitationonset/IC/build/CMakeLists.txt new file mode 100644 index 0000000..8e43b7f --- /dev/null +++ b/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}) \ No newline at end of file diff --git a/examples/cavitationonset/IC/build/compile.sh b/examples/cavitationonset/IC/build/compile.sh new file mode 100755 index 0000000..02b82c5 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/examples/cavitationonset/IC/gather_results.py b/examples/cavitationonset/IC/gather_results.py new file mode 100644 index 0000000..79929e0 --- /dev/null +++ b/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() \ No newline at end of file diff --git a/examples/cavitationonset/IC/run.apecss b/examples/cavitationonset/IC/run.apecss new file mode 100644 index 0000000..66cab93 --- /dev/null +++ b/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 diff --git a/examples/cavitationonset/IC/src/cavitationonset_apecss.c b/examples/cavitationonset/IC/src/cavitationonset_apecss.c new file mode 100644 index 0000000..1dd191a --- /dev/null +++ b/examples/cavitationonset/IC/src/cavitationonset_apecss.c @@ -0,0 +1,317 @@ +// This source file is part of APECSS, an open-source software toolbox +// for the computation of pressure-driven bubble dynamics and acoustic +// emissions in spherical symmetry. +// +// Copyright (C) 2022-2024 The APECSS Developers +// +// The APECSS Developers are listed in the README.md file available in +// the GitHub repository at https://github.com/polycfd/apecss. +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at https://mozilla.org/MPL/2.0/. + +// ------------------------------------------------------------------- +// APECSS standalone example of cavitation onset with acoustically-in- +// teracting microbubbles, based on Ida (2009), Physics of Fluids 21 +// (11), 113302, DOI : 10.1063/1.3265547 +// ------------------------------------------------------------------- + +#include +#include "apecss.h" + +// 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); + +int main(int argc, char **args) +{ + char OptionsDir[APECSS_STRINGLENGTH]; + + // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + // Interbubble time-step, defining the frequency with which the neighbor influence is updated + APECSS_FLOAT dt_interbubble = 1.0e-8; + + // Number of bubbles + const int nBubbles = 2; + + // Initialize the simulation parameters given by the execution command + double tEnd = 0.0; + double fa = 0.0; + double pa = 0.0; + double cluster_size = 0; + // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + apecss_infoscreen(); + + /* Read commandline options */ + sprintf(OptionsDir, "./run.apecss"); // This is the default + int j = 1; // First argument is the call to the executable + while (j < argc) + { + if (strcmp("-options", args[j]) == 0) + { + sprintf(OptionsDir, "%s", args[j + 1]); + j += 2; + } + else if (strcmp("-tend", args[j]) == 0) + { + sscanf(args[j + 1], "%le", &tEnd); + j += 2; + } + else if (strcmp("-freq", args[j]) == 0) + { + sscanf(args[j + 1], "%le", &fa); + j += 2; + } + else if (strcmp("-amp", args[j]) == 0) + { + sscanf(args[j + 1], "%le", &pa); + j += 2; + } + else if (strcmp("-clsize", args[j]) == 0) + { + sscanf(args[j + 1], "%le", &cluster_size); + j += 2; + } + else if (strcmp("-dt_inter", args[j]) == 0) + { + sscanf(args[j + 1], "%le", &dt_interbubble); + j += 2; + } + else if (strcmp("-h", args[j]) == 0) + { + apecss_helpscreen(); + } + else + { + char str[APECSS_STRINGLENGTH_SPRINTF]; + sprintf(str, "Unknown command line options: %s", args[j]); + apecss_erroronscreen(1, str); + ++j; + } + } + dt_interbubble = (APECSS_FLOAT) dt_interbubble; + + /* Allocate and initialize Bubble structure */ + struct APECSS_Bubble *Bubbles[nBubbles]; + for (register int i = 0; i < nBubbles; i++) Bubbles[i] = (struct APECSS_Bubble *) malloc(nBubbles * sizeof(struct APECSS_Bubble)); + for (register int i = 0; i < nBubbles; i++) apecss_bubble_initializestruct(Bubbles[i]); + + /* Set default options and read the options for the bubble */ + 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 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)); + struct APECSS_Interface *Interface = (struct APECSS_Interface *) malloc(sizeof(struct APECSS_Interface)); + struct APECSS_NumericsODE *NumericsODE = (struct APECSS_NumericsODE *) malloc(sizeof(struct APECSS_NumericsODE)); + + /* Set the default options for the fluid properties and solver parameters */ + apecss_gas_setdefaultoptions(Gas); + apecss_liquid_setdefaultoptions(Liquid); + apecss_interface_setdefaultoptions(Interface); + apecss_odesolver_setdefaultoptions(NumericsODE); + + /* Read the options file for the fluid properties and solver parameters */ + apecss_gas_readoptions(Gas, OptionsDir); + apecss_liquid_readoptions(Liquid, OptionsDir); + apecss_interface_readoptions(Interface, OptionsDir); + apecss_odesolver_readoptions(NumericsODE, OptionsDir); + + /* Associate the bubble with the relevant fluid properties and solver parameters */ + for (register int i = 0; i < nBubbles; i++) + { + Bubbles[i]->Gas = Gas; + Bubbles[i]->Liquid = Liquid; + Bubbles[i]->Interface = Interface; + Bubbles[i]->NumericsODE = NumericsODE; + } + + /* Allocate and set the excitation parameters */ + struct APECSS_Excitation *Excitation = (struct APECSS_Excitation *) malloc(sizeof(struct APECSS_Excitation)); + Excitation->type = APECSS_EXCITATION_SIN; + Excitation->f = (APECSS_FLOAT) fa; + Excitation->dp = (APECSS_FLOAT) pa; + for (register int i = 0; i < nBubbles; i++) Bubbles[i]->Excitation = Excitation; + + // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + // Create individual folders for the results of each bubble + for (register int i = 0; i < nBubbles; i++) + { + if (Bubbles[i]->Results != NULL) + { + sprintf(Bubbles[i]->Results->dir, "./Bubble_%i/", i); + struct stat st = {0}; + if (stat(Bubbles[i]->Results->dir, &st) == -1) mkdir(Bubbles[i]->Results->dir, 0700); + } + } + // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + /* Process all options */ + apecss_gas_processoptions(Gas); + apecss_liquid_processoptions(Liquid); + apecss_interface_processoptions(Interface); + apecss_odesolver_processoptions(NumericsODE); + for (register int i = 0; i < nBubbles; i++) apecss_bubble_processoptions(Bubbles[i]); + + // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + // 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 + Bubbles[0]->R0 = 2.0e-06; + Bubbles[1]->R0 = 20.0e-06; + + // Define center location for each bubble + Bubbles[0]->Interaction->location[0] = 0.0; + Bubbles[0]->Interaction->location[1] = 0.0; + Bubbles[0]->Interaction->location[2] = 0.0; + + Bubbles[1]->Interaction->location[0] = (APECSS_FLOAT) cluster_size * (Bubbles[0]->R0 + Bubbles[1]->R0); + Bubbles[1]->Interaction->location[1] = 0.0; + Bubbles[1]->Interaction->location[2] = 0.0; + // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + clock_t starttimebubble = clock(); + APECSS_FLOAT tSim = 0.0; // Simulation time of the coupled system + + for (register int i = 0; i < nBubbles; i++) + { + Bubbles[i]->tStart = tSim; + Bubbles[i]->tEnd = (APECSS_FLOAT) tEnd; + Bubbles[i]->dt = APECSS_MIN(1.0e-11, dt_interbubble); // Initial time-step + } + + /* Initialize the bubble based on the selected options */ + for (register int i = 0; i < nBubbles; i++) apecss_bubble_initialize(Bubbles[i]); + + /* Initialize */ + for (register int i = 0; i < nBubbles; i++) apecss_bubble_solver_initialize(Bubbles[i]); + + // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + // File to retrieve all valuable information for cavitation onset test case + FILE *file_ida2009; + file_ida2009 = fopen("Ida2009_results.txt", "w"); + fprintf(file_ida2009, "%d Bubbles p0(pa) %e png(Pa) %e D_multiplier(-) %e Interaction-type %d\n", nBubbles, Liquid->pref, pa, cluster_size, 1); + fprintf(file_ida2009, "Initial_radii(m)"); + for (register int i = 0; i < nBubbles; i++) fprintf(file_ida2009, " %e", Bubbles[i]->R0); + fprintf(file_ida2009, "\n"); + fprintf(file_ida2009, "#Time(s) R(m) Pt(Pa)\n"); + // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + /* Solve the bubble dynamics */ + while (tSim < (APECSS_FLOAT) tEnd) // Interaction loop, corresponding to the time-intervals at which interactions are considered + { + APECSS_FLOAT dtSim = APECSS_MIN(dt_interbubble, (APECSS_FLOAT) tEnd - tSim); + tSim += dtSim; + + for (register int i = 0; i < nBubbles; i++) apecss_bubble_solver_run(tSim, Bubbles[i]); + + // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + // Update the contribution of the neighbor bubbles + apecss_interactions_instantaneous(Bubbles); + // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + // Retrieve data + fprintf(file_ida2009, "%e", tSim); + for (register int i = 0; i < nBubbles; i++) + { + fprintf(file_ida2009, " %e", Bubbles[i]->R); + } + for (register int i = 0; i < nBubbles; i++) + { + fprintf(file_ida2009, " %e", Bubbles[i]->get_pressure_infinity(tSim, Bubbles[i])); + } + for (register int i = 0; i < nBubbles; i++) + { + fprintf(file_ida2009, " %e", Bubbles[i]->ode[0](Bubbles[i]->ODEsSol, Bubbles[i]->t, Bubbles[i])); + } + for (register int i = 0; i < nBubbles; i++) + { + fprintf(file_ida2009, " %e", Bubbles[i]->Interaction->dp_neighbor); + } + for (register int i = 0; i < nBubbles; i++) + { + fprintf(file_ida2009, " %e", Bubbles[i]->U); + } + fprintf(file_ida2009, "\n"); + // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + } + + fclose(file_ida2009); + + /* Finalize the simulation*/ + for (register int i = 0; i < nBubbles; i++) apecss_bubble_solver_finalize(Bubbles[i]); + + char str[APECSS_STRINGLENGTH_SPRINTF]; + for (register int i = 0; i < nBubbles; i++) + { + sprintf(str, "Bubble %i: Solver concluded %i time-steps and %i sub-iterations in %.3f s.", i, Bubbles[i]->dtNumber, Bubbles[i]->nSubIter, + (double) (clock() - starttimebubble) / CLOCKS_PER_SEC); + apecss_writeonscreen(str); + } + + for (register int i = 0; i < nBubbles; i++) apecss_results_rayleighplesset_write(Bubbles[i], APECSS_RESULTS_WRITE); + for (register int i = 0; i < nBubbles; i++) apecss_results_emissionsspace_write(Bubbles[i], APECSS_RESULTS_WRITE); + + /* Make sure all allocated memory is freed */ + for (register int i = 0; i < nBubbles; i++) apecss_bubble_freestruct(Bubbles[i]); + + for (register int i = 0; i < nBubbles; i++) free(Bubbles[i]); + free(Gas); + free(Liquid); + free(Interface); + free(NumericsODE); + free(Excitation); + + return (0); +} + +APECSS_FLOAT interaction_bubble_pressure_infinity(APECSS_FLOAT t, struct APECSS_Bubble *Bubble) +{ + APECSS_FLOAT T = 10.0e-06; + if (t < T) + { + return (Bubble->p0 + Bubble->Interaction->dp_neighbor); + } + else if (t > 2 * T) + { + return (Bubble->Excitation->dp + Bubble->Interaction->dp_neighbor); + } + else + { + APECSS_FLOAT W = 0.5 * (1 - APECSS_COS((t + T) * APECSS_PI / T)); + return (Bubble->p0 + W * (Bubble->Excitation->dp - Bubble->p0) + 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 T = 10.0e-06; + APECSS_FLOAT derivative = 0.0; + if ((t >= T) && (t <= 2 * T)) + { + APECSS_FLOAT inv_T = 1 / T; + derivative = 0.5 * APECSS_PI * inv_T * APECSS_SIN(APECSS_PI * (t + T) * inv_T) * (Bubble->Excitation->dp - Bubble->p0); + } + return (derivative); +} \ No newline at end of file diff --git a/examples/cavitationonset/NI/build/CMakeLists.txt b/examples/cavitationonset/NI/build/CMakeLists.txt new file mode 100644 index 0000000..8e43b7f --- /dev/null +++ b/examples/cavitationonset/NI/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}) \ No newline at end of file diff --git a/examples/cavitationonset/NI/build/compile.sh b/examples/cavitationonset/NI/build/compile.sh new file mode 100755 index 0000000..02b82c5 --- /dev/null +++ b/examples/cavitationonset/NI/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 \ No newline at end of file diff --git a/examples/cavitationonset/NI/gather_results.py b/examples/cavitationonset/NI/gather_results.py new file mode 100644 index 0000000..79929e0 --- /dev/null +++ b/examples/cavitationonset/NI/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() \ No newline at end of file diff --git a/examples/cavitationonset/NI/run.apecss b/examples/cavitationonset/NI/run.apecss new file mode 100644 index 0000000..66cab93 --- /dev/null +++ b/examples/cavitationonset/NI/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 diff --git a/examples/cavitationonset/NI/src/cavitationonset_apecss.c b/examples/cavitationonset/NI/src/cavitationonset_apecss.c new file mode 100644 index 0000000..c9971b4 --- /dev/null +++ b/examples/cavitationonset/NI/src/cavitationonset_apecss.c @@ -0,0 +1,300 @@ +// This source file is part of APECSS, an open-source software toolbox +// for the computation of pressure-driven bubble dynamics and acoustic +// emissions in spherical symmetry. +// +// Copyright (C) 2022-2024 The APECSS Developers +// +// The APECSS Developers are listed in the README.md file available in +// the GitHub repository at https://github.com/polycfd/apecss. +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at https://mozilla.org/MPL/2.0/. + +// ------------------------------------------------------------------- +// APECSS standalone example of cavitation onset with acoustically-in- +// teracting microbubbles, based on Ida (2009), Physics of Fluids 21 +// (11), 113302, DOI : 10.1063/1.3265547 +// ------------------------------------------------------------------- + +#include +#include "apecss.h" + +// 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); + +int main(int argc, char **args) +{ + char OptionsDir[APECSS_STRINGLENGTH]; + + // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + // Interbubble time-step, defining the frequency with which the neighbor influence is updated + APECSS_FLOAT dt_interbubble = 1.0e-8; + + // Number of bubbles + const int nBubbles = 2; + + // Initialize the simulation parameters given by the execution command + double tEnd = 0.0; + double fa = 0.0; + double pa = 0.0; + double cluster_size = 0; + // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + apecss_infoscreen(); + + /* Read commandline options */ + sprintf(OptionsDir, "./run.apecss"); // This is the default + int j = 1; // First argument is the call to the executable + while (j < argc) + { + if (strcmp("-options", args[j]) == 0) + { + sprintf(OptionsDir, "%s", args[j + 1]); + j += 2; + } + else if (strcmp("-tend", args[j]) == 0) + { + sscanf(args[j + 1], "%le", &tEnd); + j += 2; + } + else if (strcmp("-freq", args[j]) == 0) + { + sscanf(args[j + 1], "%le", &fa); + j += 2; + } + else if (strcmp("-amp", args[j]) == 0) + { + sscanf(args[j + 1], "%le", &pa); + j += 2; + } + else if (strcmp("-clsize", args[j]) == 0) + { + sscanf(args[j + 1], "%le", &cluster_size); + j += 2; + } + else if (strcmp("-dt_inter", args[j]) == 0) + { + sscanf(args[j + 1], "%le", &dt_interbubble); + j += 2; + } + else if (strcmp("-h", args[j]) == 0) + { + apecss_helpscreen(); + } + else + { + char str[APECSS_STRINGLENGTH_SPRINTF]; + sprintf(str, "Unknown command line options: %s", args[j]); + apecss_erroronscreen(1, str); + ++j; + } + } + dt_interbubble = (APECSS_FLOAT) dt_interbubble; + + /* Allocate and initialize Bubble structure */ + struct APECSS_Bubble *Bubbles[nBubbles]; + for (register int i = 0; i < nBubbles; i++) Bubbles[i] = (struct APECSS_Bubble *) malloc(nBubbles * sizeof(struct APECSS_Bubble)); + for (register int i = 0; i < nBubbles; i++) apecss_bubble_initializestruct(Bubbles[i]); + + /* Set default options and read the options for the bubble */ + 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 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)); + struct APECSS_Interface *Interface = (struct APECSS_Interface *) malloc(sizeof(struct APECSS_Interface)); + struct APECSS_NumericsODE *NumericsODE = (struct APECSS_NumericsODE *) malloc(sizeof(struct APECSS_NumericsODE)); + + /* Set the default options for the fluid properties and solver parameters */ + apecss_gas_setdefaultoptions(Gas); + apecss_liquid_setdefaultoptions(Liquid); + apecss_interface_setdefaultoptions(Interface); + apecss_odesolver_setdefaultoptions(NumericsODE); + + /* Read the options file for the fluid properties and solver parameters */ + apecss_gas_readoptions(Gas, OptionsDir); + apecss_liquid_readoptions(Liquid, OptionsDir); + apecss_interface_readoptions(Interface, OptionsDir); + apecss_odesolver_readoptions(NumericsODE, OptionsDir); + + /* Associate the bubble with the relevant fluid properties and solver parameters */ + for (register int i = 0; i < nBubbles; i++) + { + Bubbles[i]->Gas = Gas; + Bubbles[i]->Liquid = Liquid; + Bubbles[i]->Interface = Interface; + Bubbles[i]->NumericsODE = NumericsODE; + } + + /* Allocate and set the excitation parameters */ + struct APECSS_Excitation *Excitation = (struct APECSS_Excitation *) malloc(sizeof(struct APECSS_Excitation)); + Excitation->type = APECSS_EXCITATION_SIN; + Excitation->f = (APECSS_FLOAT) fa; + Excitation->dp = (APECSS_FLOAT) pa; + for (register int i = 0; i < nBubbles; i++) Bubbles[i]->Excitation = Excitation; + + // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + // Create individual folders for the results of each bubble + for (register int i = 0; i < nBubbles; i++) + { + if (Bubbles[i]->Results != NULL) + { + sprintf(Bubbles[i]->Results->dir, "./Bubble_%i/", i); + struct stat st = {0}; + if (stat(Bubbles[i]->Results->dir, &st) == -1) mkdir(Bubbles[i]->Results->dir, 0700); + } + } + // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + /* Process all options */ + apecss_gas_processoptions(Gas); + apecss_liquid_processoptions(Liquid); + apecss_interface_processoptions(Interface); + apecss_odesolver_processoptions(NumericsODE); + for (register int i = 0; i < nBubbles; i++) apecss_bubble_processoptions(Bubbles[i]); + + // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + // 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 + Bubbles[0]->R0 = 2.0e-06; + Bubbles[1]->R0 = 20.0e-06; + + // Define center location for each bubble + Bubbles[0]->Interaction->location[0] = 0.0; + Bubbles[0]->Interaction->location[1] = 0.0; + Bubbles[0]->Interaction->location[2] = 0.0; + + Bubbles[1]->Interaction->location[0] = (APECSS_FLOAT) cluster_size * (Bubbles[0]->R0 + Bubbles[1]->R0); + Bubbles[1]->Interaction->location[1] = 0.0; + Bubbles[1]->Interaction->location[2] = 0.0; + // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + clock_t starttimebubble = clock(); + APECSS_FLOAT tSim = 0.0; // Simulation time of the coupled system + + for (register int i = 0; i < nBubbles; i++) + { + Bubbles[i]->tStart = tSim; + Bubbles[i]->tEnd = (APECSS_FLOAT) tEnd; + Bubbles[i]->dt = APECSS_MIN(1.0e-11, dt_interbubble); // Initial time-step + } + + /* Initialize the bubble based on the selected options */ + for (register int i = 0; i < nBubbles; i++) apecss_bubble_initialize(Bubbles[i]); + + /* Initialize */ + for (register int i = 0; i < nBubbles; i++) apecss_bubble_solver_initialize(Bubbles[i]); + + // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + // File to retrieve all valuable information for cavitation onset test case + FILE *file_ida2009; + file_ida2009 = fopen("Ida2009_results.txt", "w"); + fprintf(file_ida2009, "%d Bubbles p0(pa) %e png(Pa) %e D_multiplier(-) %e Interaction-type %d\n", nBubbles, Liquid->pref, pa, cluster_size, 0); + fprintf(file_ida2009, "Initial_radii(m)"); + for (register int i = 0; i < nBubbles; i++) fprintf(file_ida2009, " %e", Bubbles[i]->R0); + fprintf(file_ida2009, "\n"); + fprintf(file_ida2009, "#Time(s) R(m) Pt(Pa)\n"); + // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + /* Solve the bubble dynamics */ + while (tSim < (APECSS_FLOAT) tEnd) // Interaction loop, corresponding to the time-intervals at which interactions are considered + { + APECSS_FLOAT dtSim = APECSS_MIN(dt_interbubble, (APECSS_FLOAT) tEnd - tSim); + tSim += dtSim; + + for (register int i = 0; i < nBubbles; i++) apecss_bubble_solver_run(tSim, Bubbles[i]); + + // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + // Retrieve data + fprintf(file_ida2009, "%e", tSim); + for (register int i = 0; i < nBubbles; i++) + { + fprintf(file_ida2009, " %e", Bubbles[i]->R); + } + for (register int i = 0; i < nBubbles; i++) + { + fprintf(file_ida2009, " %e", Bubbles[i]->get_pressure_infinity(tSim, Bubbles[i])); + } + fprintf(file_ida2009, "\n"); + // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + } + + fclose(file_ida2009); + + /* Finalize the simulation*/ + for (register int i = 0; i < nBubbles; i++) apecss_bubble_solver_finalize(Bubbles[i]); + + char str[APECSS_STRINGLENGTH_SPRINTF]; + for (register int i = 0; i < nBubbles; i++) + { + sprintf(str, "Bubble %i: Solver concluded %i time-steps and %i sub-iterations in %.3f s.", i, Bubbles[i]->dtNumber, Bubbles[i]->nSubIter, + (double) (clock() - starttimebubble) / CLOCKS_PER_SEC); + apecss_writeonscreen(str); + } + + for (register int i = 0; i < nBubbles; i++) apecss_results_rayleighplesset_write(Bubbles[i], APECSS_RESULTS_WRITE); + for (register int i = 0; i < nBubbles; i++) apecss_results_emissionsspace_write(Bubbles[i], APECSS_RESULTS_WRITE); + + /* Make sure all allocated memory is freed */ + for (register int i = 0; i < nBubbles; i++) apecss_bubble_freestruct(Bubbles[i]); + + for (register int i = 0; i < nBubbles; i++) free(Bubbles[i]); + free(Gas); + free(Liquid); + free(Interface); + free(NumericsODE); + free(Excitation); + + return (0); +} + +APECSS_FLOAT interaction_bubble_pressure_infinity(APECSS_FLOAT t, struct APECSS_Bubble *Bubble) +{ + APECSS_FLOAT T = 10.0e-06; + if (t < T) + { + return (Bubble->p0 + Bubble->Interaction->dp_neighbor); + } + else if (t > 2 * T) + { + return (Bubble->Excitation->dp + Bubble->Interaction->dp_neighbor); + } + else + { + APECSS_FLOAT W = 0.5 * (1 - APECSS_COS((t + T) * APECSS_PI / T)); + return (Bubble->p0 + W * (Bubble->Excitation->dp - Bubble->p0) + 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 T = 10.0e-06; + APECSS_FLOAT derivative = 0.0; + if ((t >= T) && (t <= 2 * T)) + { + APECSS_FLOAT inv_T = 1 / T; + derivative = 0.5 * APECSS_PI * inv_T * APECSS_SIN(APECSS_PI * (t + T) * inv_T) * (Bubble->Excitation->dp - Bubble->p0); + } + return (derivative); +} \ No newline at end of file diff --git a/examples/cavitationonset/QA/build/CMakeLists.txt b/examples/cavitationonset/QA/build/CMakeLists.txt new file mode 100644 index 0000000..8e43b7f --- /dev/null +++ b/examples/cavitationonset/QA/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}) \ No newline at end of file diff --git a/examples/cavitationonset/QA/build/compile.sh b/examples/cavitationonset/QA/build/compile.sh new file mode 100755 index 0000000..02b82c5 --- /dev/null +++ b/examples/cavitationonset/QA/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 \ No newline at end of file diff --git a/examples/cavitationonset/QA/gather_results.py b/examples/cavitationonset/QA/gather_results.py new file mode 100644 index 0000000..79929e0 --- /dev/null +++ b/examples/cavitationonset/QA/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() \ No newline at end of file diff --git a/examples/cavitationonset/QA/run.apecss b/examples/cavitationonset/QA/run.apecss new file mode 100644 index 0000000..2e45cda --- /dev/null +++ b/examples/cavitationonset/QA/run.apecss @@ -0,0 +1,36 @@ +######################################################### +# # +# APECSS Options File # +# # +######################################################### + +BUBBLE +RPModel KM +Emissions QA 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 diff --git a/examples/cavitationonset/QA/src/cavitationonset_apecss.c b/examples/cavitationonset/QA/src/cavitationonset_apecss.c new file mode 100644 index 0000000..a2bf36e --- /dev/null +++ b/examples/cavitationonset/QA/src/cavitationonset_apecss.c @@ -0,0 +1,336 @@ +// This source file is part of APECSS, an open-source software toolbox +// for the computation of pressure-driven bubble dynamics and acoustic +// emissions in spherical symmetry. +// +// Copyright (C) 2022-2024 The APECSS Developers +// +// The APECSS Developers are listed in the README.md file available in +// the GitHub repository at https://github.com/polycfd/apecss. +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at https://mozilla.org/MPL/2.0/. + +// ------------------------------------------------------------------- +// APECSS standalone example of cavitation onset with acoustically-in- +// teracting microbubbles, based on Ida (2009), Physics of Fluids 21 +// (11), 113302, DOI : 10.1063/1.3265547 +// ------------------------------------------------------------------- + +#include +#include "apecss.h" + +// 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); + +int main(int argc, char **args) +{ + char OptionsDir[APECSS_STRINGLENGTH]; + + // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + // Interbubble time-step, defining the frequency with which the neighbor influence is updated + APECSS_FLOAT dt_interbubble = 1.0e-8; + + // Number of bubbles + const int nBubbles = 2; + + // Initialize the simulation parameters given by the execution command + double tEnd = 0.0; + double fa = 0.0; + double pa = 0.0; + double cluster_size = 0; + // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + apecss_infoscreen(); + + /* Read commandline options */ + sprintf(OptionsDir, "./run.apecss"); // This is the default + int j = 1; // First argument is the call to the executable + while (j < argc) + { + if (strcmp("-options", args[j]) == 0) + { + sprintf(OptionsDir, "%s", args[j + 1]); + j += 2; + } + else if (strcmp("-tend", args[j]) == 0) + { + sscanf(args[j + 1], "%le", &tEnd); + j += 2; + } + else if (strcmp("-freq", args[j]) == 0) + { + sscanf(args[j + 1], "%le", &fa); + j += 2; + } + else if (strcmp("-amp", args[j]) == 0) + { + sscanf(args[j + 1], "%le", &pa); + j += 2; + } + else if (strcmp("-clsize", args[j]) == 0) + { + sscanf(args[j + 1], "%le", &cluster_size); + j += 2; + } + else if (strcmp("-dt_inter", args[j]) == 0) + { + sscanf(args[j + 1], "%le", &dt_interbubble); + j += 2; + } + else if (strcmp("-h", args[j]) == 0) + { + apecss_helpscreen(); + } + else + { + char str[APECSS_STRINGLENGTH_SPRINTF]; + sprintf(str, "Unknown command line options: %s", args[j]); + apecss_erroronscreen(1, str); + ++j; + } + } + dt_interbubble = (APECSS_FLOAT) dt_interbubble; + + /* Allocate and initialize Bubble structure */ + struct APECSS_Bubble *Bubbles[nBubbles]; + for (register int i = 0; i < nBubbles; i++) Bubbles[i] = (struct APECSS_Bubble *) malloc(nBubbles * sizeof(struct APECSS_Bubble)); + for (register int i = 0; i < nBubbles; i++) apecss_bubble_initializestruct(Bubbles[i]); + + /* Set default options and read the options for the bubble */ + 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 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)); + struct APECSS_Interface *Interface = (struct APECSS_Interface *) malloc(sizeof(struct APECSS_Interface)); + struct APECSS_NumericsODE *NumericsODE = (struct APECSS_NumericsODE *) malloc(sizeof(struct APECSS_NumericsODE)); + + /* Set the default options for the fluid properties and solver parameters */ + apecss_gas_setdefaultoptions(Gas); + apecss_liquid_setdefaultoptions(Liquid); + apecss_interface_setdefaultoptions(Interface); + apecss_odesolver_setdefaultoptions(NumericsODE); + + /* Read the options file for the fluid properties and solver parameters */ + apecss_gas_readoptions(Gas, OptionsDir); + apecss_liquid_readoptions(Liquid, OptionsDir); + apecss_interface_readoptions(Interface, OptionsDir); + apecss_odesolver_readoptions(NumericsODE, OptionsDir); + + /* Associate the bubble with the relevant fluid properties and solver parameters */ + for (register int i = 0; i < nBubbles; i++) + { + Bubbles[i]->Gas = Gas; + Bubbles[i]->Liquid = Liquid; + Bubbles[i]->Interface = Interface; + Bubbles[i]->NumericsODE = NumericsODE; + } + + /* Allocate and set the excitation parameters */ + struct APECSS_Excitation *Excitation = (struct APECSS_Excitation *) malloc(sizeof(struct APECSS_Excitation)); + Excitation->type = APECSS_EXCITATION_SIN; + Excitation->f = (APECSS_FLOAT) fa; + Excitation->dp = (APECSS_FLOAT) pa; + for (register int i = 0; i < nBubbles; i++) Bubbles[i]->Excitation = Excitation; + + // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + // Create individual folders for the results of each bubble + for (register int i = 0; i < nBubbles; i++) + { + if (Bubbles[i]->Results != NULL) + { + sprintf(Bubbles[i]->Results->dir, "./Bubble_%i/", i); + struct stat st = {0}; + if (stat(Bubbles[i]->Results->dir, &st) == -1) mkdir(Bubbles[i]->Results->dir, 0700); + } + } + // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + /* Process all options */ + apecss_gas_processoptions(Gas); + apecss_liquid_processoptions(Liquid); + apecss_interface_processoptions(Interface); + apecss_odesolver_processoptions(NumericsODE); + for (register int i = 0; i < nBubbles; i++) apecss_bubble_processoptions(Bubbles[i]); + + // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + // 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 + Bubbles[0]->R0 = 2.0e-06; + Bubbles[1]->R0 = 20.0e-06; + + // Define center location for each bubble + Bubbles[0]->Interaction->location[0] = 0.0; + Bubbles[0]->Interaction->location[1] = 0.0; + Bubbles[0]->Interaction->location[2] = 0.0; + + Bubbles[1]->Interaction->location[0] = (APECSS_FLOAT) cluster_size * (Bubbles[0]->R0 + Bubbles[1]->R0); + Bubbles[1]->Interaction->location[1] = 0.0; + Bubbles[1]->Interaction->location[2] = 0.0; + // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + clock_t starttimebubble = clock(); + APECSS_FLOAT tSim = 0.0; // Simulation time of the coupled system + + for (register int i = 0; i < nBubbles; i++) + { + Bubbles[i]->tStart = tSim; + Bubbles[i]->tEnd = (APECSS_FLOAT) tEnd; + Bubbles[i]->dt = APECSS_MIN(1.0e-11, dt_interbubble); // Initial time-step + } + + /* Initialize the bubble based on the selected options */ + for (register int i = 0; i < nBubbles; i++) apecss_bubble_initialize(Bubbles[i]); + + /* Initialize */ + for (register int i = 0; i < nBubbles; i++) apecss_bubble_solver_initialize(Bubbles[i]); + + // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + // File to retrieve all valuable information for cavitation onset test case + FILE *file_ida2009; + file_ida2009 = fopen("Ida2009_results.txt", "w"); + fprintf(file_ida2009, "%d Bubbles p0(pa) %e png(Pa) %e D_multiplier(-) %e Interaction-type %d\n", nBubbles, Liquid->pref, pa, cluster_size, 2); + fprintf(file_ida2009, "Initial_radii(m)"); + for (register int i = 0; i < nBubbles; i++) fprintf(file_ida2009, " %e", Bubbles[i]->R0); + fprintf(file_ida2009, "\n"); + fprintf(file_ida2009, "#Time(s) R(m) Pt(Pa)\n"); + // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + /* Solve the bubble dynamics */ + while (tSim < (APECSS_FLOAT) tEnd) // Interaction loop, corresponding to the time-intervals at which interactions are considered + { + APECSS_FLOAT dtSim = APECSS_MIN(dt_interbubble, (APECSS_FLOAT) tEnd - tSim); + tSim += dtSim; + + for (register int i = 0; i < nBubbles; i++) apecss_bubble_solver_run(tSim, Bubbles[i]); + + // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + // Update the contribution of the neighbor bubbles + apecss_interactions_quasi_acoustic(Bubbles); + // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + // Retrieve data + fprintf(file_ida2009, "%e", tSim); + for (register int i = 0; i < nBubbles; i++) + { + fprintf(file_ida2009, " %e", Bubbles[i]->R); + } + for (register int i = 0; i < nBubbles; i++) + { + fprintf(file_ida2009, " %e", Bubbles[i]->get_pressure_infinity(tSim, Bubbles[i])); + } + for (register int i = 0; i < nBubbles; i++) + { + fprintf(file_ida2009, " %e", Bubbles[i]->ode[0](Bubbles[i]->ODEsSol, Bubbles[i]->t, Bubbles[i])); + } + for (register int i = 0; i < nBubbles; i++) + { + fprintf(file_ida2009, " %e", Bubbles[i]->Interaction->dp_neighbor); + } + for (register int i = 0; i < nBubbles; i++) + { + fprintf(file_ida2009, " %e", Bubbles[i]->U); + } + fprintf(file_ida2009, "\n"); + // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + for (register int i = 0; i < nBubbles; i++) + { + Bubbles[i]->Interaction->last_t_2 = Bubbles[i]->Interaction->last_t_1; + Bubbles[i]->Interaction->last_p_2 = Bubbles[i]->Interaction->last_p_1; + + Bubbles[i]->Interaction->last_t_1 = tSim; + Bubbles[i]->Interaction->last_p_1 = Bubbles[i]->Interaction->dp_neighbor; + } + } + + fclose(file_ida2009); + + /* Finalize the simulation*/ + for (register int i = 0; i < nBubbles; i++) apecss_bubble_solver_finalize(Bubbles[i]); + + char str[APECSS_STRINGLENGTH_SPRINTF]; + for (register int i = 0; i < nBubbles; i++) + { + sprintf(str, "Bubble %i: Solver concluded %i time-steps and %i sub-iterations in %.3f s.", i, Bubbles[i]->dtNumber, Bubbles[i]->nSubIter, + (double) (clock() - starttimebubble) / CLOCKS_PER_SEC); + apecss_writeonscreen(str); + } + + for (register int i = 0; i < nBubbles; i++) apecss_results_rayleighplesset_write(Bubbles[i], APECSS_RESULTS_WRITE); + for (register int i = 0; i < nBubbles; i++) apecss_results_emissionsspace_write(Bubbles[i], APECSS_RESULTS_WRITE); + + /* Make sure all allocated memory is freed */ + for (register int i = 0; i < nBubbles; i++) apecss_bubble_freestruct(Bubbles[i]); + + for (register int i = 0; i < nBubbles; i++) free(Bubbles[i]); + free(Gas); + free(Liquid); + free(Interface); + free(NumericsODE); + free(Excitation); + + return (0); +} + +APECSS_FLOAT interaction_bubble_pressure_infinity(APECSS_FLOAT t, struct APECSS_Bubble *Bubble) +{ + APECSS_FLOAT T = 10.0e-06; + if (t < T) + { + return (Bubble->p0 + Bubble->Interaction->dp_neighbor); + } + else if (t > 2 * T) + { + return (Bubble->Excitation->dp + Bubble->Interaction->dp_neighbor); + } + else + { + APECSS_FLOAT W = 0.5 * (1 - APECSS_COS((t + T) * APECSS_PI / T)); + return (Bubble->p0 + W * (Bubble->Excitation->dp - Bubble->p0) + 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 T = 10.0e-06; + APECSS_FLOAT derivative = 0.0; + if ((t >= T) && (t <= 2 * T)) + { + APECSS_FLOAT inv_T = 1 / T; + derivative = 0.5 * APECSS_PI * inv_T * APECSS_SIN(APECSS_PI * (t + T) * inv_T) * (Bubble->Excitation->dp - Bubble->p0); + } + + APECSS_FLOAT delta_t = Bubble->Interaction->last_t_1 - Bubble->Interaction->last_t_2; + if (delta_t > Bubble->dt) + { + APECSS_FLOAT inv_delta_t = 1 / delta_t; + return (derivative + ((Bubble->Interaction->last_p_1 - Bubble->Interaction->last_p_2) * inv_delta_t)); + } + else + { + return (derivative); + } +} \ No newline at end of file diff --git a/examples/cavitationonset/README.md b/examples/cavitationonset/README.md new file mode 100644 index 0000000..7c1ed31 --- /dev/null +++ b/examples/cavitationonset/README.md @@ -0,0 +1 @@ +This example builds a standalone APECSS code for two interacting microbubbles, following the work of [Ida, M. (2009). Multibubble cavitation inception. Physics of Fluids, 21(11), 113302.](https://doi.org/10.1063/1.3265547), using the _standard incompressible model_ and the newly developped _quasi acoustic model_ to compute the acoustic emissions. There is three folders, representing each a specific interaction model ("NI" for "no interaction, "IC" for the _standard incompressible model_ and "QA" for the _quasi acoustic model_). Each folder is containing a provided [run.apecss](./IC/run.apecss) file to reproduce the bubble dynamics shown in Figs. 5, 6, 9, 11, 12 and 13 of the paper of Ida such as a [gather_results.py](./IC/gather_results.py) to properly keep the results. The [run_cavitationonset.sh](./run_cavitationonset.sh) file provides the whole execution commands needed to gather data in order to reproduce results from the source article. The plots are designed to exhibit the differences between the two interaction models when studying cavitation onset and collapse of the bubbles. \ No newline at end of file diff --git a/examples/cavitationonset/plot_results.py b/examples/cavitationonset/plot_results.py new file mode 100644 index 0000000..0e9fc78 --- /dev/null +++ b/examples/cavitationonset/plot_results.py @@ -0,0 +1,427 @@ +import os +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.colors as mcolors + +plt.rcParams['font.family']='serif' +plt.rcParams['font.serif']=['Times New Roman'] + plt.rcParams['font.serif'] +plt.rcParams['mathtext.fontset']='stix' +plt.rcParams['font.size']=25 + +color_names = list(mcolors.XKCD_COLORS) + +cm = 1/2.54 + +# File designed to recover and plot results in order to study cavitation onset and collapse with the two interaction models +# This test case is inspired by https://doi.org/10.1063/1.3265547 + +######### Step 1 : Recovering data ################################################################################################################## + +inttype_list = ["NI", "IC", "QA"] + +dic_2_bubbles = {} + +working_path = os.getcwd() + +for inttype in inttype_list : + if inttype not in list(dic_2_bubbles.keys()) : + dic_2_bubbles[inttype] = {} + + inttype_path = os.path.join(working_path, inttype) + inttype_path = os.path.join(inttype_path, "results") + + for file in os.listdir(inttype_path) : + file_path = os.path.join(inttype_path, file) + file_results = open(file_path, "r") + lines = file_results.readlines() + file_results.close() + + first_line = lines[0].split(" ") + count = int(first_line[0]) + png = float(first_line[5]) + size = float(first_line[7]) + + if png not in list(dic_2_bubbles[inttype].keys()) : + dic_2_bubbles[inttype][png] = {} + if size not in list(dic_2_bubbles[inttype][png].keys()) : + dic_2_bubbles[inttype][png][size] = [] + + dic_data = dic_2_bubbles[inttype][png][size] + + second_line = lines[1].split(" ") + for i in range(count) : + init_radius = float(second_line[i+1]) + # list format : for each bubble, [R0, t_list, R_list, Pt_list] + dic_data.append([init_radius, [], [], []]) + if (inttype == "IC" or inttype == "QA") : + dic_data[-1].append([]) + dic_data[-1].append([]) + dic_data[-1].append([]) + + for line in lines[3:] : + data = line.split(" ") + t = float(data[0]) + for i in range(count) : + r = float(data[1 + i]) + pt = float(data[1 + count + i]) + dic_data[i][1].append(t) + dic_data[i][2].append(r) + dic_data[i][3].append(pt) + if (inttype == "IC" or inttype == "QA") : + a = float(data[1 + 2 * count + i]) + dp = float(data[1 + 3 * count + i]) + u = float(data[1 + 4 * count + i]) + dic_data[i][4].append(a) + dic_data[i][5].append(dp) + dic_data[i][6].append(u) + +######### Step 2 : Plotting results ################################################################################################################# + +######### Initial parameters #################### + +T = 10.0e-06 +P0 = 0.1013e06 +sigma = 0.0728 + +######### Cavitation inception with interactions with varying distance between 2 bubbles ########## + +nrow = 1 +ncol = 2 + +fig, axs = plt.subplots(nrow, ncol, figsize=((ncol*20*cm, nrow*12.5*cm))) +plt.subplots_adjust(wspace=0.35*cm, hspace=0.5*cm) + +dist_list = [10, 12, 12.1, 12.5, 15, 20] + +axs[0].set_title(r"Incompressible interactions", y=1.025) +axs[0].set_xlabel(r"$t$ [$\mu$s]", fontsize=27.5) +axs[0].set_xlim(xmin=0.0, xmax=60.0) +axs[0].set_ylabel(r"$R$ [$\mu$m]", fontsize=27.5) +axs[0].set_ylim(ymin=-5.0, ymax=80.0) +axs[0].grid() + +for dist in dist_list : + t_list = np.array(dic_2_bubbles["IC"][-25325][dist][0][1]) * 1.0e6 + r_list = np.array(dic_2_bubbles["IC"][-25325][dist][0][2]) * 1.0e6 + + axs[0].plot(t_list, r_list, color="blue", linewidth=2.5) + +t_list = np.array(dic_2_bubbles["NI"][-25325][15.0][0][1]) * 1.0e6 +r_list = np.array(dic_2_bubbles["NI"][-25325][15.0][0][2]) * 1.0e6 +axs[0].plot(t_list, r_list, color="blue", linewidth=2.5) + +t_list = np.array(dic_2_bubbles["NI"][-25325][15.0][1][1]) * 1.0e6 +r_list = np.array(dic_2_bubbles["NI"][-25325][15.0][1][2]) * 1.0e6 +axs[0].plot(t_list, r_list, color="magenta", linestyle="dashed", linewidth=2.5) + +axs[0].text(0.5, 5.0, r"$R_{1,0}$", color="blue") +axs[0].text(0.5, 25.0, r"$R_{2,0}$", color="magenta") + +axs[0].text(36.5, 73.5, r"$\infty$", color="blue") +axs[0].text(42.5, 73.5, r"20", color="blue") +axs[0].text(54.0, 73.5, r"15", color="blue") +axs[0].text(47.0, 54.0, r"12.5", color="blue") +axs[0].text(45.0, 17.0, r"12.1", color="blue") +axs[0].text(41.0, 6.75, r"12", color="blue") +axs[0].text(25.0, -2.5, r"10", color="blue") + +dist_list = [10, 11.9, 12, 15, 20] + +axs[1].set_title(r"Quasi-acoustic interactions", y=1.025) +axs[1].set_xlabel(r"$t$ [$\mu$s]", fontsize=27.5) +axs[1].set_xlim(xmin=0.0, xmax=60.0) +axs[1].set_ylim(ymin=-5.0, ymax=80.0) +axs[1].grid() + +for dist in dist_list : + t_list = np.array(dic_2_bubbles["QA"][-25325][dist][0][1]) * 1.0e6 + r_list = np.array(dic_2_bubbles["QA"][-25325][dist][0][2]) * 1.0e6 + + axs[1].plot(t_list, r_list, color="blue", linewidth=2.5) + +t_list = np.array(dic_2_bubbles["NI"][-25325][15.0][0][1]) * 1.0e6 +r_list = np.array(dic_2_bubbles["NI"][-25325][15.0][0][2]) * 1.0e6 +axs[1].plot(t_list, r_list, color="blue", linewidth=2.5) + +t_list = np.array(dic_2_bubbles["NI"][-25325][15.0][1][1]) * 1.0e6 +r_list = np.array(dic_2_bubbles["NI"][-25325][15.0][1][2]) * 1.0e6 +axs[1].plot(t_list, r_list, color="magenta", linestyle="dashed", linewidth=2.5) + +axs[1].text(0.5, 5.0, r"$R_{1,0}$", color="blue") +axs[1].text(0.5, 25.0, r"$R_{2,0}$", color="magenta") + +axs[1].text(36.5, 73.5, r"$\infty$", color="blue") +axs[1].text(42.5, 73.5, r"20", color="blue") +axs[1].text(54.0, 73.5, r"15", color="blue") +axs[1].text(51.0, 36.0, r"12", color="blue") +axs[1].text(40.0, 10.0, r"11.9", color="blue") +axs[1].text(25.0, -2.5, r"10", color="blue") + +fig.savefig("cavitationonset_varyingdistance.pdf", bbox_inches='tight',pad_inches=0.35) + +######### Cavitation inception with interactions with varying pressure between 2 bubbles ########## + +nrow = 1 +ncol = 2 + +fig, axs = plt.subplots(nrow, ncol, figsize=((ncol*20*cm, nrow*12.5*cm))) +plt.subplots_adjust(wspace=0.35*cm, hspace=0.5*cm) + +png_list = [-25325, -27351, -27654.9, -27958.8, -29377] + +axs[0].set_title(r"Incompressible interactions", y=1.025) +axs[0].set_xlabel(r"$t$ [$\mu$s]", fontsize=27.5) +axs[0].set_xlim(xmin=0.0, xmax=60.0) +axs[0].set_ylabel(r"$R_{1}$ [$\mu$m]", fontsize=27.5) +axs[0].set_ylim(ymin=0.0, ymax=40.0) +axs[0].grid() + +for png in png_list : + t_list = np.array(dic_2_bubbles["IC"][png][10.0][0][1]) * 1.0e6 + r_list = np.array(dic_2_bubbles["IC"][png][10.0][0][2]) * 1.0e6 + + axs[0].plot(t_list, r_list, color="blue", linewidth=2.5) + +axs[0].text(0.5, 3.0, r"$R_{1,0}$", color="blue") + +axs[0].text(32.0, 37.5, r"-0.29", color="blue") +axs[0].text(45.0, 28.5, r"-0.276", color="blue") +axs[0].text(43.0, 17.5, r"-0.273", color="blue") +axs[0].text(30.0, 8.0, r"-0.27", color="blue") +axs[0].text(23.0, 0.5, r"-0.25", color="blue") + +png_list = [-25325, -27351, -27958.8, -27654.9, -29377] + +axs[1].set_title(r"Quasi-acoustic interactions", y=1.025) +axs[1].set_xlabel(r"$t$ [$\mu$s]", fontsize=27.5) +axs[1].set_xlim(xmin=0.0, xmax=60.0) +axs[1].set_ylim(ymin=0.0, ymax=40.0) +axs[1].grid() + +for png in png_list : + t_list = np.array(dic_2_bubbles["QA"][png][10.0][0][1]) * 1.0e6 + r_list = np.array(dic_2_bubbles["QA"][png][10.0][0][2]) * 1.0e6 + + axs[1].plot(t_list, r_list, color="blue", linewidth=2.5) + +axs[1].text(0.5, 3.0, r"$R_{1,0}$", color="blue") + +axs[1].text(30.5, 37.5, r"-0.29", color="blue") +axs[1].text(48.5, 37.5, r"-0.276", color="blue") +axs[1].text(50.0, 28.0, r"-0.273", color="blue") +axs[1].text(36.5, 8.0, r"-0.27", color="blue") +axs[1].text(25.0, 0.5, r"-0.25", color="blue") + +fig.savefig("cavitationonset_varyingpressure.pdf", bbox_inches='tight',pad_inches=0.35) + +##### Pressure differences between the two interactions models evolution ###### + +nrow = 2 +ncol = 1 + +fig, axs = plt.subplots(nrow, ncol, figsize=((ncol*20*cm, nrow*12.5*cm)), sharex=True) +plt.subplots_adjust(wspace=1.5*cm, hspace=0.25*cm) + +png_list = [-25325, -27654.9] +dist_list = [12.0, 10.0] + +for i in range(nrow) : + png = png_list[i] + dist = dist_list[i] + + if (i == 1) : + axs[i].set_xlabel(r"$t$ [$\mu$s]", fontsize=27.5) + axs[i].set_ylabel(r"$p_{\infty, 1} / p_{0}$", fontsize=27.5) + axs[i].set_xlim(xmin=10.0, xmax=60.0) + axs[i].grid() + + secyax = axs[i].twinx() + secyax.yaxis.label.set_color("blue") + secyax.spines["right"].set_color("blue") + secyax.tick_params("y", colors="blue") + secyax.spines["right"].set_edgecolor("blue") + secyax.set_ylabel(r"$(p_{\infty, 1, \mathrm{QA}} - p_{\infty, 1, \mathrm{IC}}) / p_{0}$", fontsize=27.5) + secyax.grid(linestyle="dashed") + + sec_bis_yaxis = axs[i].twinx() + sec_bis_yaxis.set_yticks([]) + + t_list_IC = np.array(dic_2_bubbles["IC"][png][dist][0][1]) * 1.0e6 + p_list_IC = np.array(dic_2_bubbles["IC"][png][dist][0][3]) / P0 + + t_list_QA = np.array(dic_2_bubbles["QA"][png][dist][0][1]) * 1.0e6 + p_list_QA = np.array(dic_2_bubbles["QA"][png][dist][0][3]) / P0 + + p_list_QA_new = [] + for j in range(9, len(p_list_QA), 10) : + p_list_QA_new.append(p_list_QA[j]) + + diff_p = np.array(p_list_QA_new) - p_list_IC + + axs[i].plot(t_list_IC, p_list_IC, color="red", linewidth=3.0, linestyle="solid", label="IC") + axs[i].plot(t_list_QA, p_list_QA, color="black", linewidth=3.0, linestyle="dashed", label="QA") + secyax.plot(t_list_IC, diff_p, color="blue", linewidth=3.0, linestyle="dotted", label=r"$\Delta p_{\infty, 1} / p_{0}$ (QA - IC)") + sec_bis_yaxis.plot(t_list_IC, p_list_IC, color="red", linewidth=3.0, linestyle="solid", label="IC") + sec_bis_yaxis.plot(t_list_QA, p_list_QA, color="black", linewidth=3.0, linestyle="dashed", label="QA") + + t_IC = t_list_IC[dic_2_bubbles["IC"][png][10.0][0][2].index(np.max(dic_2_bubbles["IC"][png][10.0][0][2]))] + t_QA = t_list_QA[dic_2_bubbles["QA"][png][10.0][0][2].index(np.max(dic_2_bubbles["QA"][png][10.0][0][2]))] + + secyax.set_ylim(ymin=-0.005, ymax=0.005) + secyax.set_yticks([-0.005, 0.0, 0.005]) + + if (i == 0) : + sec_bis_yaxis.legend(bbox_to_anchor=(0.5, 1.15), loc="center", ncol=2, frameon=False) + secyax.legend(bbox_to_anchor=(0.5, 1.05), loc="center", ncol=1, frameon=False) + + axs[i].text(25.0, 0.85, r"$\Delta x_{12}^{*}=$" + " {:.1f}".format(dist), horizontalalignment="center", fontsize=27.5) + axs[i].text(25.0, 0.65, r"$p_{\mathrm{ng}}^{*}=$" + " {:.3f}".format(png/P0), horizontalalignment="center", fontsize=27.5) + +fig.savefig("cavitationonset_pressuredifferences.pdf", bbox_inches='tight',pad_inches=0.35) + +#### Ambient pressure during collapse evolution #### + +png = -27351 +dist = 10.0 + +nrow = 2 +ncol = 2 + +fig, axs = plt.subplots(nrow, ncol, figsize=((ncol*18.75*cm, nrow*9.375*cm)), sharex=True) +plt.subplots_adjust(wspace=0.45*cm, hspace=0.25*cm) + +axs[0, 0].set_title("Incompressible interactions", y=1.025) +axs[0, 1].set_title("Quasi-acoustic interactions", y=1.025) + +for i in range(nrow) : + for j in range(ncol) : + axs[i, j].grid() + +axs[nrow - 1, 0].set_xlabel(r"$t$ [$\mu$s]") +axs[nrow - 1, 1].set_xlabel(r"$t$ [$\mu$s]") + +axs[0, 0].set_ylabel(r"$R/R_{0}$") +axs[nrow - 1, 0].set_ylabel(r"$p_{\infty} / p_{0}$") + +axs[1, 0].set_xlim(xmin=30.0, xmax=50.0) +axs[0, 0].set_ylim(ymin=0.0, ymax=10.0) +axs[0, 1].set_ylim(ymin=0.0, ymax=10.0) +axs[nrow - 1, 0].set_ylim(ymin=-0.3, ymax=0.0) +axs[nrow - 1, 1].set_ylim(ymin=-0.3, ymax=0.0) + +# IC +r0 = dic_2_bubbles["IC"][png][dist][0][0] +t_list = np.array(dic_2_bubbles["IC"][png][dist][0][1]) * 1.0e6 +r_list = np.array(dic_2_bubbles["IC"][png][dist][0][2]) / r0 +p_list = np.array(dic_2_bubbles["IC"][png][dist][0][3]) / P0 +dp_list = dic_2_bubbles["IC"][png][dist][0][5] + +axs[0, 0].plot(t_list, r_list, color="blue", linewidth=2.5, label=r"$R_{0}=2.0$ $\mu$m") +axs[nrow - 1, 0].plot(t_list, p_list, color="blue", linewidth=2.5) + +r0 = dic_2_bubbles["IC"][png][dist][1][0] +t_list = np.array(dic_2_bubbles["IC"][png][dist][1][1]) * 1.0e6 +r_list = np.array(dic_2_bubbles["IC"][png][dist][1][2]) / r0 +p_list = np.array(dic_2_bubbles["IC"][png][dist][1][3]) / P0 +dp_list = dic_2_bubbles["IC"][png][dist][1][5] + +axs[0, 0].plot(t_list, r_list, color="magenta", linestyle="dashed", linewidth=2.5, label=r"$R_{0}=20.0$ $\mu$m") +axs[nrow - 1, 0].plot(t_list, p_list, color="magenta", linestyle="dashed", linewidth=2.5) + +# QA +r0 = dic_2_bubbles["QA"][png][dist][0][0] +t_list = np.array(dic_2_bubbles["QA"][png][dist][0][1]) * 1.0e6 +r_list = np.array(dic_2_bubbles["QA"][png][dist][0][2]) / r0 +p_list = np.array(dic_2_bubbles["QA"][png][dist][0][3]) / P0 +dp_list = dic_2_bubbles["QA"][png][dist][0][5] + +axs[0, 1].plot(t_list, r_list, color="blue", linewidth=2.5) +axs[nrow - 1, 1].plot(t_list, p_list, color="blue", linewidth=2.5) + +r0 = dic_2_bubbles["QA"][png][dist][1][0] +t_list = np.array(dic_2_bubbles["QA"][png][dist][1][1]) * 1.0e6 +r_list = np.array(dic_2_bubbles["QA"][png][dist][1][2]) / r0 +p_list = np.array(dic_2_bubbles["QA"][png][dist][1][3]) / P0 +dp_list = dic_2_bubbles["QA"][png][dist][1][5] + +axs[0, 1].plot(t_list, r_list, color="magenta", linestyle="dashed", linewidth=2.5) +axs[nrow - 1, 1].plot(t_list, p_list, color="magenta", linestyle="dashed", linewidth=2.5) + +axs[0, 0].legend(bbox_to_anchor=(1.05, 1.05), loc="lower center", ncol=2, frameon=False, fontsize=27.5) +fig.savefig("cavitationonset_collapse.pdf", bbox_inches='tight',pad_inches=0.35) + +#### Unstable radius #### + +nrow = 2 +ncol = 2 + +fig, axs = plt.subplots(nrow, ncol, figsize=((ncol*15*cm, nrow*9.375*cm)), sharey=True, sharex=True) +plt.subplots_adjust(wspace=0.14*cm, hspace=0.25*cm) + +for i in range(nrow) : + for j in range(ncol) : + axs[i, j].grid() + if j == 0 : + axs[i, j].set_ylabel(r"$R$ [$\mu$m]", fontsize=27.5) + axs[i, j].set_ylim(ymin=0.0, ymax=40.0) + if i == nrow - 1 : + axs[i, j].set_xlabel(r"$t$ [$\mu$s]", fontsize=27.5) + axs[i, j].set_xlim(xmin=0.0, xmax=60.0) + axs[i, j].set_xticks([10, 30, 50], [10, 30, 50]) + +axs[0, 0].set_title(r"Incompressible interactions", y=1.025) +axs[0, 1].set_title(r"Quasi-acoustic interactions", y=1.025) + +png_list = [-25325, -27654.9] +dist_list = [12.0, 10.0] + +for l in range(nrow) : + png = png_list[l] + dist = dist_list[l] + + # IC + t_list = np.array(dic_2_bubbles["IC"][png][dist][0][1]) + r_list = np.array(dic_2_bubbles["IC"][png][dist][0][2]) + p_list = np.array(dic_2_bubbles["IC"][png][dist][0][3]) + + axs[l, 0].text(31.0, 30.0, r"$p_{\mathrm{ng}}^{*}=$" + " {:.3f}".format(png/P0), horizontalalignment="center", fontsize=26) + axs[l, 0].text(31.0, 35.0, r"$\Delta x_{12}^{*}=$" + " {:.1f}".format(dist), horizontalalignment="center", fontsize=26) + + R0 = r_list[0] + pG0 = P0 + (2 * sigma / R0) + + r_unstable_list = [] + for i in range(len(t_list)) : + p = np.polynomial.Polynomial([-pG0 * (R0**3), 0, 2 * sigma, p_list[i]]) + roots = p.roots() + r_unstable_list.append(float(np.real(roots[-1]))) + + if (l == 0) : + axs[l, 0].plot(t_list*1e6, r_list*1e6, color="blue", linewidth=3.0, label=r"$R_{1}$") + axs[l, 0].plot(t_list*1e6, np.array(r_unstable_list)*1e6, linestyle="solid", color="red", marker="o", markevery=500, linewidth=2.5, label=r"$R_{\mathrm{Ue}}$") + else : + axs[l, 0].plot(t_list*1e6, r_list*1e6, color="blue", linewidth=3.0) + axs[l, 0].plot(t_list*1e6, np.array(r_unstable_list)*1e6, linestyle="solid", color="red", marker="o", markevery=500, linewidth=2.5) + + # QA + t_list = np.array(dic_2_bubbles["QA"][png][dist][0][1]) + r_list = np.array(dic_2_bubbles["QA"][png][dist][0][2]) + p_list = np.array(dic_2_bubbles["QA"][png][dist][0][3]) + + axs[l, 1].text(31.0, 30.0, r"$p_{\mathrm{ng}}^{*}=$" + " {:.3f}".format(png/P0), horizontalalignment="center", fontsize=26) + axs[l, 1].text(31.0, 35.0, r"$\Delta x_{12}^{*}=$" + " {:.1f}".format(dist), horizontalalignment="center", fontsize=26) + + R0 = r_list[0] + pG0 = P0 + (2 * sigma / R0) + + r_unstable_list = [] + for i in range(len(t_list)) : + p = np.polynomial.Polynomial([-pG0 * (R0**3), 0, 2 * sigma, p_list[i]]) + roots = p.roots() + r_unstable_list.append(float(np.real(roots[-1]))) + + axs[l, 1].plot(t_list*1e6, r_list*1e6, color="blue", linewidth=3.0) + axs[l, 1].plot(t_list*1e6, np.array(r_unstable_list)*1e6, linestyle="solid", color="red", marker="o", markevery=5000, linewidth=2.5) + +axs[0, 0].legend(bbox_to_anchor=(1.05, 1.05), loc="lower center", ncol=2, frameon=False, fontsize=27.5) +fig.savefig("cavitationonset_unstableradius.pdf", bbox_inches='tight',pad_inches=0.35) \ No newline at end of file diff --git a/examples/cavitationonset/run_cavitationonset.sh b/examples/cavitationonset/run_cavitationonset.sh new file mode 100755 index 0000000..64cab53 --- /dev/null +++ b/examples/cavitationonset/run_cavitationonset.sh @@ -0,0 +1,105 @@ +######### No Interaction computations ############################################################################################################### + +cd NI/build +./compile.sh +cd .. + +######### Pressure time history & Radius evolution without interaction ############################ +./build/cavitationonset_apecss -options run.apecss -amp -25325 -tend 60.0e-6 -clsize 15 +python3 gather_results.py + +cd .. + +echo "" +echo "No interaction test cases passed" +echo "" + +######### Incompressible computations ############################################################################################################### + +cd IC/build +./compile.sh +cd .. + +######### Cavitation inception with interactions with varying distance between 2 bubbles ########## +dist_list=(10 11 12 12.05 12.1 12.5 15 20) +for d in "${dist_list[@]}" +do + ./build/cavitationonset_apecss -options run.apecss -amp -25325 -tend 60.0e-6 -clsize $d + python3 gather_results.py +done + +######### Cavitation inception with interactions with varying pressure between 2 bubbles ########## +png_list=(-25325 -27351 -27654.9 -27958.8 -29377) +for png in "${png_list[@]}" +do + ./build/cavitationonset_apecss -options run.apecss -amp $png -tend 60.0e-6 -clsize 10 + python3 gather_results.py +done + +cd .. + +echo "" +echo "Incompressible test cases passed" +echo "" + +######### Quasi acoustic computations ############################################################################################################### + +cd QA/build +./compile.sh +cd .. + +######### Cavitation inception with interactions with varying distance between 2 bubbles ########## +dist_list=(10 11.9 12 15 20) +for d in "${dist_list[@]}" +do + ./build/cavitationonset_apecss -options run.apecss -amp -25325 -tend 60.0e-6 -clsize $d -dt_inter 1.0e-09 + python3 gather_results.py +done + +######### Cavitation inception with interactions with varying pressure between 2 bubbles ########## +png_list=(-25325 -27351 -27654.9 -27958.8 -29377) +for png in "${png_list[@]}" +do + ./build/cavitationonset_apecss -options run.apecss -amp $png -tend 60.0e-6 -clsize 10 -dt_inter 1.0e-09 + python3 gather_results.py +done + +cd .. + +echo "" +echo "Quasi acoustic test cases passed" +echo "" + +######### Plotting results ########################################################################################################################## + +python3 plot_results.py + +######### Cleaning ################################################################################################################################## + +cd NI +rm -rf "Ida2009_results.txt" +for ((c=0; c<2; c++)) +do + rm -rf Bubble_$c +done +cd .. + +cd IC +rm -rf "Ida2009_results.txt" +for ((c=0; c<8; c++)) +do + rm -rf Bubble_$c +done +cd .. + +cd QA +rm -rf "Ida2009_results.txt" +for ((c=0; c<8; c++)) +do + rm -rf Bubble_$c +done +cd .. + +echo "" +echo "Cleaning completed" +echo "" \ No newline at end of file diff --git a/examples/run_all.sh b/examples/run_all.sh index 53aee26..94f0fa9 100755 --- a/examples/run_all.sh +++ b/examples/run_all.sh @@ -15,6 +15,8 @@ cd .. ./build/binaryinteraction_apecss -options run.apecss -freq 15.7e3 -amp -120e3 -tend 7.5e-4 python3 plot_result.py rm -r Bubble_0 Bubble_1 +cd ../cavitationonset +./run_cavitationonset.sh cd ../gastemperature/build ./compile.sh cd .. diff --git a/examples/sphericalclustertensionwave/IC/build/CMakeLists.txt b/examples/sphericalclustertensionwave/IC/build/CMakeLists.txt new file mode 100644 index 0000000..254078c --- /dev/null +++ b/examples/sphericalclustertensionwave/IC/build/CMakeLists.txt @@ -0,0 +1,41 @@ +cmake_minimum_required(VERSION 3.12) + +project (sphericalclustertensionwave_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(sphericalclustertensionwave_apecss ${myfiles}) +target_link_libraries(sphericalclustertensionwave_apecss ${mylibs}) \ No newline at end of file diff --git a/examples/sphericalclustertensionwave/IC/build/compile.sh b/examples/sphericalclustertensionwave/IC/build/compile.sh new file mode 100755 index 0000000..564c7a4 --- /dev/null +++ b/examples/sphericalclustertensionwave/IC/build/compile.sh @@ -0,0 +1,5 @@ +#!/bin/bash +rm sphericalclustertensionwave_apecss +cmake CMakeLists.txt -DCMAKE_BUILD_TYPE=Release +make +rm -r CMakeCache.txt CMakeFiles Makefile cmake_install.cmake \ No newline at end of file diff --git a/examples/sphericalclustertensionwave/IC/gather_results.py b/examples/sphericalclustertensionwave/IC/gather_results.py new file mode 100644 index 0000000..61f7914 --- /dev/null +++ b/examples/sphericalclustertensionwave/IC/gather_results.py @@ -0,0 +1,38 @@ +import os + +# File designed to recover data for plotting results in cavitation onset case + +file = open("tension_results.txt", "r") +lines = file.readlines() +file.close() + +count = int(lines[0].split(" ")[0]) +p1 = float(lines[0].split(" ")[5]) + +file_name = "mono_{}_{:.4E}.txt".format(count, p1) + +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") + +index = 0 +while index < len(lines) and "nan" not in lines[index] : + file_results.write(lines[index]) + index += 1 + +file_results.close() + +file_loc = open("bubble_loc.txt", "r") +lines_loc = file_loc.readlines() +file_loc.close() + +file_name_loc = "mono_{}_{:.4E}_loc.txt".format(count, p1) + +file_loc_results = open(os.path.join(results_path, file_name_loc), "w") + +for line in lines_loc : + file_loc_results.write(line) + +file_loc_results.close() \ No newline at end of file diff --git a/examples/sphericalclustertensionwave/IC/run.apecss b/examples/sphericalclustertensionwave/IC/run.apecss new file mode 100644 index 0000000..44bbe03 --- /dev/null +++ b/examples/sphericalclustertensionwave/IC/run.apecss @@ -0,0 +1,37 @@ +######################################################### +# # +# APECSS Options File # +# # +######################################################### + +BUBBLE +RPModel KM +Emissions IC 300.0e-6 +PressureAmbient 1.0e5 +END + +GAS +EoS IG +ReferenceDensity 1.2 +PolytropicExponent 1.4 +END + +LIQUID +ReferencePressure 1.0e5 +ReferenceDensity 1000.0 +ReferenceSoundSpeed 1500.0 +Viscosity 1.002e-3 +END + +INTERFACE +SurfaceTensionCoeff 0.0728 +END + +RESULTS +Bubble +END + +ODESOLVER +tolerance 1.0e-10 +MinTimeStep 1.0e-10 +END diff --git a/examples/sphericalclustertensionwave/IC/src/sphericalclustertensionwave_apecss.c b/examples/sphericalclustertensionwave/IC/src/sphericalclustertensionwave_apecss.c new file mode 100644 index 0000000..fbe0daf --- /dev/null +++ b/examples/sphericalclustertensionwave/IC/src/sphericalclustertensionwave_apecss.c @@ -0,0 +1,345 @@ +// This source file is part of APECSS, an open-source software toolbox +// for the computation of pressure-driven bubble dynamics and acoustic +// emissions in spherical symmetry. +// +// Copyright (C) 2022-2024 The APECSS Developers +// +// The APECSS Developers are listed in the README.md file available in +// the GitHub repository at https://github.com/polycfd/apecss. +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at https://mozilla.org/MPL/2.0/. + +// ------------------------------------------------------------------- +// APECSS standalone example of acoustically-interacting microbubbles +// in a spherical cluster with the goal to study cavitation onset +// ------------------------------------------------------------------- + +#include +#include "apecss.h" + +APECSS_FLOAT rand_range(double min, double max) +{ + double random = (drand48()); + double range = (max - min) * random; + double number = min + range; + + return (APECSS_FLOAT) number; +} + +// 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); + +int main(int argc, char **args) +{ + char OptionsDir[APECSS_STRINGLENGTH]; + + // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + // Set the case-dependent simulation parameters + const int nBubbles = 250; // Number of bubbles + APECSS_FLOAT bubble_bubble_dist = 20.0e-6; // Bubble-bubble minimal distance + APECSS_FLOAT cluster_radius = 232e-6; // Spherical cluster radius + + // Interbubble time-step, defining the frequency with which the neighbor influence is updated + APECSS_FLOAT dt_interbubble = 1.0e-8; + + // Initialize the simulation parameters given by the execution command + double tEnd = 0.0; + double fa = 0.0; + double pa = 0.0; + // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + apecss_infoscreen(); + + /* Read commandline options */ + sprintf(OptionsDir, "./run.apecss"); // This is the default + int j = 1; // First argument is the call to the executable + while (j < argc) + { + if (strcmp("-options", args[j]) == 0) + { + sprintf(OptionsDir, "%s", args[j + 1]); + j += 2; + } + else if (strcmp("-tend", args[j]) == 0) + { + sscanf(args[j + 1], "%le", &tEnd); + j += 2; + } + else if (strcmp("-freq", args[j]) == 0) + { + sscanf(args[j + 1], "%le", &fa); + j += 2; + } + else if (strcmp("-amp", args[j]) == 0) + { + sscanf(args[j + 1], "%le", &pa); + j += 2; + } + else if (strcmp("-h", args[j]) == 0) + { + apecss_helpscreen(); + } + else + { + char str[APECSS_STRINGLENGTH_SPRINTF]; + sprintf(str, "Unknown command line options: %s", args[j]); + apecss_erroronscreen(1, str); + ++j; + } + } + + /* Allocate and initialize Bubble structure */ + struct APECSS_Bubble *Bubbles[nBubbles]; + for (register int i = 0; i < nBubbles; i++) Bubbles[i] = (struct APECSS_Bubble *) malloc(nBubbles * sizeof(struct APECSS_Bubble)); + for (register int i = 0; i < nBubbles; i++) apecss_bubble_initializestruct(Bubbles[i]); + + /* Set default options and read the options for the bubble */ + 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 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)); + struct APECSS_Interface *Interface = (struct APECSS_Interface *) malloc(sizeof(struct APECSS_Interface)); + struct APECSS_NumericsODE *NumericsODE = (struct APECSS_NumericsODE *) malloc(sizeof(struct APECSS_NumericsODE)); + + /* Set the default options for the fluid properties and solver parameters */ + apecss_gas_setdefaultoptions(Gas); + apecss_liquid_setdefaultoptions(Liquid); + apecss_interface_setdefaultoptions(Interface); + apecss_odesolver_setdefaultoptions(NumericsODE); + + /* Read the options file for the fluid properties and solver parameters */ + apecss_gas_readoptions(Gas, OptionsDir); + apecss_liquid_readoptions(Liquid, OptionsDir); + apecss_interface_readoptions(Interface, OptionsDir); + apecss_odesolver_readoptions(NumericsODE, OptionsDir); + + /* Associate the bubble with the relevant fluid properties and solver parameters */ + for (register int i = 0; i < nBubbles; i++) + { + Bubbles[i]->Gas = Gas; + Bubbles[i]->Liquid = Liquid; + Bubbles[i]->Interface = Interface; + Bubbles[i]->NumericsODE = NumericsODE; + } + + /* Allocate and set the excitation parameters */ + struct APECSS_Excitation *Excitation = (struct APECSS_Excitation *) malloc(sizeof(struct APECSS_Excitation)); + Excitation->type = APECSS_EXCITATION_SIN; + Excitation->f = (APECSS_FLOAT) fa; + Excitation->dp = (APECSS_FLOAT) pa; + for (register int i = 0; i < nBubbles; i++) Bubbles[i]->Excitation = Excitation; + + // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + // Create individual folders for the results of each bubble + for (register int i = 0; i < nBubbles; i++) + { + if (Bubbles[i]->Results != NULL) + { + sprintf(Bubbles[i]->Results->dir, "./Bubble_%i/", i); + struct stat st = {0}; + if (stat(Bubbles[i]->Results->dir, &st) == -1) mkdir(Bubbles[i]->Results->dir, 0700); + } + } + // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + /* Process all options */ + apecss_gas_processoptions(Gas); + apecss_liquid_processoptions(Liquid); + apecss_interface_processoptions(Interface); + apecss_odesolver_processoptions(NumericsODE); + for (register int i = 0; i < nBubbles; i++) apecss_bubble_processoptions(Bubbles[i]); + + // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + // 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++) + { + Bubbles[i]->R0 = 2.0e-6; + } + + // Define center location for each bubble + for (register int i = 0; i < nBubbles; i++) + { + if (i == 0) + { + Bubbles[i]->Interaction->location[0] = 0.0; + Bubbles[i]->Interaction->location[1] = 0.0; + Bubbles[i]->Interaction->location[2] = 0.0; + } + else + { + APECSS_FLOAT x = rand_range((double) -cluster_radius, (double) cluster_radius); + APECSS_FLOAT y = rand_range((double) -cluster_radius, (double) cluster_radius); + APECSS_FLOAT z = rand_range((double) -cluster_radius, (double) cluster_radius); + + APECSS_FLOAT radius = APECSS_SQRT(APECSS_POW2(x) + APECSS_POW2(y) + APECSS_POW2(z)); + + while (radius > cluster_radius) + { + x = rand_range((double) -cluster_radius, (double) cluster_radius); + y = rand_range((double) -cluster_radius, (double) cluster_radius); + z = rand_range((double) -cluster_radius, (double) cluster_radius); + + radius = APECSS_SQRT(APECSS_POW2(x) + APECSS_POW2(y) + APECSS_POW2(z)); + } + + Bubbles[i]->Interaction->location[0] = x; + Bubbles[i]->Interaction->location[1] = y; + Bubbles[i]->Interaction->location[2] = z; + + for (register int k = 0; k < i; k++) + { + APECSS_FLOAT bubbledist = APECSS_SQRT(APECSS_POW2(Bubbles[i]->Interaction->location[0] - Bubbles[k]->Interaction->location[0]) + + APECSS_POW2(Bubbles[i]->Interaction->location[1] - Bubbles[k]->Interaction->location[1]) + + APECSS_POW2(Bubbles[i]->Interaction->location[2] - Bubbles[k]->Interaction->location[2])); + if (bubbledist < bubble_bubble_dist) + { + i--; + break; + } + } + } + } + // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + clock_t starttimebubble = clock(); + APECSS_FLOAT tSim = 0.0; // Simulation time of the coupled system + + for (register int i = 0; i < nBubbles; i++) + { + Bubbles[i]->tStart = tSim; + Bubbles[i]->tEnd = (APECSS_FLOAT) tEnd; + Bubbles[i]->dt = APECSS_MIN(1.0e-11, dt_interbubble); // Initial time-step + } + + /* Initialize the bubble based on the selected options */ + for (register int i = 0; i < nBubbles; i++) apecss_bubble_initialize(Bubbles[i]); + + /* Initialize */ + for (register int i = 0; i < nBubbles; i++) apecss_bubble_solver_initialize(Bubbles[i]); + + // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + // Files to retrieve all valuable information for cavitation onset test case + FILE *file_loc; + file_loc = fopen("bubble_loc.txt", "w"); + fprintf(file_loc, "#number x(m) y(m) z(m)\n"); + for (register int i = 0; i < nBubbles; i++) + { + fprintf(file_loc, "%d %e %e %e\n", i, Bubbles[i]->Interaction->location[0], Bubbles[i]->Interaction->location[1], Bubbles[i]->Interaction->location[2]); + } + fclose(file_loc); + + FILE *file_tension; + file_tension = fopen("tension_results.txt", "w"); + fprintf(file_tension, "%d Bubbles p0(pa) %e p1(Pa) %e\n", nBubbles, Liquid->pref, pa); + fprintf(file_tension, "Initial_radii(m)"); + for (register int i = 0; i < nBubbles; i++) fprintf(file_tension, " %e", Bubbles[i]->R0); + fprintf(file_tension, "\n"); + fprintf(file_tension, "#Time(s) R(m) Pt(Pa)\n"); + // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + /* Solve the bubble dynamics */ + while (tSim < (APECSS_FLOAT) tEnd) // Interaction loop, corresponding to the time-intervals at which interactions are considered + { + APECSS_FLOAT dtSim = APECSS_MIN(dt_interbubble, (APECSS_FLOAT) tEnd - tSim); + tSim += dtSim; + + for (register int i = 0; i < nBubbles; i++) apecss_bubble_solver_run(tSim, Bubbles[i]); + + // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + // Update the contribution of the neighbor bubble + apecss_interactions_instantaneous(Bubbles); + // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + // Retrieve data + fprintf(file_tension, "%e", tSim); + for (register int i = 0; i < nBubbles; i++) + { + fprintf(file_tension, " %e", Bubbles[i]->R); + } + for (register int i = 0; i < nBubbles; i++) + { + fprintf(file_tension, " %e", Bubbles[i]->get_pressure_infinity(Bubbles[i]->t, Bubbles[i])); + } + fprintf(file_tension, "\n"); + // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + } + + fclose(file_tension); + + /* Finalize the simulation*/ + for (register int i = 0; i < nBubbles; i++) apecss_bubble_solver_finalize(Bubbles[i]); + + char str[APECSS_STRINGLENGTH_SPRINTF]; + for (register int i = 0; i < nBubbles; i++) + { + sprintf(str, "Bubble %i: Solver concluded %i time-steps and %i sub-iterations in %.3f s.", i, Bubbles[i]->dtNumber, Bubbles[i]->nSubIter, + (double) (clock() - starttimebubble) / CLOCKS_PER_SEC); + apecss_writeonscreen(str); + } + + for (register int i = 0; i < nBubbles; i++) apecss_results_rayleighplesset_write(Bubbles[i], APECSS_RESULTS_WRITE); + for (register int i = 0; i < nBubbles; i++) apecss_results_emissionsspace_write(Bubbles[i], APECSS_RESULTS_WRITE); + + /* Make sure all allocated memory is freed */ + for (register int i = 0; i < nBubbles; i++) apecss_bubble_freestruct(Bubbles[i]); + + for (register int i = 0; i < nBubbles; i++) free(Bubbles[i]); + free(Gas); + free(Liquid); + free(Interface); + free(NumericsODE); + free(Excitation); + + return (0); +} + +APECSS_FLOAT interaction_bubble_pressure_infinity(APECSS_FLOAT t, struct APECSS_Bubble *Bubble) +{ + APECSS_FLOAT tau = 1.75e-6; + if (t < tau) + { + return (Bubble->p0 - (Bubble->p0 - Bubble->Excitation->dp) * APECSS_POW2(APECSS_SIN(APECSS_PI * t / tau)) + Bubble->Interaction->dp_neighbor); + } + else + { + return (Bubble->p0 + Bubble->Interaction->dp_neighbor); + } +} + +APECSS_FLOAT interaction_bubble_pressurederivative_infinity(APECSS_FLOAT t, struct APECSS_Bubble *Bubble) +{ + APECSS_FLOAT tau = 1.75e-6; + if (t < tau) + { + APECSS_FLOAT inv_tau = 1 / tau; + return (-2.0 * APECSS_PI * inv_tau * (Bubble->p0 - Bubble->Excitation->dp) * APECSS_COS(APECSS_PI * t * inv_tau) * APECSS_SIN(APECSS_PI * t * inv_tau)); + } + else + { + return (0.0); + } +} \ No newline at end of file diff --git a/examples/sphericalclustertensionwave/NI/build/CMakeLists.txt b/examples/sphericalclustertensionwave/NI/build/CMakeLists.txt new file mode 100644 index 0000000..254078c --- /dev/null +++ b/examples/sphericalclustertensionwave/NI/build/CMakeLists.txt @@ -0,0 +1,41 @@ +cmake_minimum_required(VERSION 3.12) + +project (sphericalclustertensionwave_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(sphericalclustertensionwave_apecss ${myfiles}) +target_link_libraries(sphericalclustertensionwave_apecss ${mylibs}) \ No newline at end of file diff --git a/examples/sphericalclustertensionwave/NI/build/compile.sh b/examples/sphericalclustertensionwave/NI/build/compile.sh new file mode 100755 index 0000000..564c7a4 --- /dev/null +++ b/examples/sphericalclustertensionwave/NI/build/compile.sh @@ -0,0 +1,5 @@ +#!/bin/bash +rm sphericalclustertensionwave_apecss +cmake CMakeLists.txt -DCMAKE_BUILD_TYPE=Release +make +rm -r CMakeCache.txt CMakeFiles Makefile cmake_install.cmake \ No newline at end of file diff --git a/examples/sphericalclustertensionwave/NI/gather_results.py b/examples/sphericalclustertensionwave/NI/gather_results.py new file mode 100644 index 0000000..61f7914 --- /dev/null +++ b/examples/sphericalclustertensionwave/NI/gather_results.py @@ -0,0 +1,38 @@ +import os + +# File designed to recover data for plotting results in cavitation onset case + +file = open("tension_results.txt", "r") +lines = file.readlines() +file.close() + +count = int(lines[0].split(" ")[0]) +p1 = float(lines[0].split(" ")[5]) + +file_name = "mono_{}_{:.4E}.txt".format(count, p1) + +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") + +index = 0 +while index < len(lines) and "nan" not in lines[index] : + file_results.write(lines[index]) + index += 1 + +file_results.close() + +file_loc = open("bubble_loc.txt", "r") +lines_loc = file_loc.readlines() +file_loc.close() + +file_name_loc = "mono_{}_{:.4E}_loc.txt".format(count, p1) + +file_loc_results = open(os.path.join(results_path, file_name_loc), "w") + +for line in lines_loc : + file_loc_results.write(line) + +file_loc_results.close() \ No newline at end of file diff --git a/examples/sphericalclustertensionwave/NI/run.apecss b/examples/sphericalclustertensionwave/NI/run.apecss new file mode 100644 index 0000000..44bbe03 --- /dev/null +++ b/examples/sphericalclustertensionwave/NI/run.apecss @@ -0,0 +1,37 @@ +######################################################### +# # +# APECSS Options File # +# # +######################################################### + +BUBBLE +RPModel KM +Emissions IC 300.0e-6 +PressureAmbient 1.0e5 +END + +GAS +EoS IG +ReferenceDensity 1.2 +PolytropicExponent 1.4 +END + +LIQUID +ReferencePressure 1.0e5 +ReferenceDensity 1000.0 +ReferenceSoundSpeed 1500.0 +Viscosity 1.002e-3 +END + +INTERFACE +SurfaceTensionCoeff 0.0728 +END + +RESULTS +Bubble +END + +ODESOLVER +tolerance 1.0e-10 +MinTimeStep 1.0e-10 +END diff --git a/examples/sphericalclustertensionwave/NI/src/sphericalclustertensionwave_apecss.c b/examples/sphericalclustertensionwave/NI/src/sphericalclustertensionwave_apecss.c new file mode 100644 index 0000000..0f1cde0 --- /dev/null +++ b/examples/sphericalclustertensionwave/NI/src/sphericalclustertensionwave_apecss.c @@ -0,0 +1,340 @@ +// This source file is part of APECSS, an open-source software toolbox +// for the computation of pressure-driven bubble dynamics and acoustic +// emissions in spherical symmetry. +// +// Copyright (C) 2022-2024 The APECSS Developers +// +// The APECSS Developers are listed in the README.md file available in +// the GitHub repository at https://github.com/polycfd/apecss. +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at https://mozilla.org/MPL/2.0/. + +// ------------------------------------------------------------------- +// APECSS standalone example of acoustically-interacting microbubbles +// in a spherical cluster with the goal to study cavitation onset +// ------------------------------------------------------------------- + +#include +#include "apecss.h" + +APECSS_FLOAT rand_range(double min, double max) +{ + double random = (drand48()); + double range = (max - min) * random; + double number = min + range; + + return (APECSS_FLOAT) number; +} + +// 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); + +int main(int argc, char **args) +{ + char OptionsDir[APECSS_STRINGLENGTH]; + + // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + // Set the case-dependent simulation parameters + const int nBubbles = 250; // Number of bubbles + APECSS_FLOAT bubble_bubble_dist = 20.0e-6; // Bubble-bubble minimal distance + APECSS_FLOAT cluster_radius = 232e-6; // Spherical cluster radius + + // Interbubble time-step, defining the frequency with which the neighbor influence is updated + APECSS_FLOAT dt_interbubble = 1.0e-8; + + // Initialize the simulation parameters given by the execution command + double tEnd = 0.0; + double fa = 0.0; + double pa = 0.0; + // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + apecss_infoscreen(); + + /* Read commandline options */ + sprintf(OptionsDir, "./run.apecss"); // This is the default + int j = 1; // First argument is the call to the executable + while (j < argc) + { + if (strcmp("-options", args[j]) == 0) + { + sprintf(OptionsDir, "%s", args[j + 1]); + j += 2; + } + else if (strcmp("-tend", args[j]) == 0) + { + sscanf(args[j + 1], "%le", &tEnd); + j += 2; + } + else if (strcmp("-freq", args[j]) == 0) + { + sscanf(args[j + 1], "%le", &fa); + j += 2; + } + else if (strcmp("-amp", args[j]) == 0) + { + sscanf(args[j + 1], "%le", &pa); + j += 2; + } + else if (strcmp("-h", args[j]) == 0) + { + apecss_helpscreen(); + } + else + { + char str[APECSS_STRINGLENGTH_SPRINTF]; + sprintf(str, "Unknown command line options: %s", args[j]); + apecss_erroronscreen(1, str); + ++j; + } + } + + /* Allocate and initialize Bubble structure */ + struct APECSS_Bubble *Bubbles[nBubbles]; + for (register int i = 0; i < nBubbles; i++) Bubbles[i] = (struct APECSS_Bubble *) malloc(nBubbles * sizeof(struct APECSS_Bubble)); + for (register int i = 0; i < nBubbles; i++) apecss_bubble_initializestruct(Bubbles[i]); + + /* Set default options and read the options for the bubble */ + 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 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)); + struct APECSS_Interface *Interface = (struct APECSS_Interface *) malloc(sizeof(struct APECSS_Interface)); + struct APECSS_NumericsODE *NumericsODE = (struct APECSS_NumericsODE *) malloc(sizeof(struct APECSS_NumericsODE)); + + /* Set the default options for the fluid properties and solver parameters */ + apecss_gas_setdefaultoptions(Gas); + apecss_liquid_setdefaultoptions(Liquid); + apecss_interface_setdefaultoptions(Interface); + apecss_odesolver_setdefaultoptions(NumericsODE); + + /* Read the options file for the fluid properties and solver parameters */ + apecss_gas_readoptions(Gas, OptionsDir); + apecss_liquid_readoptions(Liquid, OptionsDir); + apecss_interface_readoptions(Interface, OptionsDir); + apecss_odesolver_readoptions(NumericsODE, OptionsDir); + + /* Associate the bubble with the relevant fluid properties and solver parameters */ + for (register int i = 0; i < nBubbles; i++) + { + Bubbles[i]->Gas = Gas; + Bubbles[i]->Liquid = Liquid; + Bubbles[i]->Interface = Interface; + Bubbles[i]->NumericsODE = NumericsODE; + } + + /* Allocate and set the excitation parameters */ + struct APECSS_Excitation *Excitation = (struct APECSS_Excitation *) malloc(sizeof(struct APECSS_Excitation)); + Excitation->type = APECSS_EXCITATION_SIN; + Excitation->f = (APECSS_FLOAT) fa; + Excitation->dp = (APECSS_FLOAT) pa; + for (register int i = 0; i < nBubbles; i++) Bubbles[i]->Excitation = Excitation; + + // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + // Create individual folders for the results of each bubble + for (register int i = 0; i < nBubbles; i++) + { + if (Bubbles[i]->Results != NULL) + { + sprintf(Bubbles[i]->Results->dir, "./Bubble_%i/", i); + struct stat st = {0}; + if (stat(Bubbles[i]->Results->dir, &st) == -1) mkdir(Bubbles[i]->Results->dir, 0700); + } + } + // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + /* Process all options */ + apecss_gas_processoptions(Gas); + apecss_liquid_processoptions(Liquid); + apecss_interface_processoptions(Interface); + apecss_odesolver_processoptions(NumericsODE); + for (register int i = 0; i < nBubbles; i++) apecss_bubble_processoptions(Bubbles[i]); + + // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + // 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++) + { + Bubbles[i]->R0 = 2.0e-6; + } + + // Define center location for each bubble + for (register int i = 0; i < nBubbles; i++) + { + if (i == 0) + { + Bubbles[i]->Interaction->location[0] = 0.0; + Bubbles[i]->Interaction->location[1] = 0.0; + Bubbles[i]->Interaction->location[2] = 0.0; + } + else + { + APECSS_FLOAT x = rand_range((double) -cluster_radius, (double) cluster_radius); + APECSS_FLOAT y = rand_range((double) -cluster_radius, (double) cluster_radius); + APECSS_FLOAT z = rand_range((double) -cluster_radius, (double) cluster_radius); + + APECSS_FLOAT radius = APECSS_SQRT(APECSS_POW2(x) + APECSS_POW2(y) + APECSS_POW2(z)); + + while (radius > cluster_radius) + { + x = rand_range((double) -cluster_radius, (double) cluster_radius); + y = rand_range((double) -cluster_radius, (double) cluster_radius); + z = rand_range((double) -cluster_radius, (double) cluster_radius); + + radius = APECSS_SQRT(APECSS_POW2(x) + APECSS_POW2(y) + APECSS_POW2(z)); + } + + Bubbles[i]->Interaction->location[0] = x; + Bubbles[i]->Interaction->location[1] = y; + Bubbles[i]->Interaction->location[2] = z; + + for (register int k = 0; k < i; k++) + { + APECSS_FLOAT bubbledist = APECSS_SQRT(APECSS_POW2(Bubbles[i]->Interaction->location[0] - Bubbles[k]->Interaction->location[0]) + + APECSS_POW2(Bubbles[i]->Interaction->location[1] - Bubbles[k]->Interaction->location[1]) + + APECSS_POW2(Bubbles[i]->Interaction->location[2] - Bubbles[k]->Interaction->location[2])); + if (bubbledist < bubble_bubble_dist) + { + i--; + break; + } + } + } + } + // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + clock_t starttimebubble = clock(); + APECSS_FLOAT tSim = 0.0; // Simulation time of the coupled system + + for (register int i = 0; i < nBubbles; i++) + { + Bubbles[i]->tStart = tSim; + Bubbles[i]->tEnd = (APECSS_FLOAT) tEnd; + Bubbles[i]->dt = APECSS_MIN(1.0e-11, dt_interbubble); // Initial time-step + } + + /* Initialize the bubble based on the selected options */ + for (register int i = 0; i < nBubbles; i++) apecss_bubble_initialize(Bubbles[i]); + + /* Initialize */ + for (register int i = 0; i < nBubbles; i++) apecss_bubble_solver_initialize(Bubbles[i]); + + // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + // Files to retrieve all valuable information for cavitation onset test case + FILE *file_loc; + file_loc = fopen("bubble_loc.txt", "w"); + fprintf(file_loc, "#number x(m) y(m) z(m)\n"); + for (register int i = 0; i < nBubbles; i++) + { + fprintf(file_loc, "%d %e %e %e\n", i, Bubbles[i]->Interaction->location[0], Bubbles[i]->Interaction->location[1], Bubbles[i]->Interaction->location[2]); + } + fclose(file_loc); + + FILE *file_tension; + file_tension = fopen("tension_results.txt", "w"); + fprintf(file_tension, "%d Bubbles p0(pa) %e p1(Pa) %e\n", nBubbles, Liquid->pref, pa); + fprintf(file_tension, "Initial_radii(m)"); + for (register int i = 0; i < nBubbles; i++) fprintf(file_tension, " %e", Bubbles[i]->R0); + fprintf(file_tension, "\n"); + fprintf(file_tension, "#Time(s) R(m) Pt(Pa)\n"); + // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + /* Solve the bubble dynamics */ + while (tSim < (APECSS_FLOAT) tEnd) // Interaction loop, corresponding to the time-intervals at which interactions are considered + { + APECSS_FLOAT dtSim = APECSS_MIN(dt_interbubble, (APECSS_FLOAT) tEnd - tSim); + tSim += dtSim; + + for (register int i = 0; i < nBubbles; i++) apecss_bubble_solver_run(tSim, Bubbles[i]); + + // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + // Retrieve data + fprintf(file_tension, "%e", tSim); + for (register int i = 0; i < nBubbles; i++) + { + fprintf(file_tension, " %e", Bubbles[i]->R); + } + for (register int i = 0; i < nBubbles; i++) + { + fprintf(file_tension, " %e", Bubbles[i]->get_pressure_infinity(Bubbles[i]->t, Bubbles[i])); + } + fprintf(file_tension, "\n"); + // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + } + + fclose(file_tension); + + /* Finalize the simulation*/ + for (register int i = 0; i < nBubbles; i++) apecss_bubble_solver_finalize(Bubbles[i]); + + char str[APECSS_STRINGLENGTH_SPRINTF]; + for (register int i = 0; i < nBubbles; i++) + { + sprintf(str, "Bubble %i: Solver concluded %i time-steps and %i sub-iterations in %.3f s.", i, Bubbles[i]->dtNumber, Bubbles[i]->nSubIter, + (double) (clock() - starttimebubble) / CLOCKS_PER_SEC); + apecss_writeonscreen(str); + } + + for (register int i = 0; i < nBubbles; i++) apecss_results_rayleighplesset_write(Bubbles[i], APECSS_RESULTS_WRITE); + for (register int i = 0; i < nBubbles; i++) apecss_results_emissionsspace_write(Bubbles[i], APECSS_RESULTS_WRITE); + + /* Make sure all allocated memory is freed */ + for (register int i = 0; i < nBubbles; i++) apecss_bubble_freestruct(Bubbles[i]); + + for (register int i = 0; i < nBubbles; i++) free(Bubbles[i]); + free(Gas); + free(Liquid); + free(Interface); + free(NumericsODE); + free(Excitation); + + return (0); +} + +APECSS_FLOAT interaction_bubble_pressure_infinity(APECSS_FLOAT t, struct APECSS_Bubble *Bubble) +{ + APECSS_FLOAT tau = 1.75e-6; + if (t < tau) + { + return (Bubble->p0 - (Bubble->p0 - Bubble->Excitation->dp) * APECSS_POW2(APECSS_SIN(APECSS_PI * t / tau)) + Bubble->Interaction->dp_neighbor); + } + else + { + return (Bubble->p0 + Bubble->Interaction->dp_neighbor); + } +} + +APECSS_FLOAT interaction_bubble_pressurederivative_infinity(APECSS_FLOAT t, struct APECSS_Bubble *Bubble) +{ + APECSS_FLOAT tau = 1.75e-6; + if (t < tau) + { + APECSS_FLOAT inv_tau = 1 / tau; + return (-2.0 * APECSS_PI * inv_tau * (Bubble->p0 - Bubble->Excitation->dp) * APECSS_COS(APECSS_PI * t * inv_tau) * APECSS_SIN(APECSS_PI * t * inv_tau)); + } + else + { + return (0.0); + } +} \ No newline at end of file diff --git a/examples/sphericalclustertensionwave/QA/build/CMakeLists.txt b/examples/sphericalclustertensionwave/QA/build/CMakeLists.txt new file mode 100644 index 0000000..3c801ed --- /dev/null +++ b/examples/sphericalclustertensionwave/QA/build/CMakeLists.txt @@ -0,0 +1,41 @@ +cmake_minimum_required(VERSION 3.12) + +project (sphericalclustertensionwave_apecss) +set(CMAKE_C_COMPILER mpicc) + +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(sphericalclustertensionwave_apecss ${myfiles}) +target_link_libraries(sphericalclustertensionwave_apecss ${mylibs}) \ No newline at end of file diff --git a/examples/sphericalclustertensionwave/QA/build/compile.sh b/examples/sphericalclustertensionwave/QA/build/compile.sh new file mode 100755 index 0000000..564c7a4 --- /dev/null +++ b/examples/sphericalclustertensionwave/QA/build/compile.sh @@ -0,0 +1,5 @@ +#!/bin/bash +rm sphericalclustertensionwave_apecss +cmake CMakeLists.txt -DCMAKE_BUILD_TYPE=Release +make +rm -r CMakeCache.txt CMakeFiles Makefile cmake_install.cmake \ No newline at end of file diff --git a/examples/sphericalclustertensionwave/QA/execmpi.sh b/examples/sphericalclustertensionwave/QA/execmpi.sh new file mode 100755 index 0000000..c1bbf3c --- /dev/null +++ b/examples/sphericalclustertensionwave/QA/execmpi.sh @@ -0,0 +1,18 @@ +nbubbles=250 +ncores=10 + +cd build +./compile.sh +cd .. +mpiexec -n $ncores ./build/sphericalclustertensionwave_apecss -options run.apecss -amp -3.0e4 -tend 25.0e-06 -cldistrib 0 +python3 gather_results.py + +rm -rf bubble_loc.txt +for ((c=0; c<$ncores; c++)) +do + rm -rf tension_results_$c.txt +done +for ((c=0; c<$nbubbles; c++)) +do + rm -rf Bubble_$c +done \ No newline at end of file diff --git a/examples/sphericalclustertensionwave/QA/gather_results.py b/examples/sphericalclustertensionwave/QA/gather_results.py new file mode 100644 index 0000000..280efa3 --- /dev/null +++ b/examples/sphericalclustertensionwave/QA/gather_results.py @@ -0,0 +1,74 @@ +import os + +# File designed to recover data for plotting results in cavitation onset case + +count_core = 0 +for file in os.listdir() : + if "tension_results_" in file : + count_core += 1 + +lines_tension = [] +for i in range(count_core) : + file_tension = open("tension_results_{}.txt".format(i), "r") + lines = file_tension.readlines() + file_tension.close() + lines_tension.append(lines) + +firstline = lines_tension[0][0].split(" ") +count = int(firstline[0]) +p1 = float(firstline[5]) + +file_name = "mono_{}_{:.4E}.txt".format(count, p1) + +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") + +results_tension = open(os.path.join(results_path, file_name), "w") +# Header +for line in lines_tension[0][:3] : + results_tension.write(line) + +# Combining with the good format the values located in the files associated with the cores used for computation +for t in range(3, len(lines_tension[0])) : + time = lines_tension[0][t].split(" ")[0] + r_list = [] + p_list = [] + + for i in range(len(lines_tension)) : + data = lines_tension[i][t].split(" ") + nbubble_core = int((len(data) - 1) / 2) + r_list_file = data[1:nbubble_core+1] + p_list_file = data[nbubble_core+1:2*nbubble_core+1] + + for r in r_list_file : + r_list.append(r) + + for p in p_list_file : + if "\n" in p : + p = p.split("\n")[0] + p_list.append(p) + + results_tension.write(time) + for rad in r_list : + results_tension.write(" {}".format(rad)) + for pres in p_list : + results_tension.write(" {}".format(pres)) + results_tension.write("\n") + +results_tension.close() + +file_loc = open("bubble_loc.txt", "r") +lines_loc = file_loc.readlines() +file_loc.close() + +file_name_loc = "mono_{}_{:.4E}_loc.txt".format(count, p1) + +file_loc_results = open(os.path.join(results_path, file_name_loc), "w") + +for line in lines_loc : + file_loc_results.write(line) + +file_loc_results.close() \ No newline at end of file diff --git a/examples/sphericalclustertensionwave/QA/run.apecss b/examples/sphericalclustertensionwave/QA/run.apecss new file mode 100644 index 0000000..5a2c1d6 --- /dev/null +++ b/examples/sphericalclustertensionwave/QA/run.apecss @@ -0,0 +1,37 @@ +######################################################### +# # +# APECSS Options File # +# # +######################################################### + +BUBBLE +RPModel KM +Emissions QA 250.0e-6 +PressureAmbient 1.0e5 +END + +GAS +EoS IG +ReferenceDensity 1.2 +PolytropicExponent 1.4 +END + +LIQUID +ReferencePressure 1.0e5 +ReferenceDensity 1000.0 +ReferenceSoundSpeed 1500.0 +Viscosity 1.002e-3 +END + +INTERFACE +SurfaceTensionCoeff 0.0728 +END + +RESULTS +Bubble +END + +ODESOLVER +tolerance 1.0e-10 +MinTimeStep 1.0e-10 +END diff --git a/examples/sphericalclustertensionwave/QA/src/sphericalclustertensionwave_apecss.c b/examples/sphericalclustertensionwave/QA/src/sphericalclustertensionwave_apecss.c new file mode 100644 index 0000000..bf1d0df --- /dev/null +++ b/examples/sphericalclustertensionwave/QA/src/sphericalclustertensionwave_apecss.c @@ -0,0 +1,595 @@ +// This source file is part of APECSS, an open-source software toolbox +// for the computation of pressure-driven bubble dynamics and acoustic +// emissions in spherical symmetry. +// +// Copyright (C) 2022-2024 The APECSS Developers +// +// The APECSS Developers are listed in the README.md file available in +// the GitHub repository at https://github.com/polycfd/apecss. +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at https://mozilla.org/MPL/2.0/. + +// ------------------------------------------------------------------- +// APECSS standalone example of acoustically-interacting microbubbles +// in a spherical cluster with the goal to study cavitation onset +// ------------------------------------------------------------------- + +#include +#include +#include "apecss.h" + +APECSS_FLOAT rand_range(double min, double max) +{ + double random = (drand48()); + double range = (max - min) * random; + double number = min + range; + + return (APECSS_FLOAT) number; +} + +struct APECSS_Parallel_Cluster +{ + int rank, size; + int nBubbles_local, nBubbles_global; + + int *bubblerank; // max id of bubbles for each rank + APECSS_FLOAT *bubbleglobal_R, *bubbleglobal_x, *bubbleglobal_y, *bubbleglobal_z; // Radius and x,y,z location of all bubbles + APECSS_FLOAT *sumGU_rank; // Contributions from local bubbles to all bubbles +}; + +// Declaration of additional case-dependent functions +APECSS_FLOAT parallel_interactions_bubble_pressure_infinity(APECSS_FLOAT t, struct APECSS_Bubble *Bubble); +APECSS_FLOAT parallel_interactions_bubble_pressurederivative_infinity(APECSS_FLOAT t, struct APECSS_Bubble *Bubble); + +int parallel_interactions_quasi_acoustic(struct APECSS_Bubble *Bubble[], struct APECSS_Parallel_Cluster *RankInfo); +int parallel_interactions_proper_cutoffdistance(struct APECSS_Bubble *Bubbles[], struct APECSS_Parallel_Cluster *RankInfo); + +int main(int argc, char **args) +{ + /* Initialize MPI */ + int mpi_rank, mpi_size; + MPI_Init(&argc, &args); + MPI_Comm_size(MPI_COMM_WORLD, &mpi_size); + MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank); + + char OptionsDir[APECSS_STRINGLENGTH]; + + // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + // Set the case-dependent simulation parameters + const int nBubbles = 250; // Number of bubbles + APECSS_FLOAT bubble_bubble_dist = 20.0e-6; // Bubble-bubble minimal distance + APECSS_FLOAT cluster_radius = 232e-6; // Spherical cluster radius + + // Interbubble time-step, defining the frequency with which the neighbor influence is updated + APECSS_FLOAT dt_interbubble = 1.0e-9; + + // Initialize the simulation parameters given by the execution command + double tEnd = 0.0; + double fa = 0.0; + double pa = 0.0; + // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + apecss_infoscreen(); + + /* Read commandline options */ + sprintf(OptionsDir, "./run.apecss"); // This is the default + int j = 1; // First argument is the call to the executable + while (j < argc) + { + if (strcmp("-options", args[j]) == 0) + { + sprintf(OptionsDir, "%s", args[j + 1]); + j += 2; + } + else if (strcmp("-tend", args[j]) == 0) + { + sscanf(args[j + 1], "%le", &tEnd); + j += 2; + } + else if (strcmp("-freq", args[j]) == 0) + { + sscanf(args[j + 1], "%le", &fa); + j += 2; + } + else if (strcmp("-amp", args[j]) == 0) + { + sscanf(args[j + 1], "%le", &pa); + j += 2; + } + else if (strcmp("-h", args[j]) == 0) + { + apecss_helpscreen(); + } + else + { + char str[APECSS_STRINGLENGTH_SPRINTF]; + sprintf(str, "Unknown command line options: %s", args[j]); + apecss_erroronscreen(1, str); + ++j; + } + } + + /* Allocate structure for parallel data */ + struct APECSS_Parallel_Cluster *RankInfo = (struct APECSS_Parallel_Cluster *) malloc(sizeof(struct APECSS_Parallel_Cluster)); + RankInfo->rank = mpi_rank; + RankInfo->size = mpi_size; + + /* Determine the number of bubbles per rank */ + int max_per_rank = ceil((double) nBubbles / (double) mpi_size); + RankInfo->nBubbles_local = APECSS_MAX(0, APECSS_MIN(max_per_rank, nBubbles - mpi_rank * max_per_rank)); + + /* Share the parallel distribution of bubbles with all ranks */ + RankInfo->bubblerank = malloc((RankInfo->size + 1) * sizeof(int)); + + RankInfo->nBubbles_global = 0; + RankInfo->bubblerank[0] = 0; + for (int r = 0; r < RankInfo->size; r++) + { + int temp = RankInfo->nBubbles_local; + MPI_Bcast(&temp, 1, MPI_INT, r, MPI_COMM_WORLD); + RankInfo->bubblerank[r + 1] = RankInfo->bubblerank[r] + temp; + RankInfo->nBubbles_global += temp; + } + + /* Allocate and initialize Bubble structure */ + struct APECSS_Bubble *Bubbles[RankInfo->nBubbles_local]; + for (register int i = 0; i < RankInfo->nBubbles_local; i++) Bubbles[i] = (struct APECSS_Bubble *) malloc(sizeof(struct APECSS_Bubble)); + for (register int i = 0; i < RankInfo->nBubbles_local; i++) apecss_bubble_initializestruct(Bubbles[i]); + + /* Set default options and read the options for the bubble */ + for (register int i = 0; i < RankInfo->nBubbles_local; i++) apecss_bubble_setdefaultoptions(Bubbles[i]); + for (register int i = 0; i < RankInfo->nBubbles_local; i++) apecss_bubble_readoptions(Bubbles[i], OptionsDir); + + /* 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)); + struct APECSS_Interface *Interface = (struct APECSS_Interface *) malloc(sizeof(struct APECSS_Interface)); + struct APECSS_NumericsODE *NumericsODE = (struct APECSS_NumericsODE *) malloc(sizeof(struct APECSS_NumericsODE)); + + /* Set the default options for the fluid properties and solver parameters */ + apecss_gas_setdefaultoptions(Gas); + apecss_liquid_setdefaultoptions(Liquid); + apecss_interface_setdefaultoptions(Interface); + apecss_odesolver_setdefaultoptions(NumericsODE); + + /* Read the options file for the fluid properties and solver parameters */ + apecss_gas_readoptions(Gas, OptionsDir); + apecss_liquid_readoptions(Liquid, OptionsDir); + apecss_interface_readoptions(Interface, OptionsDir); + apecss_odesolver_readoptions(NumericsODE, OptionsDir); + + /* Associate the bubble with the relevant fluid properties and solver parameters */ + for (register int i = 0; i < RankInfo->nBubbles_local; i++) + { + Bubbles[i]->Gas = Gas; + Bubbles[i]->Liquid = Liquid; + Bubbles[i]->Interface = Interface; + Bubbles[i]->NumericsODE = NumericsODE; + } + + /* Allocate and set the excitation parameters */ + struct APECSS_Excitation *Excitation = (struct APECSS_Excitation *) malloc(sizeof(struct APECSS_Excitation)); + Excitation->type = APECSS_EXCITATION_SIN; + Excitation->f = (APECSS_FLOAT) fa; + Excitation->dp = (APECSS_FLOAT) pa; + for (register int i = 0; i < RankInfo->nBubbles_local; i++) Bubbles[i]->Excitation = Excitation; + + // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + // Create individual folders for the results of each bubble + for (register int i = 0; i < RankInfo->nBubbles_local; i++) + { + if (Bubbles[i]->Results != NULL) + { + sprintf(Bubbles[i]->Results->dir, "./Bubble_%i/", RankInfo->bubblerank[RankInfo->rank] + i); + struct stat st = {0}; + if (stat(Bubbles[i]->Results->dir, &st) == -1) mkdir(Bubbles[i]->Results->dir, 0700); + } + } + // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + /* Process all options */ + apecss_gas_processoptions(Gas); + apecss_liquid_processoptions(Liquid); + apecss_interface_processoptions(Interface); + apecss_odesolver_processoptions(NumericsODE); + for (register int i = 0; i < RankInfo->nBubbles_local; i++) apecss_bubble_processoptions(Bubbles[i]); + + // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + // Use the revised pressure at infinity, including neighbor contributions + for (register int i = 0; i < RankInfo->nBubbles_local; i++) Bubbles[i]->get_pressure_infinity = parallel_interactions_bubble_pressure_infinity; + for (register int i = 0; i < RankInfo->nBubbles_local; i++) + Bubbles[i]->get_pressurederivative_infinity = parallel_interactions_bubble_pressurederivative_infinity; + + // Allocate interaction structure + for (register int i = 0; i < RankInfo->nBubbles_local; i++) + { + Bubbles[i]->Interaction = (struct APECSS_Interaction *) malloc(sizeof(struct APECSS_Interaction)); + 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; + } + + // Update interaction structure + for (register int i = 0; i < RankInfo->nBubbles_local; i++) Bubbles[i]->Interaction->nBubbles = nBubbles; + + // Define the size of each bubble + for (register int i = 0; i < RankInfo->nBubbles_local; i++) + { + Bubbles[i]->R0 = 2.0e-6; + Bubbles[i]->R = Bubbles[i]->R0; + } + + // Define center location for each bubble + APECSS_FLOAT Bubble_Center[nBubbles][3]; + + for (register int i = 0; i < nBubbles; i++) + { + if (i == 0) + { + Bubble_Center[i][0] = 0.0; + Bubble_Center[i][1] = 0.0; + Bubble_Center[i][2] = 0.0; + } + else + { + APECSS_FLOAT x = rand_range((double) -cluster_radius, (double) cluster_radius); + APECSS_FLOAT y = rand_range((double) -cluster_radius, (double) cluster_radius); + APECSS_FLOAT z = rand_range((double) -cluster_radius, (double) cluster_radius); + + APECSS_FLOAT radius = APECSS_SQRT(APECSS_POW2(x) + APECSS_POW2(y) + APECSS_POW2(z)); + + while (radius > cluster_radius) + { + x = rand_range((double) -cluster_radius, (double) cluster_radius); + y = rand_range((double) -cluster_radius, (double) cluster_radius); + z = rand_range((double) -cluster_radius, (double) cluster_radius); + + radius = APECSS_SQRT(APECSS_POW2(x) + APECSS_POW2(y) + APECSS_POW2(z)); + } + + Bubble_Center[i][0] = x; + Bubble_Center[i][1] = y; + Bubble_Center[i][2] = z; + + for (register int k = 0; k < i; k++) + { + APECSS_FLOAT bubbledist = APECSS_SQRT(APECSS_POW2(Bubble_Center[i][0] - Bubble_Center[k][0]) + APECSS_POW2(Bubble_Center[i][1] - Bubble_Center[k][1]) + + APECSS_POW2(Bubble_Center[i][2] - Bubble_Center[k][2])); + if (bubbledist < bubble_bubble_dist) + { + i--; + break; + } + } + } + } + + for (register int n = 0; n < RankInfo->nBubbles_local; n++) + { + Bubbles[n]->Interaction->location[0] = Bubble_Center[RankInfo->bubblerank[RankInfo->rank] + n][0]; + Bubbles[n]->Interaction->location[1] = Bubble_Center[RankInfo->bubblerank[RankInfo->rank] + n][1]; + Bubbles[n]->Interaction->location[2] = Bubble_Center[RankInfo->bubblerank[RankInfo->rank] + n][2]; + } + + // Share the location of each bubble with all ranks + + RankInfo->bubbleglobal_R = malloc(RankInfo->nBubbles_global * sizeof(APECSS_FLOAT)); + RankInfo->bubbleglobal_x = malloc(RankInfo->nBubbles_global * sizeof(APECSS_FLOAT)); + RankInfo->bubbleglobal_y = malloc(RankInfo->nBubbles_global * sizeof(APECSS_FLOAT)); + RankInfo->bubbleglobal_z = malloc(RankInfo->nBubbles_global * sizeof(APECSS_FLOAT)); + + for (int r = 0; r < RankInfo->size; r++) + { + int temp = RankInfo->nBubbles_local; + MPI_Bcast(&temp, 1, MPI_INT, r, MPI_COMM_WORLD); + + APECSS_FLOAT temp_array[4]; + + for (int i = 0; i < temp; i++) + { + if (r == RankInfo->rank) + { + temp_array[0] = Bubbles[i]->Interaction->location[0]; + temp_array[1] = Bubbles[i]->Interaction->location[1]; + temp_array[2] = Bubbles[i]->Interaction->location[2]; + temp_array[3] = Bubbles[i]->R; + } + + MPI_Bcast(temp_array, 4, MPI_DOUBLE, r, MPI_COMM_WORLD); + + RankInfo->bubbleglobal_x[RankInfo->bubblerank[r] + i] = temp_array[0]; + RankInfo->bubbleglobal_y[RankInfo->bubblerank[r] + i] = temp_array[1]; + RankInfo->bubbleglobal_z[RankInfo->bubblerank[r] + i] = temp_array[2]; + RankInfo->bubbleglobal_R[RankInfo->bubblerank[r] + i] = temp_array[3]; + } + } + + // Update cut off distance for each bubble + parallel_interactions_proper_cutoffdistance(Bubbles, RankInfo); + // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + clock_t starttimebubble = clock(); + APECSS_FLOAT tSim = 0.0; // Simulation time of the coupled system + + for (register int i = 0; i < RankInfo->nBubbles_local; i++) + { + Bubbles[i]->tStart = tSim; + Bubbles[i]->tEnd = (APECSS_FLOAT) tEnd; + Bubbles[i]->dt = APECSS_MIN(1.0e-11, dt_interbubble); // Initial time-step + } + + /* Initialize the bubble based on the selected options */ + for (register int i = 0; i < RankInfo->nBubbles_local; i++) apecss_bubble_initialize(Bubbles[i]); + + /* Initialize */ + for (register int i = 0; i < RankInfo->nBubbles_local; i++) apecss_bubble_solver_initialize(Bubbles[i]); + + // Allocate the pressure contribution array + RankInfo->sumGU_rank = malloc(2 * RankInfo->nBubbles_global * sizeof(APECSS_FLOAT)); + + // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + // Files to retrieve all valuable information for cavitation onset test case + // Bubbles locations + FILE *file_loc = fopen("bubble_loc.txt", "w"); + fprintf(file_loc, "#number x(m) y(m) z(m)\n"); + for (register int i = 0; i < RankInfo->nBubbles_global; i++) + { + fprintf(file_loc, "%d %e %e %e\n", i, RankInfo->bubbleglobal_x[i], RankInfo->bubbleglobal_y[i], RankInfo->bubbleglobal_z[i]); + } + fclose(file_loc); + + // Radius and pressure evolution for each bubble during computation + FILE *file_tension; + char file_name[APECSS_STRINGLENGTH_SPRINTF_LONG]; + sprintf(file_name, "tension_results_%d.txt", RankInfo->rank); + file_tension = fopen(file_name, "w"); + + fprintf(file_tension, "%d Bubbles p0(pa) %e p1(Pa) %e\n", nBubbles, Liquid->pref, pa); + fprintf(file_tension, "Initial_radii(m)"); + for (register int i = 0; i < RankInfo->nBubbles_global; i++) + { + fprintf(file_tension, " %e", RankInfo->bubbleglobal_R[i]); + } + fprintf(file_tension, "\n"); + fprintf(file_tension, "#Time(s) R(m) Pt(Pa)\n"); + // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + /* Solve the bubble dynamics */ + while (tSim < (APECSS_FLOAT) tEnd) // Interaction loop, corresponding to the time-intervals at which interactions are considered + { + APECSS_FLOAT dtSim = APECSS_MIN(dt_interbubble, (APECSS_FLOAT) tEnd - tSim); + tSim += dtSim; + + for (register int i = 0; i < RankInfo->nBubbles_local; i++) apecss_bubble_solver_run(tSim, Bubbles[i]); + + // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + // Update the contribution of the neighbor bubble + parallel_interactions_quasi_acoustic(Bubbles, RankInfo); + // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + // Retrieve data + fprintf(file_tension, "%e", tSim); + for (register int i = 0; i < RankInfo->nBubbles_local; i++) + { + fprintf(file_tension, " %e", Bubbles[i]->R); + } + for (register int i = 0; i < RankInfo->nBubbles_local; i++) + { + fprintf(file_tension, " %e", Bubbles[i]->get_pressure_infinity(Bubbles[i]->t, Bubbles[i])); + } + fprintf(file_tension, "\n"); + // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + for (register int i = 0; i < RankInfo->nBubbles_local; i++) + { + Bubbles[i]->Interaction->last_t_2 = Bubbles[i]->Interaction->last_t_1; + Bubbles[i]->Interaction->last_p_2 = Bubbles[i]->Interaction->last_p_1; + + Bubbles[i]->Interaction->last_t_1 = tSim; + Bubbles[i]->Interaction->last_p_1 = Bubbles[i]->Interaction->dp_neighbor; + } + } + + fclose(file_tension); + + /* Finalize the simulation*/ + for (register int i = 0; i < RankInfo->nBubbles_local; i++) apecss_bubble_solver_finalize(Bubbles[i]); + + char str[APECSS_STRINGLENGTH_SPRINTF]; + for (register int i = 0; i < RankInfo->nBubbles_local; i++) + { + sprintf(str, "Bubble %i: Solver concluded %i time-steps and %i sub-iterations in %.3f s.", RankInfo->bubblerank[RankInfo->rank] + i, Bubbles[i]->dtNumber, + Bubbles[i]->nSubIter, (double) (clock() - starttimebubble) / CLOCKS_PER_SEC); + apecss_writeonscreen(str); + } + + for (register int i = 0; i < RankInfo->nBubbles_local; i++) apecss_results_rayleighplesset_write(Bubbles[i], APECSS_RESULTS_WRITE); + for (register int i = 0; i < RankInfo->nBubbles_local; i++) apecss_results_emissionsspace_write(Bubbles[i], APECSS_RESULTS_WRITE); + + /* Make sure all allocated memory is freed */ + for (register int i = 0; i < RankInfo->nBubbles_local; i++) apecss_bubble_freestruct(Bubbles[i]); + + for (register int i = 0; i < RankInfo->nBubbles_local; i++) free(Bubbles[i]); + free(Gas); + free(Liquid); + free(Interface); + free(NumericsODE); + free(Excitation); + + free(RankInfo->bubblerank); + RankInfo->bubblerank = NULL; + free(RankInfo->bubbleglobal_R); + RankInfo->bubbleglobal_R = NULL; + free(RankInfo->bubbleglobal_x); + RankInfo->bubbleglobal_x = NULL; + free(RankInfo->bubbleglobal_y); + RankInfo->bubbleglobal_y = NULL; + free(RankInfo->bubbleglobal_z); + RankInfo->bubbleglobal_z = NULL; + free(RankInfo->sumGU_rank); + RankInfo->sumGU_rank = NULL; + free(RankInfo); + + MPI_Finalize(); + + return (0); +} + +APECSS_FLOAT parallel_interactions_bubble_pressure_infinity(APECSS_FLOAT t, struct APECSS_Bubble *Bubble) +{ + APECSS_FLOAT tau = 1.75e-6; + if (t < tau) + { + return (Bubble->p0 - (Bubble->p0 - Bubble->Excitation->dp) * APECSS_POW2(APECSS_SIN(APECSS_PI * t / tau)) + Bubble->Interaction->dp_neighbor); + } + else + { + return (Bubble->p0 + Bubble->Interaction->dp_neighbor); + } +} + +APECSS_FLOAT parallel_interactions_bubble_pressurederivative_infinity(APECSS_FLOAT t, struct APECSS_Bubble *Bubble) +{ + // Approximate numerical computation of p_infinity derivative + APECSS_FLOAT tau = 1.75e-6; + APECSS_FLOAT derivative = 0.0; + if (t < tau) + { + APECSS_FLOAT inv_tau = 1 / tau; + derivative = -2.0 * APECSS_PI * inv_tau * (Bubble->p0 - Bubble->Excitation->dp) * APECSS_COS(APECSS_PI * t * inv_tau) * APECSS_SIN(APECSS_PI * t * inv_tau); + } + APECSS_FLOAT delta_t = Bubble->Interaction->last_t_1 - Bubble->Interaction->last_t_2; + if (delta_t > Bubble->dt) + { + APECSS_FLOAT inv_delta_t = 1 / delta_t; + return (derivative + ((Bubble->Interaction->last_p_1 - Bubble->Interaction->last_p_2) * inv_delta_t)); + } + else + { + return (derivative); + } +} + +int parallel_interactions_quasi_acoustic(struct APECSS_Bubble *Bubbles[], struct APECSS_Parallel_Cluster *RankInfo) +{ + // All bubbles are supposed to be in the same liquid + struct APECSS_Liquid *Liquid = Bubbles[0]->Liquid; + + // Update bubble radii info + APECSS_FLOAT *tempR = malloc(RankInfo->nBubbles_global * sizeof(APECSS_FLOAT)); + for (register int j = 0; j < RankInfo->nBubbles_global; j++) tempR[j] = 0.0; + for (register int i = 0; i < RankInfo->nBubbles_local; i++) tempR[RankInfo->bubblerank[RankInfo->rank] + i] = Bubbles[i]->R; + MPI_Allreduce(tempR, RankInfo->bubbleglobal_R, RankInfo->nBubbles_global, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + free(tempR); + + // Reset pressure contributions of the neighours + for (register int i = 0; i < RankInfo->nBubbles_local; i++) Bubbles[i]->Interaction->dp_neighbor = 0.0; + + // Locally compute the contribution to each bubble + for (register int j = 0; j < RankInfo->nBubbles_global; j++) + { + APECSS_FLOAT x_j = RankInfo->bubbleglobal_x[j]; + APECSS_FLOAT y_j = RankInfo->bubbleglobal_y[j]; + APECSS_FLOAT z_j = RankInfo->bubbleglobal_z[j]; + + RankInfo->sumGU_rank[j] = 0.0; + RankInfo->sumGU_rank[j + RankInfo->nBubbles_global] = 0.0; + + double R_j = RankInfo->bubbleglobal_R[j]; + + for (register int i = 0; i < RankInfo->nBubbles_local; i++) + { + if (j != RankInfo->bubblerank[RankInfo->rank] + i) + { + APECSS_FLOAT sumU_bubble = 0.0; + APECSS_FLOAT sumG_bubble = 0.0; + int nodecount_bubble = 0; + + APECSS_FLOAT x_i = Bubbles[i]->Interaction->location[0]; + APECSS_FLOAT y_i = Bubbles[i]->Interaction->location[1]; + APECSS_FLOAT z_i = Bubbles[i]->Interaction->location[2]; + + APECSS_FLOAT interbubble_dist = APECSS_SQRT(APECSS_POW2(x_i - x_j) + APECSS_POW2(y_i - y_j) + APECSS_POW2(z_i - z_j)); + + // The loop over the emission nodes begins with the most far away from the emitting bubble + struct APECSS_EmissionNode *Current = Bubbles[i]->Emissions->LastNode; + while (Current != NULL) + { + if (Current->r < interbubble_dist - R_j) + { + // The following emission nodes have not yet reached the bubble of interest + Current = NULL; + } + else if (Current->r > interbubble_dist + R_j) + { + // The bubble of interest is still not reached + Current = Current->backward; + } + else + { + // The current emission node is located inside the bubble of interest + APECSS_FLOAT gr = Current->g / Current->r; + sumU_bubble += (Current->f / APECSS_POW2(Current->r)) + gr / Liquid->cref; + sumG_bubble += gr; + nodecount_bubble++; + Current = Current->backward; + } + } + + if (nodecount_bubble) + { + APECSS_FLOAT inv_nodecount = 1.0 / (APECSS_FLOAT) nodecount_bubble; + RankInfo->sumGU_rank[j] += sumG_bubble * inv_nodecount; + RankInfo->sumGU_rank[j + RankInfo->nBubbles_global] += sumU_bubble * inv_nodecount; + } + } + } + } + + APECSS_FLOAT *sumGU_all = malloc(2 * RankInfo->nBubbles_global * sizeof(APECSS_FLOAT)); + MPI_Allreduce(RankInfo->sumGU_rank, sumGU_all, 2 * RankInfo->nBubbles_global, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + + // Compute the total pressure contributions of the neighbours for all local bubbles + for (register int i = 0; i < RankInfo->nBubbles_local; i++) + { + Bubbles[i]->Interaction->dp_neighbor = + Liquid->rhoref * (sumGU_all[RankInfo->bubblerank[RankInfo->rank] + i] - + 0.5 * APECSS_POW2(sumGU_all[RankInfo->nBubbles_global + RankInfo->bubblerank[RankInfo->rank] + i])); + } + + free(sumGU_all); + + return (0); +} + +int parallel_interactions_proper_cutoffdistance(struct APECSS_Bubble *Bubbles[], struct APECSS_Parallel_Cluster *RankInfo) +{ + for (register int i = 0; i < RankInfo->nBubbles_local; i++) + { + APECSS_FLOAT x_i = Bubbles[i]->Interaction->location[0]; + APECSS_FLOAT y_i = Bubbles[i]->Interaction->location[1]; + APECSS_FLOAT z_i = Bubbles[i]->Interaction->location[2]; + APECSS_FLOAT dist = 0.0; + + for (register int j = 0; j < RankInfo->nBubbles_global; j++) + { + APECSS_FLOAT x_j = RankInfo->bubbleglobal_x[j]; + APECSS_FLOAT y_j = RankInfo->bubbleglobal_y[j]; + APECSS_FLOAT z_j = RankInfo->bubbleglobal_z[j]; + + APECSS_FLOAT dist_ij = APECSS_SQRT(APECSS_POW2(x_i - x_j) + APECSS_POW2(y_i - y_j) + APECSS_POW2(z_i - z_j)); + + if (dist_ij > dist) + { + dist = dist_ij; + } + } + if (Bubbles[i]->Emissions != NULL) Bubbles[i]->Emissions->CutOffDistance = 1.25 * dist; + } + return (0); +} \ No newline at end of file diff --git a/examples/sphericalclustertensionwave/README.md b/examples/sphericalclustertensionwave/README.md new file mode 100644 index 0000000..127679e --- /dev/null +++ b/examples/sphericalclustertensionwave/README.md @@ -0,0 +1,3 @@ +This example builds a standalone APECSS code for a monodispersed cluster of 250 interacting microbubbles. The cluster has a spherical shape with radius $R_{\mathrm{C}} = 232 \ \mu\mathrm{m}$, with the bubbles randomly distributed inside. The bubbles initial radius is $R_{0} = 2 \mu\mathrm{m}$. The incident wave is a single tension pulse defined by $p_{\mathrm{a}}(t) = -(p_{0} - p_{\mathrm{ng}})\sin^{2}\left(\pi t / \tau \right)$, with $\tau = 1.75 \ \mu\mathrm{s}$ and $p_{\mathrm{ng}} = -3.0 \times 10^{4} \ \mathrm{Pa}$. The computations are done with several interaction models : *no interactions* ("NI"), *incompressible interactions* ("IC") or *quasi-acoustic interactions* ("QA"). Each interaction model has its own repository for gathering data. The [run_sphericalclustertensionwave.sh](run_sphericalclustertensionwave.sh) file provides the whole execution command to run all computations, plot the results and clean the folders. + +The quasi-acoustic model computation is done using parallelization, and takes around **30 minutes** to finish when running on 12 cores of an Intel 13th-Gen i7-13700K x 24 processor. The number of cores used can be changed in the header section of [run_sphericalclustertensionwave.sh](run_sphericalclustertensionwave.sh). \ No newline at end of file diff --git a/examples/sphericalclustertensionwave/plot_results.py b/examples/sphericalclustertensionwave/plot_results.py new file mode 100644 index 0000000..94d8f3f --- /dev/null +++ b/examples/sphericalclustertensionwave/plot_results.py @@ -0,0 +1,353 @@ +import os +import numpy as np +import matplotlib.pyplot as plt +from math import sqrt + +fontsize = 15 + +plt.rcParams['font.family']='serif' +plt.rcParams['font.serif']=['Times New Roman'] + plt.rcParams['font.serif'] +plt.rcParams['mathtext.fontset']='stix' +plt.rcParams['font.size']=fontsize + +cm = 1/2.54 + +# File designed to plot the results for the cavitation onset test case with a spherical cluster + +######### Step 1 : Retrieving data ################################################################################################################## +interaction_types = ["NI", "IC", "QA"] +cluster_types = ["mono"] + +dic_bubbles = {} +dic_bubbles_loc = {} + +for inttype in interaction_types : + if inttype not in list(dic_bubbles.keys()) : + dic_bubbles[inttype] = {} + if inttype not in list(dic_bubbles_loc.keys()) : + dic_bubbles_loc[inttype] = {} + for cltype in cluster_types : + if cltype not in list(dic_bubbles[inttype].keys()) : + dic_bubbles[inttype][cltype] = {} + if cltype not in list(dic_bubbles_loc[inttype].keys()) : + dic_bubbles_loc[inttype][cltype] = {} + +working_path = os.getcwd() +for inttype in interaction_types : + interaction_path = os.path.join(working_path, inttype) + interaction_path = os.path.join(interaction_path, "results") + for file in os.listdir(interaction_path) : + if "_loc.txt" in file : + file_name = file.split("_loc.txt")[0] + + file_loc_path = os.path.join(interaction_path, file) + file_loc = open(file_loc_path, "r") + lines_loc = file_loc.readlines() + file_loc.close() + + file_path = os.path.join(interaction_path, file_name + ".txt") + file_res = open(file_path, "r") + lines_res = file_res.readlines() + file_res.close() + + firstline = lines_res[0].split(" ") + count = int(firstline[0]) + p1 = float(firstline[5]) + + dic_res = dic_bubbles[inttype]["mono"] + dic_loc = dic_bubbles_loc[inttype]["mono"] + + if count not in list(dic_res.keys()) : + dic_res[count] = {} + if count not in list(dic_loc.keys()) : + dic_loc[count] = {} + + if p1 not in list(dic_res[count].keys()) : + dic_res[count][p1] = [] + if p1 not in list(dic_loc[count].keys()) : + dic_loc[count][p1] = [] + + dic_res[count][p1].append([]) + for i in range(count) : + init_radius = float(lines_res[1].split(" ")[1 + i]) + x = float(lines_loc[1 + i].split(" ")[1]) + y = float(lines_loc[1 + i].split(" ")[2]) + z = float(lines_loc[1 + i].split(" ")[3]) + + dic_loc[count][p1].append([x, y, z]) + dic_res[count][p1].append([init_radius, [], []]) + + for line in lines_res[3:] : + data = line.split(" ") + t = float(data[0]) + dic_res[count][p1][0].append(t) + + for i in range(count) : + r = float(data[1 + i]) + p = float(data[1 + count + i]) + dic_res[count][p1][i+1][1].append(r) + dic_res[count][p1][i+1][2].append(p) + +######### Step 2 : Define bubble distributions ###################################################################################################### + +# For monodispersed cluster, results are plotted based on bubbles locations +cluster_radius = 232.0e-6 +interval_size = cluster_radius / 10 + +dic_loc_distrib_global = {} + +for inttype in list(dic_bubbles_loc.keys()) : + if inttype not in list(dic_loc_distrib_global.keys()) : + dic_loc_distrib_global[inttype] = {} + + for count in list(dic_bubbles_loc[inttype]["mono"].keys()) : + if count not in list(dic_loc_distrib_global[inttype].keys()) : + dic_loc_distrib_global[inttype][count] = {} + + for p1 in list(dic_bubbles_loc[inttype]["mono"][count].keys()) : + if p1 not in list(dic_loc_distrib_global[inttype][count].keys()) : + dic_loc_distrib_global[inttype][count][p1] = {0.0 : [], 0.25 : [], 0.5 : [], 0.75 : [], 1.0 : []} + + for i in range(count) : + radius_to_center = sqrt(dic_bubbles_loc[inttype]["mono"][count][p1][i][0]**2 + dic_bubbles_loc[inttype]["mono"][count][p1][i][1]**2 + dic_bubbles_loc[inttype]["mono"][count][p1][i][2]**2) + if radius_to_center < 0.0 * cluster_radius + interval_size : + dic_loc_distrib_global[inttype][count][p1][0.0].append(i) + elif 0.25 * cluster_radius - 0.5 * interval_size < radius_to_center < 0.25 * cluster_radius + 0.5 * interval_size : + dic_loc_distrib_global[inttype][count][p1][0.25].append(i) + elif 0.5 * (cluster_radius - interval_size) < radius_to_center < 0.5 * (cluster_radius + interval_size) : + dic_loc_distrib_global[inttype][count][p1][0.5].append(i) + elif 0.75 * cluster_radius - 0.5 * interval_size < radius_to_center < 0.75 * cluster_radius + 0.5 * interval_size : + dic_loc_distrib_global[inttype][count][p1][0.75].append(i) + elif 1.0 * cluster_radius - interval_size < radius_to_center : + dic_loc_distrib_global[inttype][count][p1][1.0].append(i) + +dic_loc_label = {1.0 : "{:.1f}".format(1.0 - interval_size/cluster_radius) + r" $\leq$ $r/R_{C}$ $\leq$ " + "{:.1f}".format(1.0), + 0.75 : "{:.2f}".format(0.75 - 0.5 * interval_size/cluster_radius) + r" $\leq$ $r/R_{C}$ $\leq$ " + "{:.2f}".format(0.75 + 0.5 * interval_size/cluster_radius), + 0.5 : "{:.2f}".format(0.5 - 0.5 * interval_size/cluster_radius) + r" $\leq$ $r/R_{C}$ $\leq$ " + "{:.2f}".format(0.5 + 0.5 * interval_size/cluster_radius), + 0.25 : "{:.2f}".format(0.25 - 0.5 * interval_size/cluster_radius) + r" $\leq$ $r/R_{C}$ $\leq$ " + "{:.2f}".format(0.25 + 0.5 * interval_size/cluster_radius), + 0.0 : "{:.1f}".format(0.0) + r" $\leq$ $r/R_{C}$ $\leq$ " + "{:.1f}".format(interval_size/cluster_radius)} + +######### Step 3 : Plot results ##################################################################################################################### + +#### General parameters ##### +rho_l = 1000.0 +r_ref = 2.0e-6 +p0 = 1.0e5 + +######### Radius and minimum distance between bubbles in both configurations ###################### +p1 = -3.0e4 + +for cluster in cluster_types : + dic_loc = dic_bubbles_loc["NI"][cluster] + + for count in list(dic_loc.keys()) : + min_dist = cluster_radius + loc_list = dic_loc[count][p1] + for i in range(count) : + x_i = loc_list[i][0] + y_i = loc_list[i][1] + z_i = loc_list[i][2] + min_dist_bubble = cluster_radius + + for j in range(count) : + if (i != j) : + x_j = loc_list[j][0] + y_j = loc_list[j][1] + z_j = loc_list[j][2] + dist = sqrt((x_i - x_j)**2 + (y_i - y_j)**2 + (z_i - z_j)**2) + if dist < min_dist_bubble : + min_dist_bubble = dist + + if min_dist_bubble < min_dist : + min_dist = min_dist_bubble + + print("{} cluster with N={} bubbles, minimum distance between bubbles equals {:.2f} microns".format(cluster, count, min_dist*1.0e6)) + +######### Averaged radius evolution ############### +count = 250 +p1 = -3.0e4 +t_i = r_ref * sqrt(rho_l / (p0 - p1)) + +nrow = 1 +ncol = 3 +fig, axs = plt.subplots(nrow,ncol,figsize=(ncol*17.5*cm, nrow*12.5*cm)) +for i in range(ncol) : + axs[i].set_xlabel(r"$ t$ [$\mu$s]", fontsize=27.5) + if i == 0 : axs[i].set_ylabel(r"$ / R_{0}$", fontsize=27.5) + axs[i].set_xlim(xmin=0.0, xmax=15.0) + axs[i].tick_params(axis="both", labelsize=25) + axs[i].grid() + +plt.subplots_adjust(wspace=0.35*cm) + +axs[0].set_title(r"No interactions", fontsize=27.5, y=1.025) +axs[1].set_title(r"Incompressible interactions", fontsize=27.5, y=1.025) +axs[2].set_title(r"Quasi-acoustic interactions", fontsize=27.5, y=1.025) + +dic_color_loc = {0.0 : "blue", 0.25 : "magenta", 0.5 : "red", 0.75 : "green", 1.0 : "black"} +dic_lines_loc = {0.0 : "dotted", 0.25 : "dashed", 0.5 : "dashed", 0.75 : "dashed", 1.0 : "solid"} + +for k in [0.0, 0.5, 1.0] : + # No interactions + index_list = dic_loc_distrib_global["NI"][count][p1][k] + t_list = 1.0e06 * np.array(dic_bubbles["NI"]["mono"][count][p1][0]) + + r_list_0 = np.array(dic_bubbles["NI"]["mono"][count][p1][index_list[0] + 1][1]) + init_r_0 = dic_bubbles["NI"]["mono"][count][p1][index_list[0] + 1][0] + avg_radius = r_list_0 / init_r_0 + + for i in range(1, len(index_list)) : + index = index_list[i] + + r_list = np.array(dic_bubbles["NI"]["mono"][count][p1][index + 1][1]) + init_r = dic_bubbles["NI"]["mono"][count][p1][index + 1][0] + avg_radius = avg_radius + np.array(r_list) / init_r + + avg_radius = (1 / len(index_list)) * avg_radius + + axs[0].plot(t_list, avg_radius, linewidth=3.0, linestyle=dic_lines_loc[k], color=dic_color_loc[k]) + + # Incompressible interactions + index_list = dic_loc_distrib_global["IC"][count][p1][k] + t_list = 1.0e06 * np.array(dic_bubbles["IC"]["mono"][count][p1][0]) + + r_list_0 = np.array(dic_bubbles["IC"]["mono"][count][p1][index_list[0] + 1][1]) + init_r_0 = dic_bubbles["IC"]["mono"][count][p1][index_list[0] + 1][0] + avg_radius = r_list_0 / init_r_0 + + for i in range(1, len(index_list)) : + index = index_list[i] + + r_list = np.array(dic_bubbles["IC"]["mono"][count][p1][index + 1][1]) + init_r = dic_bubbles["IC"]["mono"][count][p1][index + 1][0] + avg_radius = avg_radius + np.array(r_list) / init_r + + avg_radius = (1 / len(index_list)) * avg_radius + + axs[1].plot(t_list, avg_radius, linewidth=3.0, linestyle=dic_lines_loc[k], color=dic_color_loc[k], label=dic_loc_label[k]) + + # Quasi-acoustic interactions + index_list = dic_loc_distrib_global["QA"][count][p1][k] + t_list = 1.0e06 * np.array(dic_bubbles["QA"]["mono"][count][p1][0]) + + r_list_0 = np.array(dic_bubbles["QA"]["mono"][count][p1][index_list[0] + 1][1]) + init_r_0 = dic_bubbles["QA"]["mono"][count][p1][index_list[0] + 1][0] + avg_radius = r_list_0 / init_r_0 + + for i in range(1, len(index_list)) : + index = index_list[i] + + r_list = np.array(dic_bubbles["QA"]["mono"][count][p1][index + 1][1]) + init_r = dic_bubbles["QA"]["mono"][count][p1][index + 1][0] + avg_radius = avg_radius + np.array(r_list) / init_r + + avg_radius = (1 / len(index_list)) * avg_radius + + axs[2].plot(t_list, avg_radius, linewidth=3.0, linestyle=dic_lines_loc[k], color=dic_color_loc[k]) + + +axs[1].legend(bbox_to_anchor=(0.5, 1.175), loc="center", fontsize=26.5, ncol=3, frameon=False) +fig.savefig("sphericalclustertensionwave_radiusevolution.pdf", bbox_inches='tight',pad_inches=0.35) + +######### Averaged pressure infinity evolution ### + +count = 250 +p1 = -3.0e4 +p0 = 1.0e5 +t_i = r_ref * sqrt(rho_l / (p0 - p1)) + +nrow = 1 +ncol = 3 +fig, axs = plt.subplots(nrow,ncol,figsize=(ncol*17.5*cm, nrow*12.5*cm)) +for i in range(ncol) : + axs[i].set_xlabel(r"$ t$ [$\mu$s]", fontsize=27.5) + if i == 0 : axs[i].set_ylabel(r"$ / p_{0}$", fontsize=27.5) + axs[i].set_xlim(xmin=0.0, xmax=15.0) + axs[i].tick_params(axis="both", labelsize=25) + axs[i].grid() + +plt.subplots_adjust(wspace=0.35*cm) + +axs[0].set_title(r"No interactions", fontsize=27.5, y=1.025) +axs[1].set_title(r"Incompressible interactions", fontsize=27.5, y=1.025) +axs[2].set_title(r"Quasi-acoustic interactions", fontsize=27.5, y=1.025) + +dic_color_loc = {0.0 : "blue", 0.25 : "magenta", 0.5 : "red", 0.75 : "green", 1.0 : "black"} +dic_lines_loc = {0.0 : "dotted", 0.25 : "dashed", 0.5 : "dashed", 0.75 : "dashed", 1.0 : "solid"} + +zm_IC = axs[1].inset_axes([0.35, 0.35, 0.60, 0.60]) +zm_IC.grid() +zm_IC.set_xlim(xmin=0.0,xmax=5.0) +zm_IC.set_ylim(ymin=-0.5,ymax=2.0) +zm_IC.set_xticks([0, 2, 4]) +zm_IC.set_yticks([-0.267, 0.0, 1.0, 2.0]) +zm_IC.set_yticklabels([r"$p_{\mathrm{C}}$", 0, 1, 2]) +zm_IC.tick_params(axis="both", labelsize=25) + +zm_QA = axs[2].inset_axes([0.30, 0.4, 0.65, 0.55]) +zm_QA.grid() +zm_QA.set_xlim(xmin=0.0,xmax=5.0) +zm_QA.set_ylim(ymin=-0.5,ymax=2.0) +zm_QA.set_xticks([0, 2, 4]) +zm_QA.set_yticks([-0.267, 0.0, 1.0, 2.0]) +zm_QA.set_yticklabels([r"$p_{\mathrm{C}}$", 0, 1, 2]) +zm_QA.tick_params(axis="both", labelsize=25) + +for k in [0.0, 0.5, 1.0] : + # No interactions + index_list = dic_loc_distrib_global["NI"][count][p1][k] + t_list = 1.0e06 * np.array(dic_bubbles["NI"]["mono"][count][p1][0]) + + p_list_0 = np.array(dic_bubbles["NI"]["mono"][count][p1][index_list[0] + 1][2]) + avg_pressure = p_list_0 / p0 + + for i in range(1, len(index_list)) : + index = index_list[i] + + p_list = np.array(dic_bubbles["NI"]["mono"][count][p1][index + 1][2]) + avg_pressure = avg_pressure + np.array(p_list) / p0 + + avg_pressure = (1 / len(index_list)) * avg_pressure + + axs[0].plot(t_list, avg_pressure, linewidth=3.0, linestyle=dic_lines_loc[k], color=dic_color_loc[k]) + + # Incompressible interactions + index_list = dic_loc_distrib_global["IC"][count][p1][k] + t_list = 1.0e06 * np.array(dic_bubbles["IC"]["mono"][count][p1][0]) + + p_list_0 = np.array(dic_bubbles["IC"]["mono"][count][p1][index_list[0] + 1][2]) + avg_pressure = p_list_0 / p0 + + for i in range(1, len(index_list)) : + index = index_list[i] + + p_list = np.array(dic_bubbles["IC"]["mono"][count][p1][index + 1][2]) + avg_pressure = avg_pressure + np.array(p_list) / p0 + + avg_pressure = (1 / len(index_list)) * avg_pressure + + axs[1].plot(t_list, avg_pressure, linewidth=3.0, linestyle=dic_lines_loc[k], color=dic_color_loc[k], label=dic_loc_label[k]) + + zm_IC.plot(t_list, avg_pressure, linewidth=3.0, linestyle=dic_lines_loc[k], color=dic_color_loc[k]) + + # Quasi-acoustic interactions + index_list = dic_loc_distrib_global["QA"][count][p1][k] + t_list = 1.0e06 * np.array(dic_bubbles["QA"]["mono"][count][p1][0]) + + p_list_0 = np.array(dic_bubbles["QA"]["mono"][count][p1][index_list[0] + 1][2]) + avg_pressure = p_list_0 / p0 + + for i in range(1, len(index_list)) : + index = index_list[i] + + p_list = np.array(dic_bubbles["QA"]["mono"][count][p1][index + 1][2]) + avg_pressure = avg_pressure + np.array(p_list) / p0 + + avg_pressure = (1 / len(index_list)) * avg_pressure + + axs[2].plot(t_list, avg_pressure, linewidth=3.0, linestyle=dic_lines_loc[k], color=dic_color_loc[k]) + + zm_QA.plot(t_list, avg_pressure, linewidth=3.0, linestyle=dic_lines_loc[k], color=dic_color_loc[k]) + +axs[1].legend(bbox_to_anchor=(0.5, 1.175), loc="center", fontsize=26.5, ncol=3, frameon=False) +fig.savefig("sphericalclustertensionwave_pressureevolution.pdf", bbox_inches='tight',pad_inches=0.35) \ No newline at end of file diff --git a/examples/sphericalclustertensionwave/run_sphericalclustertensionwave.sh b/examples/sphericalclustertensionwave/run_sphericalclustertensionwave.sh new file mode 100755 index 0000000..6a017dd --- /dev/null +++ b/examples/sphericalclustertensionwave/run_sphericalclustertensionwave.sh @@ -0,0 +1,89 @@ +######### Parameters ################################################################################################################################ + +ncores=12 +nbubbles=250 + +pressure=-3.0e4 + +######### No interaction computations ############################################################################################################### + +cd NI/build +./compile.sh +cd .. + +./build/sphericalclustertensionwave_apecss -options run.apecss -amp $pressure -tend 15.0e-6 +python3 gather_results.py + +cd .. + +echo "" +echo "No interaction test cases passed" +echo "" + +######### Incompressible computations ############################################################################################################### + +cd IC/build +./compile.sh +cd .. + +./build/sphericalclustertensionwave_apecss -options run.apecss -amp $pressure -tend 15.0e-6 +python3 gather_results.py + +cd .. + +echo "" +echo "Incompressible test cases passed" +echo "" + +######### Quasi acoustic computations ############################################################################################################### + +cd QA/build +./compile.sh +cd .. + +mpiexec -n $ncores ./build/sphericalclustertensionwave_apecss -options run.apecss -amp $pressure -tend 25.0e-06 +python3 gather_results.py + +cd .. + +echo "" +echo "Quasi acoustic test cases passed" +echo "" + +######### Plot results ############################################################################################################################## + +python3 plot_results.py + +######### Cleaning ################################################################################################################################## + +cd NI +rm -rf "tension_results.txt" "bubble_loc.txt" +for ((c=0; c<$nbubbles; c++)) +do + rm -rf Bubble_$c +done +cd .. + +cd IC +rm -rf "tension_results.txt" "bubble_loc.txt" +for ((c=0; c<$nbubbles; c++)) +do + rm -rf Bubble_$c +done +cd .. + +cd QA +rm -rf "bubble_loc.txt" +for ((c=0; c<$ncores; c++)) +do + rm -rf tension_results_$c.txt +done +for ((c=0; c<$nbubbles; c++)) +do + rm -rf Bubble_$c +done +cd .. + +echo "" +echo "Cleaning completed" +echo "" \ No newline at end of file diff --git a/include/apecss.h b/include/apecss.h index 7f5b9d1..d702fd0 100644 --- a/include/apecss.h +++ b/include/apecss.h @@ -311,6 +311,17 @@ struct APECSS_Excitation APECSS_FLOAT dp; // Maximum pressure amplitude [Pa] }; +struct APECSS_Interaction +{ + int nBubbles; // Number of bubbles in the whole cluster considered + APECSS_FLOAT location[3]; // 3D coordinates of the considered bubble + APECSS_FLOAT dp_neighbor; // Pressure induced by interactions with neighboring bubbles + APECSS_FLOAT last_t_1; // Previous timestep when interaction where considered + APECSS_FLOAT last_t_2; // Previous timestep when interaction where considered before last_t_1 + APECSS_FLOAT last_p_1; // Pressure induced by neighboring bubbles or total pressure in the far field at last_t_1 + APECSS_FLOAT last_p_2; // Pressure induced by neighboring bubbles or total pressure in the far field at last_t_2 +}; + struct APECSS_ResultsBubble { int freq; // Frequency with which the results are stored (with respect to the time-step number) @@ -444,6 +455,9 @@ struct APECSS_Bubble // Excitation of the bubble (if applicable) struct APECSS_Excitation *Excitation; + // Interactions with neighboring bubbles (if applicable) + struct APECSS_Interaction *Interaction; + // Pointers to the functions describing the liquid pressure and its derivative at infinity APECSS_FLOAT (*get_pressure_infinity)(APECSS_FLOAT t, struct APECSS_Bubble *Bubble); APECSS_FLOAT (*get_pressurederivative_infinity)(APECSS_FLOAT t, struct APECSS_Bubble *Bubble); @@ -577,6 +591,13 @@ APECSS_FLOAT apecss_interface_pressurederivative_viscous_marmottantexpl(APECSS_F APECSS_FLOAT apecss_interface_pressurederivative_viscous_cleanimpl(APECSS_FLOAT R, struct APECSS_Interface *Interface); APECSS_FLOAT apecss_interface_pressurederivative_viscous_marmottantimpl(APECSS_FLOAT R, struct APECSS_Interface *Interface); +// --------------------- +// interactions.c + +int apecss_interactions_instantaneous(struct APECSS_Bubble *Bubbles[]); +int apecss_interactions_quasi_acoustic(struct APECSS_Bubble *Bubbles[]); +int apecss_interactions_cutoffdistance(struct APECSS_Bubble *Bubbles[]); + // --------------------- // liquid.c diff --git a/src/bubble.c b/src/bubble.c index af564f2..9ba6b7b 100644 --- a/src/bubble.c +++ b/src/bubble.c @@ -39,6 +39,7 @@ int apecss_bubble_initializestruct(struct APECSS_Bubble *Bubble) Bubble->Liquid = NULL; Bubble->Interface = NULL; Bubble->Excitation = NULL; + Bubble->Interaction = NULL; Bubble->Emissions = NULL; Bubble->Results = NULL; @@ -867,6 +868,12 @@ int apecss_bubble_freestruct(struct APECSS_Bubble *Bubble) Bubble->Emissions = NULL; } + if (Bubble->Interaction != NULL) + { + free(Bubble->Interaction); + Bubble->Interaction = NULL; + } + if (Bubble->Results != NULL) { if (Bubble->Results->RayleighPlesset != NULL) diff --git a/src/emissions.c b/src/emissions.c index c5a5cd4..1ded191 100644 --- a/src/emissions.c +++ b/src/emissions.c @@ -147,6 +147,7 @@ int apecss_emissions_prunelist(struct APECSS_Bubble *Bubble) Current->backward->forward = Current->forward; Current->forward->backward = Current->backward; Current = Current->backward; + Bubble->Emissions->nNodes -= 1; free(Obsolete); } else diff --git a/src/interactions.c b/src/interactions.c new file mode 100644 index 0000000..0bf92bf --- /dev/null +++ b/src/interactions.c @@ -0,0 +1,164 @@ +// This source file is part of APECSS, an open-source software toolbox +// for the computation of pressure-driven bubble dynamics and acoustic +// emissions in spherical symmetry. +// +// Copyright (C) 2022-2024 The APECSS Developers +// +// The APECSS Developers are listed in the README.md file available in +// the GitHub repository at https://github.com/polycfd/apecss. +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at https://mozilla.org/MPL/2.0/. + +#include "apecss.h" + +// ------------------------------------------------------------------- +// OPTIONS +// ------------------------------------------------------------------- +// Functions computing different types of interactions between cavita- +// tion bubbles inside a cluster +// ------------------------------------------------------------------- + +int apecss_interactions_instantaneous(struct APECSS_Bubble *Bubbles[]) +{ + // Compute interactions by considering them instantaneous (fully incompressible approach) + int nBubbles = Bubbles[0]->Interaction->nBubbles; + + // All bubbles are supposed to be in the same liquid + struct APECSS_Liquid *Liquid = Bubbles[0]->Liquid; + + for (register int i = 0; i < nBubbles; i++) + { + APECSS_FLOAT x_i = Bubbles[i]->Interaction->location[0]; + APECSS_FLOAT y_i = Bubbles[i]->Interaction->location[1]; + APECSS_FLOAT z_i = Bubbles[i]->Interaction->location[2]; + + Bubbles[i]->Interaction->dp_neighbor = 0.0; + + for (register int j = 0; j < nBubbles; j++) + { + if (j != i) + { + APECSS_FLOAT x_j = Bubbles[j]->Interaction->location[0]; + APECSS_FLOAT y_j = Bubbles[j]->Interaction->location[1]; + APECSS_FLOAT z_j = Bubbles[j]->Interaction->location[2]; + + APECSS_FLOAT interbubble_dist = APECSS_SQRT(APECSS_POW2(x_i - x_j) + APECSS_POW2(y_i - y_j) + APECSS_POW2(z_i - z_j)); + + APECSS_FLOAT dp = Liquid->rhoref * 2.0 * Bubbles[j]->R * APECSS_POW2(Bubbles[j]->U) * (1 / interbubble_dist); + dp += Liquid->rhoref * APECSS_POW2(Bubbles[j]->R) * Bubbles[j]->ode[0](Bubbles[j]->ODEsSol, Bubbles[j]->t, Bubbles[j]) * (1 / interbubble_dist); + dp += -Liquid->rhoref * APECSS_POW4(Bubbles[j]->R) * APECSS_POW2(Bubbles[j]->U) * (1 / (2 * APECSS_POW4(interbubble_dist))); + + Bubbles[i]->Interaction->dp_neighbor += dp; + } + } + } + + return (0); +} + +int apecss_interactions_quasi_acoustic(struct APECSS_Bubble *Bubbles[]) +{ + // Compute interactions using emission nodes in the quasi acoustic assumption + int nBubbles = Bubbles[0]->Interaction->nBubbles; + + // All bubbles are supposed to be in the same liquid + struct APECSS_Liquid *Liquid = Bubbles[0]->Liquid; + + // Reset pressure contributions of the neighbours + for (register int i = 0; i < nBubbles; i++) + { + Bubbles[i]->Interaction->dp_neighbor = 0.0; + } + + for (register int i = 0; i < nBubbles; i++) + { + APECSS_FLOAT sumU = 0.0; + APECSS_FLOAT sumG = 0.0; + + APECSS_FLOAT x_i = Bubbles[i]->Interaction->location[0]; + APECSS_FLOAT y_i = Bubbles[i]->Interaction->location[1]; + APECSS_FLOAT z_i = Bubbles[i]->Interaction->location[2]; + + for (register int j = 0; j < nBubbles; j++) + { + if (j != i) + { + APECSS_FLOAT sumU_bubble = 0.0; + APECSS_FLOAT sumG_bubble = 0.0; + int nodecount_bubble = 0; + + APECSS_FLOAT x_j = Bubbles[j]->Interaction->location[0]; + APECSS_FLOAT y_j = Bubbles[j]->Interaction->location[1]; + APECSS_FLOAT z_j = Bubbles[j]->Interaction->location[2]; + + APECSS_FLOAT interbubble_dist = APECSS_SQRT(APECSS_POW2(x_i - x_j) + APECSS_POW2(y_i - y_j) + APECSS_POW2(z_i - z_j)); + + APECSS_FLOAT r_b = Bubbles[i]->R; + + // The loop over the emission nodes begins with the most far away from the emitting bubble + struct APECSS_EmissionNode *Current = Bubbles[j]->Emissions->LastNode; + while (Current != NULL) + { + if (Current->r < interbubble_dist - r_b) + // The following emission nodes have not yet reached the bubble of interest + Current = NULL; + else if (Current->r > interbubble_dist + r_b) + // The bubble of interest is still not reached + Current = Current->backward; + else + { + // The current emission node is located inside the bubble of interest + APECSS_FLOAT gr = Current->g / Current->r; + sumU_bubble += (Current->f / APECSS_POW2(Current->r)) + (gr) / Liquid->cref; + sumG_bubble += gr; + nodecount_bubble++; + Current = Current->backward; + } + } + + if (nodecount_bubble != 0) + { + APECSS_FLOAT inv_nodecount = 1 / (APECSS_FLOAT) nodecount_bubble; + APECSS_FLOAT sumU_temp = sumU_bubble * inv_nodecount; + APECSS_FLOAT sumG_temp = sumG_bubble * inv_nodecount; + + sumU += (sumU_temp); + sumG += (sumG_temp); + } + } + } + Bubbles[i]->Interaction->dp_neighbor += Liquid->rhoref * (sumG - 0.5 * APECSS_POW2(sumU)); + } + + return (0); +} + +int apecss_interactions_cutoffdistance(struct APECSS_Bubble *Bubbles[]) +{ + // Compute the proper cut off distance for each bubbles in a cluster depending on the location of their neighbors + int nBubbles = Bubbles[0]->Interaction->nBubbles; + + for (register int i = 0; i < nBubbles; i++) + { + APECSS_FLOAT x_i = Bubbles[i]->Interaction->location[0]; + APECSS_FLOAT y_i = Bubbles[i]->Interaction->location[1]; + APECSS_FLOAT z_i = Bubbles[i]->Interaction->location[2]; + APECSS_FLOAT dist = 0.0; + + for (register int j = 0; j < nBubbles; j++) + { + APECSS_FLOAT x_j = Bubbles[j]->Interaction->location[0]; + APECSS_FLOAT y_j = Bubbles[j]->Interaction->location[1]; + APECSS_FLOAT z_j = Bubbles[j]->Interaction->location[2]; + + APECSS_FLOAT dist_ij = APECSS_SQRT(APECSS_POW2(x_i - x_j) + APECSS_POW2(y_i - y_j) + APECSS_POW2(z_i - z_j)); + + if (dist_ij > dist) dist = dist_ij; + } + if (Bubbles[i]->Emissions != NULL) Bubbles[i]->Emissions->CutOffDistance = 1.05 * dist; + } + + return (0); +} \ No newline at end of file