-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathSetFluidGasSolidViscosity.cpp
188 lines (165 loc) · 7.8 KB
/
SetFluidGasSolidViscosity.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
// ---------------------------------------------------------------------
//
// Copyright (c) 2019 - 2022 by the IBAMR developers
// All rights reserved.
//
//
// IBAMR is free software and is distributed under the 3-clause BSD
// license. The full text of the license can be found in the file
// COPYRIGHT at the top level directory of IBAMR.
//
// ---------------------------------------------------------------------
#include <ibamr/app_namespaces.h>
#include <ibtk/HierarchyMathOps.h>
#include "SetFluidGasSolidViscosity.h"
#include <CartesianGridGeometry.h>
/////////////////////////////// STATIC ///////////////////////////////////////
void
callSetFluidGasSolidViscosityCallbackFunction(int mu_idx,
Pointer<Variable<NDIM> > mu_var,
Pointer<IBTK::HierarchyMathOps> hier_math_ops,
const int cycle_num,
const double time,
const double current_time,
const double new_time,
void* ctx)
{
// Set the density from the level set information
static SetFluidGasSolidViscosity* ptr_SetFluidGasSolidViscosity = static_cast<SetFluidGasSolidViscosity*>(ctx);
ptr_SetFluidGasSolidViscosity->setViscosityPatchData(
mu_idx, mu_var, hier_math_ops, cycle_num, time, current_time, new_time);
return;
} // callSetFluidGasSolidViscosityCallBackFunction
/////////////////////////////// PUBLIC //////////////////////////////////////
SetFluidGasSolidViscosity::SetFluidGasSolidViscosity(const std::string& object_name,
Pointer<AdvDiffHierarchyIntegrator> adv_diff_solver,
Pointer<CellVariable<NDIM, double> > ls_solid_var,
Pointer<CellVariable<NDIM, double> > ls_gas_var,
const double mu_fluid,
const double mu_gas,
const double mu_solid,
const int ls_reinit_interval,
const double num_solid_interface_cells,
const double num_gas_interface_cells,
const bool set_mu_solid)
: d_object_name(object_name),
d_adv_diff_solver(adv_diff_solver),
d_ls_solid_var(ls_solid_var),
d_ls_gas_var(ls_gas_var),
d_mu_fluid(mu_fluid),
d_mu_gas(mu_gas),
d_mu_solid(mu_solid),
d_ls_reinit_interval(ls_reinit_interval),
d_num_solid_interface_cells(num_solid_interface_cells),
d_num_gas_interface_cells(num_gas_interface_cells),
d_set_mu_solid(set_mu_solid)
{
// intentionally left blank
return;
} // SetFluidGasSolidViscosity
void
SetFluidGasSolidViscosity::setViscosityPatchData(int mu_idx,
Pointer<Variable<NDIM> > mu_var,
Pointer<HierarchyMathOps> hier_math_ops,
const int /*cycle_num*/,
const double time,
const double current_time,
const double new_time)
{
// Get the current level set information
VariableDatabase<NDIM>* var_db = VariableDatabase<NDIM>::getDatabase();
int ls_solid_idx = -1, ls_gas_idx = -1;
if (MathUtilities<double>::equalEps(time, current_time))
{
ls_solid_idx = var_db->mapVariableAndContextToIndex(d_ls_solid_var, d_adv_diff_solver->getCurrentContext());
}
else if (MathUtilities<double>::equalEps(time, new_time))
{
ls_solid_idx = var_db->mapVariableAndContextToIndex(d_ls_solid_var, d_adv_diff_solver->getNewContext());
}
else
{
TBOX_ERROR("This statement should not be reached");
}
if (MathUtilities<double>::equalEps(time, current_time))
{
ls_gas_idx = var_db->mapVariableAndContextToIndex(d_ls_gas_var, d_adv_diff_solver->getCurrentContext());
}
else if (MathUtilities<double>::equalEps(time, new_time))
{
ls_gas_idx = var_db->mapVariableAndContextToIndex(d_ls_gas_var, d_adv_diff_solver->getNewContext());
}
else
{
TBOX_ERROR("This statement should not be reached");
}
Pointer<PatchHierarchy<NDIM> > patch_hierarchy = hier_math_ops->getPatchHierarchy();
const int coarsest_ln = 0;
const int finest_ln = patch_hierarchy->getFinestLevelNumber();
Pointer<CellVariable<NDIM, double> > mu_cc_var = mu_var;
if (mu_cc_var)
{
for (int ln = coarsest_ln; ln <= finest_ln; ++ln)
{
Pointer<PatchLevel<NDIM> > level = patch_hierarchy->getPatchLevel(ln);
for (PatchLevel<NDIM>::Iterator p(level); p; p++)
{
Pointer<Patch<NDIM> > patch = level->getPatch(p());
Pointer<CartesianPatchGeometry<NDIM> > patch_geom = patch->getPatchGeometry();
const double* const patch_dx = patch_geom->getDx();
const double alpha = d_num_solid_interface_cells * patch_dx[0];
const double beta = d_num_gas_interface_cells * patch_dx[1];
const Box<NDIM>& patch_box = patch->getBox();
const Pointer<CellData<NDIM, double> > ls_solid_data = patch->getPatchData(ls_solid_idx);
const Pointer<CellData<NDIM, double> > ls_gas_data = patch->getPatchData(ls_gas_idx);
Pointer<CellData<NDIM, double> > mu_data = patch->getPatchData(mu_idx);
// Calderer et al, 2014
for (Box<NDIM>::Iterator it(patch_box); it; it++)
{
CellIndex<NDIM> ci(it());
const double phi_s = (*ls_solid_data)(ci);
const double phi_g = (*ls_gas_data)(ci);
double Hphi_s, Hphi_g;
if (phi_s < -alpha)
{
Hphi_s = 0.0;
}
else if (std::abs(phi_s) <= alpha)
{
Hphi_s = 0.5 + 0.5 * phi_s / alpha + 1.0 / (2.0 * M_PI) * std::sin(M_PI * phi_s / alpha);
}
else
{
Hphi_s = 1.0;
}
if (phi_g < -beta)
{
Hphi_g = 0.0;
}
else if (std::abs(phi_g) <= beta)
{
Hphi_g = 0.5 + 0.5 * phi_g / beta + 1.0 / (2.0 * M_PI) * std::sin(M_PI * phi_g / beta);
}
else
{
Hphi_g = 1.0;
}
// First, compute the viscosity of the "flowing" phases
const double mu_flow = (d_mu_fluid - d_mu_gas) * Hphi_g + d_mu_gas;
// Next, set the viscosity of the solid phase in the usual way
if (d_set_mu_solid)
(*mu_data)(ci) = (mu_flow - d_mu_solid) * Hphi_s + d_mu_solid;
else
(*mu_data)(ci) = mu_flow;
}
}
}
}
else
{
// Erroring out if any other centered is used for mu
TBOX_ERROR("This statement should not have been reached");
}
return;
} // setViscosityPatchData
/////////////////////////////// PRIVATE //////////////////////////////////////