Skip to content

Commit ce186f5

Browse files
ESeNonFossiIoESeNonFossiIo
authored
Merge pull request #209 from ESeNonFossiIo/stokes_main
Stokes main
2 parents 2a2df45 + 639f653 commit ce186f5

35 files changed

+29294
-28027
lines changed

‎include/interfaces/navier_stokes.h‎

Lines changed: 40 additions & 72 deletions
Original file line numberDiff line numberDiff line change
@@ -23,25 +23,6 @@
2323
* \f]
2424
* can be recoverd setting @p dynamic = false
2525
*
26-
* Dynamic Stokes Equation:
27-
* \f[
28-
* \begin{cases}
29-
* \partial_t u - \nu\textrm{div} \epsilon(u)
30-
* + \frac{1}{\rho}\nabla p = f \\
31-
* \textrm{div}u=0
32-
* \end{cases}
33-
* \f]
34-
* can be recoverd setting @p convection = false
35-
*
36-
* Stokes Equation:
37-
* \f[
38-
* \begin{cases}
39-
* - \nu\textrm{div} \epsilon(u)
40-
* + \frac{1}{\rho}\nabla p = f \\
41-
* \textrm{div}u=0
42-
* \end{cases}
43-
* \f]
44-
* can be recoverd setting @p dynamic = false and @p convection = false
4526
*
4627
* In the code we adopt the following notations:
4728
* - Mp := block resulting from \f$ ( \partial_t p, q ) \f$
@@ -119,7 +100,7 @@ class NavierStokes
119100

120101
public:
121102
~NavierStokes () {};
122-
NavierStokes (bool dynamic, bool convection);
103+
NavierStokes (bool dynamic);
123104

124105
void declare_parameters (ParameterHandler &prm);
125106
void parse_parameters_call_back ();
@@ -186,12 +167,6 @@ class NavierStokes
186167
*/
187168
bool dynamic;
188169

189-
/**
190-
* Enable convection term: \f$ (\nabla u)u \f$.
191-
* This term introduce the non-linearity of Navier Stokes Equations.
192-
*/
193-
bool convection;
194-
195170
/**
196171
* TODO:
197172
*/
@@ -282,7 +257,7 @@ class NavierStokes
282257

283258
template <int dim, int spacedim, typename LAC>
284259
NavierStokes<dim,spacedim, LAC>::
285-
NavierStokes(bool dynamic, bool convection)
260+
NavierStokes(bool dynamic)
286261
:
287262
PDESystemInterface<dim,spacedim,NavierStokes<dim,spacedim,LAC>, LAC>(
288263
"Navier Stokes Interface",
@@ -295,7 +270,6 @@ NavierStokes(bool dynamic, bool convection)
295270
AMG_Ap("AMG for Ap"),
296271
jac_Mp("Jacobi for Mp", 1.4),
297272
dynamic(dynamic),
298-
convection(convection),
299273
compute_Mp(false),
300274
compute_Ap(false),
301275
is_parallel(Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD) > 1)
@@ -333,10 +307,6 @@ declare_parameters (ParameterHandler &prm)
333307
"Enable dynamic term (\\partial_t u)", "true",
334308
Patterns::Bool(),
335309
"Enable the dynamic term of the equation.");
336-
this->add_parameter(prm, &convection,
337-
"Enable convection term ((\\nabla u)u)", "true",
338-
Patterns::Bool(),
339-
"Enable the convection term of the equation. Set it false if you want to solve Stokes Equation.");
340310
this->add_parameter(prm, &non_linear_term, "Non linear term","linear",
341311
Patterns::Selection("fully_non_linear|linear|RHS"),
342312
"Available options: \n"
@@ -415,8 +385,8 @@ solution_preprocessing(FEValuesCache<dim,spacedim> &fe_cache) const
415385
{
416386
this->reinit(dummy, cell, face, fe_cache);
417387
// Velocity:
418-
auto &sym_grad_u_ = fe_cache.get_symmetric_gradients( "explicit_solution", "grad_u", velocity, dummy);
419-
auto &p_ = fe_cache.get_values( "explicit_solution", "p", pressure, dummy);
388+
auto &sym_grad_u_ = fe_cache.get_symmetric_gradients( "explicit_solution", "f_grad_u", velocity, dummy);
389+
auto &p_ = fe_cache.get_values( "explicit_solution", "f_p", pressure, dummy);
420390

421391
auto &fev = fe_cache.get_current_fe_values();
422392
auto &q_points = fe_cache.get_quadrature_points();
@@ -565,14 +535,14 @@ energies_and_residuals(const typename DoFHandler<dim,spacedim>::active_cell_iter
565535
auto &grad_ps = fe_cache.get_gradients("solution", "grad_p", pressure,et);
566536

567537
// Previous time step solution:
568-
auto &u_olds = fe_cache.get_values("explicit_solution", "u", velocity, dummy);
569-
auto &grad_u_olds = fe_cache.get_gradients("explicit_solution", "grad_u", velocity, dummy);
538+
auto &u_olds = fe_cache.get_values("explicit_solution", "ue", velocity, dummy);
539+
auto &grad_u_olds = fe_cache.get_gradients("explicit_solution", "grad_ue", velocity, dummy);
570540

571541
// Previous Jacobian step solution:
572542
fe_cache.cache_local_solution_vector("prev_solution",
573543
this->get_locally_relevant_solution(), dummy);
574-
auto &u_prevs = fe_cache.get_values("prev_solution", "u", velocity, dummy);
575-
auto &grad_u_prevs = fe_cache.get_gradients("prev_solution", "grad_u", velocity, dummy);
544+
auto &u_prevs = fe_cache.get_values("prev_solution", "up", velocity, dummy);
545+
auto &grad_u_prevs = fe_cache.get_gradients("prev_solution", "grad_up", velocity, dummy);
576546

577547
const unsigned int n_quad_points = us.size();
578548
auto &JxW = fe_cache.get_JxW_values();
@@ -620,42 +590,40 @@ energies_and_residuals(const typename DoFHandler<dim,spacedim>::active_cell_iter
620590
if (dynamic)
621591
res += rho * u_dot * v;
622592

623-
// Convection:
624-
if (convection)
593+
Tensor<2, dim, ResidualType> gradoldu;
594+
Tensor<1, dim, ResidualType> oldu;
595+
596+
if (linearize_in_time)
625597
{
626-
Tensor<2, dim, ResidualType> gradoldu;
627-
Tensor<1, dim, ResidualType> oldu;
628-
629-
if (linearize_in_time)
630-
{
631-
gradoldu=grad_u_old;
632-
oldu=u_old;
633-
}
634-
else
635-
{
636-
gradoldu=grad_u_prev;
637-
oldu=u_prev;
638-
}
639-
640-
if (non_linear_term=="fully_non_linear")
641-
nl_u = grad_u * u;
642-
else if (non_linear_term=="linear")
643-
nl_u = grad_u * oldu;
644-
else if (non_linear_term=="RHS")
645-
nl_u = gradoldu * oldu;
646-
647-
ResidualType non_linear_term = 0;
648-
649-
if (use_skew_symmetric_advection)
650-
res += 0.5 * ( scalar_product( nl_u, v) + scalar_product(grad_v * oldu, u) );
651-
else
652-
res += scalar_product( nl_u, v);
653-
654-
double norm = std::sqrt(SacadoTools::to_double(oldu*oldu));
655-
656-
if (norm > 0 && SUPG_alpha > 0)
657-
res += scalar_product( nl_u, SUPG_alpha * (h/norm) * grad_v * oldu);
598+
gradoldu=grad_u_old;
599+
oldu=u_old;
658600
}
601+
else
602+
{
603+
gradoldu=grad_u_prev;
604+
oldu=u_prev;
605+
}
606+
607+
if (non_linear_term=="fully_non_linear")
608+
nl_u = grad_u * u;
609+
else if (non_linear_term=="linear")
610+
nl_u = grad_u * oldu;
611+
else if (non_linear_term=="RHS")
612+
nl_u = gradoldu * oldu;
613+
614+
ResidualType non_linear_term = 0;
615+
616+
if (use_skew_symmetric_advection)
617+
res += 0.5 * ( scalar_product( nl_u, v) + scalar_product(grad_v * oldu, u) );
618+
else
619+
res += scalar_product( nl_u, v);
620+
621+
double norm = std::sqrt(SacadoTools::to_double(oldu*oldu));
622+
623+
if (norm > 0 && SUPG_alpha > 0)
624+
res += scalar_product( nl_u, SUPG_alpha * (h/norm) * grad_v * oldu);
625+
626+
659627
// grad-div stabilization term:
660628
if (gamma!=0.0)
661629
res += gamma * div_u * div_v;

‎include/interfaces/stokes.h‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -6,8 +6,8 @@
66

77
#include <deal.II/lac/solver_cg.h>
88
#include <deal2lkit/sacado_tools.h>
9-
#include <deal2lkit/parsed_preconditioner_amg.h>
10-
#include <deal2lkit/parsed_preconditioner_jacobi.h>
9+
#include <deal2lkit/parsed_preconditioner/amg.h>
10+
#include <deal2lkit/parsed_preconditioner/jacobi.h>
1111

1212
using namespace SacadoTools;
1313

‎main/fluid-dynamics.cc‎

Lines changed: 40 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,6 @@
11
#include "interfaces/ALE_navier_stokes.h"
22
#include "interfaces/navier_stokes.h"
3+
#include "interfaces/stokes.h"
34
#include "pidomus.h"
45

56
#include "deal.II/base/numbers.h"
@@ -24,8 +25,16 @@
2425
ParameterAcceptor::initialize(prm_file, pde_name+"_used.prm"); \
2526
equation.run ();
2627

27-
#define problem_NS(dynamic,stokes,dim,spacedim,LAC) \
28-
NavierStokes<dim,spacedim,LAC> energy(dynamic, stokes); \
28+
#define problem_stokes(dim,spacedim,LAC) \
29+
StokesInterface<dim,spacedim,LAC> energy; \
30+
piDoMUS<dim,spacedim,LAC> equation ( \
31+
"piDoMUS", \
32+
energy); \
33+
ParameterAcceptor::initialize(prm_file, pde_name+"_used.prm"); \
34+
equation.run ();
35+
36+
#define problem_NS(dynamic,dim,spacedim,LAC) \
37+
NavierStokes<dim,spacedim,LAC> energy(dynamic); \
2938
piDoMUS<dim,spacedim,LAC> equation ( \
3039
"piDoMUS", \
3140
energy); \
@@ -109,22 +118,18 @@ int main (int argc, char *argv[])
109118
if (dynamic)
110119
string_dynamic="Dynamic";
111120

112-
bool stokes;
113121
std::string string_pde_name="";
114122
if (pde_name == "NS")
115123
{
116124
string_pde_name = "Navier Stokes";
117-
stokes = false;
118125
}
119126
else if (pde_name == "ALE")
120127
{
121128
string_pde_name = "ALE Navier Stokes";
122-
stokes = false;
123129
}
124130
else if (pde_name == "stokes")
125131
{
126132
string_pde_name = "Stokes";
127-
stokes = false;
128133
}
129134
else
130135
{
@@ -170,28 +175,53 @@ int main (int argc, char *argv[])
170175
}
171176
}
172177
}
178+
else if ( pde_name == "stokes" )
179+
{
180+
if (dim==2)
181+
{
182+
if (trilinos)
183+
{
184+
problem_stokes(2,2,LATrilinos);
185+
}
186+
else
187+
{
188+
problem_stokes(2,2,LADealII);
189+
}
190+
}
191+
else
192+
{
193+
if (trilinos)
194+
{
195+
problem_stokes(3,3,LATrilinos);
196+
}
197+
else
198+
{
199+
problem_stokes(3,3,LADealII);
200+
}
201+
}
202+
}
173203
else
174204
{
175205
if (dim==2)
176206
{
177207
if (trilinos)
178208
{
179-
problem_NS(dynamic,stokes,2,2,LATrilinos);
209+
problem_NS(dynamic,2,2,LATrilinos);
180210
}
181211
else
182212
{
183-
problem_NS(dynamic,stokes,2,2,LADealII);
213+
problem_NS(dynamic,2,2,LADealII);
184214
}
185215
}
186216
else
187217
{
188218
if (trilinos)
189219
{
190-
problem_NS(dynamic,stokes,3,3,LATrilinos);
220+
problem_NS(dynamic,3,3,LATrilinos);
191221
}
192222
else
193223
{
194-
problem_NS(dynamic,stokes,3,3,LADealII);
224+
problem_NS(dynamic,3,3,LADealII);
195225
}
196226
}
197227
}

‎tests/dynamic_stokes_00.cc‎

Lines changed: 0 additions & 42 deletions
This file was deleted.

0 commit comments

Comments
 (0)