Skip to content

Commit

Permalink
Merge pull request #26 from fabiandenner/main
Browse files Browse the repository at this point in the history 
Invariant f; prune emissions; updated examples.
  • Loading branch information
polycfd authored Apr 16, 2024
2 parents 2986803 + 421ea66 commit 5a0d0ae
Show file tree
Hide file tree
Showing 6 changed files with 105 additions and 25 deletions.
Binary file modified documentation/APECSS-Documentation.pdf
View file
Binary file not shown.
1 change: 1 addition & 0 deletions documentation/chapters/emissions.tex
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,7 @@ \chapter{Acoustic emissions}
& {\tt EmissionIntegration Euler} & Integrates the radial position and, if applicable, the velocity using an Euler scheme.\\
& {\tt EmissionIntegration RK4} & Integrates the radial position and, if applicable, the velocity using a conventional fourth-order Runge-Kutta scheme. This is the default.\\
& {\tt KBIterTolerance <float>} & Tolerance $\eta$ for the evaluation of the pressure using a model based on the Kirkwood-Bethe hypothesis in conjunction with the NASG EoS.\\
& {\tt PruneEmissions} & Prune the list of emission nodes according to a test function provided by the user. This test function must be hooked up to the function pointer {\tt Bubble->Emissions->prune\_test()}.\\
\hline
\end{tabular} \vspace{0.2em}

Expand Down
7 changes: 3 additions & 4 deletions examples/acousticemitter/plot_result_spherical.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,16 +28,15 @@
ax1.set_yticks([-0.6,-0.3,0,0.3])
ax1.grid(color='gainsboro', linestyle='-', linewidth=0.5)

ax2.plot((kbres[:,1] - R0)/L, ((R0/kbres[:,1])*np.cos(12*np.pi-k*(kbres[:,1]-R0)-0.5*np.pi-(0.5*np.pi-np.arctan(k*kbres[:,1])))/(np.cos(0.5*np.pi-np.arctan(k*kbres[:,1])))),linestyle='solid', linewidth=2, color='goldenrod',label='Analytical')
ax2.plot((kbres[:,1] - R0)/L, ((kbres[:,2]-1e5)),linestyle='dashed', linewidth=1.5, color='navy',label=r'$p_1/(\rho c)$')
ax2.plot((kbres[:,1] - R0)/L, ((R0*R0/(kbres[:,1]*kbres[:,1]))*(np.cos(12*np.pi-k*(kbres[:,1]-R0)-0.5*np.pi)-np.cos(12*np.pi-k*(kbres[:,1]-R0)-0.5*np.pi-(0.5*np.pi-np.arctan(k*kbres[:,1])))/(np.cos(0.5*np.pi-np.arctan(k*kbres[:,1])))))+((R0/kbres[:,1])*np.cos(12*np.pi-k*(kbres[:,1]-R0)-0.5*np.pi-(0.5*np.pi-np.arctan(k*kbres[:,1])))/(np.cos(0.5*np.pi-np.arctan(k*kbres[:,1])))), linewidth=2, color='goldenrod',label='Analytical')
ax2.plot((kbres[:,1] - R0)/L, (kbres[:,3]*997*1478),linestyle='solid', linewidth=1, color='black',label='Kirkwood-Bethe')
ax2.set(xlabel=r'${(r-R)}/{\lambda_\mathrm{a}}$',ylabel=r'$u_1 \rho c/\Delta p_\mathrm{a}$')
ax2.set(xlabel=r'${(r-R)}/{\lambda_\mathrm{a}}$',ylabel=r'$u \rho c/\Delta p_\mathrm{a}$')
ax2.set_xlim(xmin=0,xmax=4)
ax2.set_ylim(ymin=-0.6,ymax=0.3)
ax2.set_yticks([-0.6,-0.3,0,0.3])
ax2.grid(color='gainsboro', linestyle='-', linewidth=0.5)

ax2.legend(fontsize='small', loc=(-0.41,1.04), ncol=3, frameon=False)
ax2.legend(fontsize='small', loc=(0,1.04), ncol=3, frameon=False)

ax1.xaxis.set_label_coords(0.5,-0.19)
ax2.xaxis.set_label_coords(0.5,-0.19)
Expand Down
16 changes: 11 additions & 5 deletions include/apecss.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -264,11 +264,15 @@ struct APECSS_Emissions
struct APECSS_EmissionNode *FirstNode; // First node of the linked list
struct APECSS_EmissionNode *LastNode; // Last node of the linked list

// Pruning the list of emission nodes
int pruneList; // Flag indicating whether the linked list is pruned
int (*prune_test)(struct APECSS_EmissionNode *Node);

// Pointer to the function advancing the emission nodes
int (*advance)(struct APECSS_Bubble *Bubble);

// Pointer to the function that computes the invariant f
APECSS_FLOAT (*compute_f)(struct APECSS_Bubble *Bubble, APECSS_FLOAT g, APECSS_FLOAT pL, APECSS_FLOAT rhoL);
APECSS_FLOAT (*compute_f)(struct APECSS_Bubble *Bubble, struct APECSS_EmissionNode *Node);

// Pointer to the function integrating the radial position and velocity along the outgoing characteristic
int (*integrate_along_characteristic)(struct APECSS_Bubble *Bubble, struct APECSS_EmissionNode *Current, APECSS_FLOAT hinf);
Expand Down Expand Up @@ -514,6 +518,8 @@ int apecss_emissions_freelinkedlist(struct APECSS_Bubble *Bubble);
int apecss_emissions_updatenone(struct APECSS_Bubble *Bubble);
int apecss_emissions_updatelinkedlist(struct APECSS_Bubble *Bubble);
int apecss_emissions_addnode(struct APECSS_Bubble *Bubble);
int apecss_emissions_prunelist(struct APECSS_Bubble *Bubble);
int apecss_emissions_prune_no_node(struct APECSS_EmissionNode *Node);
int apecss_emissions_removenode(struct APECSS_Bubble *Bubble);
int apecss_emissions_advance_finitespeedincompressible(struct APECSS_Bubble *Bubble);
int apecss_emissions_advance_quasiacoustic(struct APECSS_Bubble *Bubble);
Expand All @@ -527,10 +533,10 @@ int apecss_emissions_integrate_ev_general_euler(struct APECSS_Bubble *Bubble, st
int apecss_emissions_integrate_ev_general_rk4(struct APECSS_Bubble *Bubble, struct APECSS_EmissionNode *Current, APECSS_FLOAT hinf);
int apecss_emissions_integrate_tiv_general_euler(struct APECSS_Bubble *Bubble, struct APECSS_EmissionNode *Current, APECSS_FLOAT hinf);
int apecss_emissions_integrate_tiv_general_rk4(struct APECSS_Bubble *Bubble, struct APECSS_EmissionNode *Current, APECSS_FLOAT hinf);
APECSS_FLOAT apecss_emissions_f_zero(struct APECSS_Bubble *Bubble, APECSS_FLOAT g, APECSS_FLOAT pL, APECSS_FLOAT rhoL);
APECSS_FLOAT apecss_emissions_f_finitespeedincompressible(struct APECSS_Bubble *Bubble, APECSS_FLOAT g, APECSS_FLOAT pL, APECSS_FLOAT rhoL);
APECSS_FLOAT apecss_emissions_f_quasiacoustic(struct APECSS_Bubble *Bubble, APECSS_FLOAT g, APECSS_FLOAT pL, APECSS_FLOAT rhoL);
APECSS_FLOAT apecss_emissions_f_kirkwoodbethe(struct APECSS_Bubble *Bubble, APECSS_FLOAT g, APECSS_FLOAT pL, APECSS_FLOAT rhoL);
APECSS_FLOAT apecss_emissions_f_zero(struct APECSS_Bubble *Bubble, struct APECSS_EmissionNode *Node);
APECSS_FLOAT apecss_emissions_f_finitespeedincompressible(struct APECSS_Bubble *Bubble, struct APECSS_EmissionNode *Node);
APECSS_FLOAT apecss_emissions_f_quasiacoustic(struct APECSS_Bubble *Bubble, struct APECSS_EmissionNode *Node);
APECSS_FLOAT apecss_emissions_f_kirkwoodbethe(struct APECSS_Bubble *Bubble, struct APECSS_EmissionNode *Node);

// ---------------------
// gas.c
Expand Down
15 changes: 14 additions & 1 deletion src/bubble.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -233,6 +233,10 @@ int apecss_bubble_readoptions(struct APECSS_Bubble *Bubble, char *OptionsDir)
}
}
}
else if (strncasecmp(option2, "pruneemissions", 14) == 0)
{
Bubble->Emissions->pruneList = 1;
}
else if (strncasecmp(option2, "kbitertolerance", 15) == 0)
{
if (Bubble->Emissions == NULL) apecss_emissions_initializestruct(Bubble);
Expand Down Expand Up @@ -814,7 +818,16 @@ int apecss_bubble_initialize(struct APECSS_Bubble *Bubble)
// ---------------------------------------
// Emissions

if (Bubble->Emissions != NULL) Bubble->Emissions->CutOffDistance = APECSS_MAX(Bubble->Emissions->CutOffDistance, Bubble->R0);
if (Bubble->Emissions != NULL)
{
Bubble->Emissions->CutOffDistance = APECSS_MAX(Bubble->Emissions->CutOffDistance, Bubble->R0);

if (Bubble->Emissions->pruneList && Bubble->Emissions->prune_test == NULL)
{
apecss_erroronscreen(0, "Prune emissions option invoked but no test function defined. No nodes will be pruned!");
Bubble->Emissions->prune_test = apecss_emissions_prune_no_node;
}
}

// ---------------------------------------
// Results
Expand Down
91 changes: 76 additions & 15 deletions src/emissions.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,8 @@ int apecss_emissions_initializestruct(struct APECSS_Bubble *Bubble)
Bubble->Emissions->CutOffDistance = 1.0e-3;
Bubble->Emissions->KB_IterTolerance = 1.0;
Bubble->Emissions->nNodes = 0;
Bubble->Emissions->pruneList = 0;
Bubble->Emissions->prune_test = NULL;
Bubble->Emissions->FirstNode = NULL;
Bubble->Emissions->LastNode = NULL;
Bubble->Emissions->advance = NULL;
Expand Down Expand Up @@ -81,6 +83,7 @@ int apecss_emissions_updatelinkedlist(struct APECSS_Bubble *Bubble)
if (Bubble->Emissions->nNodes) Bubble->Emissions->advance(Bubble);
apecss_emissions_addnode(Bubble);
if (Bubble->Emissions->LastNode->r > Bubble->Emissions->CutOffDistance) apecss_emissions_removenode(Bubble);
if (Bubble->Emissions->pruneList) apecss_emissions_prunelist(Bubble);

return (0);
}
Expand All @@ -98,16 +101,16 @@ int apecss_emissions_addnode(struct APECSS_Bubble *Bubble)
APECSS_FLOAT pinf = Bubble->get_pressure_infinity(Bubble->t, Bubble);
APECSS_FLOAT hinf = Bubble->Liquid->get_enthalpy(pinf, Bubble->Liquid->get_density(pinf, Bubble->Liquid), Bubble->Liquid);

// Compute the invariants f and g
New->g = Bubble->get_dimensionalradius(Bubble->R) * (hL - hinf + 0.5 * APECSS_POW2(Bubble->U));
New->f = Bubble->Emissions->compute_f(Bubble, New->g, pL, rhoL);

// Apply the boundary conditions
New->r = Bubble->R;
New->u = Bubble->U;
New->h = hL;
New->p = pL;

// Compute the invariants f and g
New->g = Bubble->get_dimensionalradius(Bubble->R) * (hL - hinf + 0.5 * APECSS_POW2(Bubble->U));
New->f = Bubble->Emissions->compute_f(Bubble, New); // Requires knowledge of invariant g

// Set the neighbor information
if (Bubble->Emissions->nNodes)
{
Expand All @@ -130,6 +133,32 @@ int apecss_emissions_addnode(struct APECSS_Bubble *Bubble)
return (0);
}

int apecss_emissions_prunelist(struct APECSS_Bubble *Bubble)
{
if (Bubble->Emissions->nNodes > 2) // The list needs to consist of at least 3 nodes.
{
struct APECSS_EmissionNode *Current = Bubble->Emissions->LastNode->backward;

while (Current->backward != NULL)
{
if (Bubble->Emissions->prune_test(Current))
{
struct APECSS_EmissionNode *Obsolete = Current;
Current->backward->forward = Current->forward;
Current->forward->backward = Current->backward;
Current = Current->backward;
free(Obsolete);
}
else
{
Current = Current->backward; // Move to the next node
}
}
}

return (0);
}

int apecss_emissions_removenode(struct APECSS_Bubble *Bubble)
{
struct APECSS_EmissionNode *Obsolete;
Expand Down Expand Up @@ -271,6 +300,8 @@ int apecss_emissions_advance_kirkwoodbethe_tait(struct APECSS_Bubble *Bubble)
// Shock treatment (if necessary)

int check_for_shocks = 1;
int shocks_detected = 0;

while (check_for_shocks)
{
Current = Bubble->Emissions->FirstNode;
Expand All @@ -281,14 +312,14 @@ int apecss_emissions_advance_kirkwoodbethe_tait(struct APECSS_Bubble *Bubble)
if (Current->forward != NULL && Current->r > Current->forward->r) // Check for shock formation
{
check_for_shocks = 1; // Requires another complete subsequent iteration to make sure all multivalued solutions have been or are treated
shocks_detected = 1; // A multivalued solution associated with a shock has been detected.

// Define new properties of the forward node
Current->forward->r = 0.5 * (Current->r + Current->forward->r);
Current->forward->u = 0.5 * (Current->u + Current->forward->u);
Current->forward->g = 0.5 * (Current->g + Current->forward->g);
Current->forward->h = hinf + Current->forward->g / Bubble->get_dimensionalradius(Current->forward->r) - 0.5 * APECSS_POW2(Current->forward->u);
Current->forward->p = APECSS_POW(h_fac * Current->h, h_exp) - B;
Current->forward->f = Bubble->Emissions->compute_f(Bubble, Current->g, Current->p, Bubble->Liquid->get_density(Current->p, Bubble->Liquid));

// Current node is obsolete and will be deleted
struct APECSS_EmissionNode *Obsolete = Current;
Expand All @@ -315,6 +346,20 @@ int apecss_emissions_advance_kirkwoodbethe_tait(struct APECSS_Bubble *Bubble)
}
}

// ---------------------------------------
// Update invariant f for the explicit velocity expression (if applicable)

if (shocks_detected && Bubble->Emissions->Type == APECSS_EMISSION_EV)
{
Current = Bubble->Emissions->FirstNode;

while (Current != NULL)
{
Current->f = Bubble->Emissions->compute_f(Bubble, Current);
Current = Current->forward;
}
}

// ---------------------------------------
// Store data (if applicable)

Expand Down Expand Up @@ -376,6 +421,8 @@ int apecss_emissions_advance_kirkwoodbethe_general(struct APECSS_Bubble *Bubble)
// Shock treatment (if necessary)

int check_for_shocks = 1;
int shocks_detected = 0;

while (check_for_shocks)
{
Current = Bubble->Emissions->FirstNode;
Expand All @@ -386,6 +433,7 @@ int apecss_emissions_advance_kirkwoodbethe_general(struct APECSS_Bubble *Bubble)
if (Current->forward != NULL && Current->r > Current->forward->r) // Check for shock formation
{
check_for_shocks = 1; // Requires another complete subsequent iteration to make sure all multivalued solutions have been or are treated
shocks_detected = 1; // A multivalued solution associated with a shock has been detected.

// Define new properties of the forward node
Current->forward->r = 0.5 * (Current->r + Current->forward->r);
Expand Down Expand Up @@ -437,7 +485,8 @@ int apecss_emissions_advance_kirkwoodbethe_general(struct APECSS_Bubble *Bubble)
Current->p = ((Gamma - 1.0) * rho * Current->h - (1.0 - b * rho) * Gamma * B) / (Gamma - b * rho);
} while (APECSS_ABS((p_prev - Current->p)) > tol * APECSS_ABS(Current->p));

Current->f = Bubble->Emissions->compute_f(Bubble, Current->g, Current->p, Bubble->Liquid->get_density(Current->p, Bubble->Liquid));
// Update invariant f for the explicit velocity expression (if applicable)
if (shocks_detected && Bubble->Emissions->Type == APECSS_EMISSION_EV) Current->f = Bubble->Emissions->compute_f(Bubble, Current);

Current = Current->forward;
}
Expand Down Expand Up @@ -862,22 +911,34 @@ int apecss_emissions_integrate_tiv_general_rk4(struct APECSS_Bubble *Bubble, str
// Bubble->Emissions->compute_f().
// -------------------------------------------------------------------

APECSS_FLOAT apecss_emissions_f_zero(struct APECSS_Bubble *Bubble, APECSS_FLOAT g, APECSS_FLOAT pL, APECSS_FLOAT rhoL) { return (0.0); }
APECSS_FLOAT apecss_emissions_f_zero(struct APECSS_Bubble *Bubble, struct APECSS_EmissionNode *Node) { return (0.0); }

APECSS_FLOAT apecss_emissions_f_finitespeedincompressible(struct APECSS_Bubble *Bubble, APECSS_FLOAT g, APECSS_FLOAT pL, APECSS_FLOAT rhoL)
APECSS_FLOAT apecss_emissions_f_finitespeedincompressible(struct APECSS_Bubble *Bubble, struct APECSS_EmissionNode *Node)
{
return (APECSS_POW2(Bubble->R) * Bubble->U);
return (APECSS_POW2(Node->r) * Node->u);
}

APECSS_FLOAT apecss_emissions_f_quasiacoustic(struct APECSS_Bubble *Bubble, APECSS_FLOAT g, APECSS_FLOAT pL, APECSS_FLOAT rhoL)
APECSS_FLOAT apecss_emissions_f_quasiacoustic(struct APECSS_Bubble *Bubble, struct APECSS_EmissionNode *Node)
{
return (APECSS_POW2(Bubble->R) * Bubble->U - Bubble->R * g / Bubble->Liquid->cref);
return (APECSS_POW2(Node->r) * Node->u - Node->r * Node->g / Bubble->Liquid->cref);
}

APECSS_FLOAT apecss_emissions_f_kirkwoodbethe(struct APECSS_Bubble *Bubble, APECSS_FLOAT g, APECSS_FLOAT pL, APECSS_FLOAT rhoL)
APECSS_FLOAT apecss_emissions_f_kirkwoodbethe(struct APECSS_Bubble *Bubble, struct APECSS_EmissionNode *Node)
{
return (2.0 *
(Bubble->get_dimensionalradius(Bubble->R) * Bubble->R * Bubble->U -
Bubble->R * g / (Bubble->Liquid->get_soundspeed(pL, rhoL, Bubble->Liquid) + Bubble->U)) /
(Bubble->get_dimensionalradius(Node->r) * Node->r * Node->u -
Node->r * Node->g / (Bubble->Liquid->get_soundspeed(Node->p, Bubble->Liquid->get_density(Node->p, Bubble->Liquid), Bubble->Liquid) + Node->u)) /
(Bubble->dimensionality + APECSS_SMALL));
}
}

// -------------------------------------------------------------------
// PRUNING OF EMISSION NODES
// -------------------------------------------------------------------
// Dummy test function for the pruning of emission nodes.
// -------------------------------------------------------------------
// This function is hooked up to the function pointer
// Bubble->Emissions->prune_test() if no user-defined test function
// has been defined. No nodes will be pruned.
// -------------------------------------------------------------------

int apecss_emissions_prune_no_node(struct APECSS_EmissionNode *Node) { return (0); }

0 comments on commit 5a0d0ae

Please sign in to comment.