From b9a4d5433b9c07ea500d51e9c3e6335730655779 Mon Sep 17 00:00:00 2001 From: bangerth Date: Mon, 12 May 2008 19:45:18 +0000 Subject: [PATCH] Make it work now. git-svn-id: https://svn.dealii.org/trunk@16079 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-33/step-33.cc | 168 ++++++++++++++-------------- 1 file changed, 83 insertions(+), 85 deletions(-) diff --git a/deal.II/examples/step-33/step-33.cc b/deal.II/examples/step-33/step-33.cc index 9726a40b87..3a24155736 100644 --- a/deal.II/examples/step-33/step-33.cc +++ b/deal.II/examples/step-33/step-33.cc @@ -213,6 +213,7 @@ struct EulerEquations // is the basic Lax-Friedrich's flux with a stabilization parameter // $\alpha$. template + static void LFNumFlux(std::vector > > &nflux, const std::vector > &points, const std::vector > &normals, @@ -280,8 +281,8 @@ class InitialCondition : public FunctionParser template InitialCondition::InitialCondition () : - FunctionParser (n_components()), - expressions(n_components(), "0.0") + FunctionParser (EulerEquations::n_components), + expressions(EulerEquations::n_components, "0.0") {} // Here we set up x,y,z as the variables that one should use in the input @@ -577,11 +578,6 @@ void ConsLaw::assemble_cell_term( unsigned int /*cell_no*/ ) { - // The residual for each row (i) will be accumulating - // into this fad variable. At the end of the assembly - // for this row, we will query for the sensitivities - // to this variable and add them into the Jacobian. - Sacado::Fad::DFad F_i; unsigned int dofs_per_cell = fe_v.dofs_per_cell; unsigned int n_q_points = fe_v.n_quadrature_points; @@ -590,18 +586,18 @@ void ConsLaw::assemble_cell_term( // Values of the conservative variables at the quadrature points. std::vector > > W (n_q_points, - std::vector >(n_components())); + std::vector >(EulerEquations::n_components)); // Values at the last time step of the conservative variables. // Note that these do not use fad variables, since they do // not depend on the 'variables to be sought'=DOFS. std::vector > Wl (n_q_points, - std::vector(n_components())); + std::vector(EulerEquations::n_components)); // Here we will hold the averaged values of the conservative // variables that we will linearize around (cn=Crank Nicholson). std::vector > > Wcn (n_q_points, - std::vector >(n_components())); + std::vector >(EulerEquations::n_components)); // Gradients of the current variables. It is a // bit of a shame that we have to compute these; we almost don't. @@ -609,13 +605,10 @@ void ConsLaw::assemble_cell_term( // the flux doesn't generally involve any gradients. We do // need these, however, for the diffusion stabilization. std::vector > > > Wgrads (n_q_points, - std::vector > >(n_components(), + std::vector > >(EulerEquations::n_components, std::vector >(dim))); - const std::vector &JxW = fe_v.get_JxW_values (); - - // Here is the magical point where we declare a subset // of the fad variables as degrees of freedom. All // calculations that reference these variables (either @@ -633,7 +626,7 @@ void ConsLaw::assemble_cell_term( // fad types, only the local cell variables, we explicitly // code this loop; for (unsigned int q = 0; q < n_q_points; q++) { - for (unsigned int di = 0; di < n_components(); di++) { + for (unsigned int di = 0; di < EulerEquations::n_components; di++) { W[q][di] = 0; Wl[q][di] = 0; Wcn[q][di] = 0; @@ -666,10 +659,10 @@ void ConsLaw::assemble_cell_term( // could be a rather large object, but for now it // seems to work just fine. typedef Sacado::Fad::DFad FluxMatrix[EulerEquations::n_components][dim]; - std::vector flux(n_q_points); - + FluxMatrix *flux = new FluxMatrix[n_q_points]; + for (unsigned int q=0; q < n_q_points; ++q) - flux_matrix(flux[q], Wcn[q]); + EulerEquations::flux_matrix(flux[q], Wcn[q]); // We now have all of the function values/grads/fluxes, @@ -689,57 +682,63 @@ void ConsLaw::assemble_cell_term( const unsigned int component_i = fe_v.get_fe().system_to_component_index(i).first; - // Initialize the fad residual to zero (removes - // any previous sensitivities. - F_i = 0; - - // Loop quadrature points. - for (unsigned int point=0; point F_i; - Sacado::Fad::DFad fdotgv = 0; - - // Integrate the flux times gradient of the test function - for (unsigned int d = 0; d < dim; d++) - fdotgv += flux[point][component_i][d]*fe_v.shape_grad_component(i, point, component_i)[d]; - - F_i -= fdotgv*JxW[point]; - - // The mass term (if the simulation is non-stationary). - Sacado::Fad::DFad delta_t= 1.0/dT*(W[point][component_i] - Wl[point][component_i]); - if (!is_stationary) F_i += delta_t* - fe_v.shape_value_component(i, point, component_i)*JxW[point]; - - // Stabilization (cell wise diffusion) - Sacado::Fad::DFad guv = 0; - for (unsigned int d = 0; d < dim; d++) { - guv += fe_v.shape_grad_component(i, point, component_i)[d]* - Wgrads[point][component_i][d]; - } - - F_i += 1.0*std::pow(cell_diameter, diffusion_power)*guv*JxW[point]; + for (unsigned int point=0; point()]*fe_v.shape_value_component(i,point, component_i)*JxW[point]; - } else if (component_i == energy_component()) { - F_i += gravity*Wcn[point][density_component()]*Wcn[point][dim-1]* - fe_v.shape_value_component(i,point, component_i)*JxW[point]; + // The gravity component only enters into the energy + // equation and into the vertical component of the + // velocity. + if (component_i == dim - 1) + F_i += gravity * + Wcn[point][EulerEquations::density_component] * + fe_v.shape_value_component(i,point, component_i) * + fe_v.JxW(point); + else if (component_i == EulerEquations::energy_component) + F_i += gravity * + Wcn[point][EulerEquations::density_component] * + Wcn[point][dim-1] * + fe_v.shape_value_component(i,point, component_i) * + fe_v.JxW(point); } - } // for q // Here we gain access to the array of sensitivities // of the residual. We then sum these into the // Epetra matrix. double *values = &(F_i.fastAccessDx(0)); Matrix->SumIntoGlobalValues(dofs[i], - dofs_per_cell, &values[0], reinterpret_cast(&dofs[0])); + dofs_per_cell, + values, + reinterpret_cast(&dofs[0])); // Add minus the residual to the right hand side. right_hand_side(dofs[i]) -= F_i.val(); + } - } // for i + delete[] flux; } // @sect4{%Function: assemble_face_term} // These are either @@ -771,12 +770,11 @@ void ConsLaw::assemble_face_term( // The conservative variables for this cell, // and for std::vector > > Wplus (n_q_points, - std::vector >(n_components())); + std::vector >(EulerEquations::n_components)); std::vector > > Wminus (n_q_points, - std::vector >(n_components())); + std::vector >(EulerEquations::n_components)); - const std::vector &JxW = fe_v.get_JxW_values (); const std::vector > &normals = fe_v.get_normal_vectors (); @@ -800,7 +798,7 @@ void ConsLaw::assemble_face_term( // Set the values of the local conservative variables. // Initialize all variables to zero. for (unsigned int q = 0; q < n_q_points; q++) { - for (unsigned int di = 0; di < n_components(); di++) { + for (unsigned int di = 0; di < EulerEquations::n_components; di++) { Wplus[q][di] = 0; Wminus[q][di] = 0; } @@ -840,13 +838,13 @@ void ConsLaw::assemble_face_term( // and implicit values. If a particular component is not // prescribed, the values evaluate to zero and are // ignored, below. - std::vector > bvals(n_q_points, Vector(n_components())); + std::vector > bvals(n_q_points, Vector(EulerEquations::n_components)); bme->second.second->vector_value_list(fe_v.get_quadrature_points(), bvals); // We loop the quadrature points, and we treat each // component individualy. for (unsigned int q = 0; q < n_q_points; q++) { - for (unsigned int di = 0; di < n_components(); di++) { + for (unsigned int di = 0; di < EulerEquations::n_components; di++) { // An inflow/dirichlet type of boundary condition if (bme->second.first[di] == INFLOW_BC) { @@ -863,8 +861,8 @@ void ConsLaw::assemble_face_term( Sacado::Fad::DFad rho_vel_sqr = 0; Sacado::Fad::DFad dens; - dens = bme->second.first[density_component()] == INFLOW_BC ? bvals[q](density_component()) : - Wplus[q][density_component()]; + dens = bme->second.first[EulerEquations::density_component] == INFLOW_BC ? bvals[q](EulerEquations::density_component) : + Wplus[q][EulerEquations::density_component]; for (unsigned int d=0; d < dim; d++) { if (bme->second.first[d] == INFLOW_BC) @@ -875,7 +873,7 @@ void ConsLaw::assemble_face_term( rho_vel_sqr /= dens; // Finally set the energy value as determined by the // prescribed pressure and the other variables. - Wminus[q][di] = bvals[q](di)/(gas_gamma-1.0) + + Wminus[q][di] = bvals[q](di)/(EulerEquations::gas_gamma-1.0) + 0.5*rho_vel_sqr; } else if (bme->second.first[di] == OUTFLOW_BC) { @@ -902,7 +900,7 @@ void ConsLaw::assemble_face_term( // Determine the Lax-Friedrich's stability parameter, // and evaluate the numerical flux function at the quadrature points - std::vector > > nflux(n_q_points, std::vector >(n_components(), 0)); + std::vector > > nflux(n_q_points, std::vector >(EulerEquations::n_components, 0)); double alpha = 1; switch(flux_params.LF_stab) { @@ -914,7 +912,7 @@ void ConsLaw::assemble_face_term( break; } - LFNumFlux, dim>(nflux, fe_v.get_quadrature_points(), normals, Wplus, Wminus, + EulerEquations::LFNumFlux(nflux, fe_v.get_quadrature_points(), normals, Wplus, Wminus, alpha); // Now assemble the face term @@ -926,7 +924,7 @@ void ConsLaw::assemble_face_term( const unsigned int component_i = fe_v.get_fe().system_to_component_index(i).first; - F_i += nflux[point][component_i]*fe_v.shape_value_component(i, point, component_i)*JxW[point]; + F_i += nflux[point][component_i]*fe_v.shape_value_component(i, point, component_i)*fe_v.JxW(point); } @@ -1216,7 +1214,7 @@ ConsLaw::ConsLaw () // continuous finite elements. The choice was made here. template void ConsLaw::build_fe() { - fe_ptr = new FESystem(FE_Q(1), n_components()); + fe_ptr = new FESystem(FE_Q(1), EulerEquations::n_components); } // Bye bye Conservation law. @@ -1462,11 +1460,11 @@ void ConsLaw::postprocess() { mapping, *fe_ptr, unit_support, update_flags1); std::vector > U(n_uq_points, - Vector(n_components())); + Vector(EulerEquations::n_components)); std::vector > UU(n_q_points, - Vector(n_components())); + Vector(EulerEquations::n_components)); std::vector > > dU(n_uq_points, - std::vector >(n_components())); + std::vector >(EulerEquations::n_components)); typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), @@ -1483,8 +1481,8 @@ void ConsLaw::postprocess() { fe_v.get_function_values(solution, UU); for (unsigned int q = 0; q < fe_v.get_fe().base_element(0).n_dofs_per_cell(); q++) { - unsigned int didx = fe_v.get_fe().component_to_system_index(density_component(), q); - unsigned int eidx = fe_v.get_fe().component_to_system_index(energy_component(), q); + unsigned int didx = fe_v.get_fe().component_to_system_index(EulerEquations::density_component, q); + unsigned int eidx = fe_v.get_fe().component_to_system_index(EulerEquations::energy_component, q); double rho_normVsqr = 0; for (unsigned int d = 0; d < dim; d++) { unsigned int vidx = fe_v.get_fe().component_to_system_index(d, q); @@ -1493,7 +1491,7 @@ void ConsLaw::postprocess() { } rho_normVsqr /= solution(dofs[didx]); // Pressure - ppsolution(dofs[eidx]) = (gas_gamma-1.0)*(solution(dofs[eidx]) - 0.5*rho_normVsqr); + ppsolution(dofs[eidx]) = (EulerEquations::gas_gamma-1.0)*(solution(dofs[eidx]) - 0.5*rho_normVsqr); // Either output density or gradient squared of density, // depending on what the user wants. @@ -1501,7 +1499,7 @@ void ConsLaw::postprocess() { ppsolution(dofs[didx]) = solution(dofs[didx]); } else { double ng = 0; - for (unsigned int i = 0; i < dim; i++) ng += dU[q][density_component()][i]*dU[q][density_component()][i]; + for (unsigned int i = 0; i < dim; i++) ng += dU[q][EulerEquations::density_component][i]*dU[q][EulerEquations::density_component][i]; ng = std::sqrt(ng); ppsolution(dofs[didx]) = ng; } @@ -1532,9 +1530,9 @@ void ConsLaw::estimate() { mapping, *fe_ptr, quadrature_formula, update_flags); std::vector > U(n_q_points, - Vector(n_components())); + Vector(EulerEquations::n_components)); std::vector > > dU(n_q_points, - std::vector >(n_components())); + std::vector >(EulerEquations::n_components)); typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), @@ -1548,7 +1546,7 @@ void ConsLaw::estimate() { indicator(cell_no) = 0; for (unsigned int q = 0; q < n_q_points; q++) { double ng = 0; - for (unsigned int d = 0; d < dim; d++) ng += dU[q][density_component()][d]*dU[q][density_component()][d]; + for (unsigned int d = 0; d < dim; d++) ng += dU[q][EulerEquations::density_component][d]*dU[q][EulerEquations::density_component][d]; indicator(cell_no) += std::log(1+std::sqrt(ng)); @@ -1710,7 +1708,7 @@ void ConsLaw::declare_parameters() { ""); // declare a slot for each of the conservative // variables. - for (unsigned int di = 0; di < n_components(); di++) { + for (unsigned int di = 0; di < EulerEquations::n_components; di++) { char var[512]; std::sprintf(var, "w_%d", di); prm.declare_entry(var, "outflow", @@ -1730,7 +1728,7 @@ void ConsLaw::declare_parameters() { // Initial condition block. prm.enter_subsection("initial condition"); - for (unsigned int di = 0; di < n_components(); di++) { + for (unsigned int di = 0; di < EulerEquations::n_components; di++) { char var[512]; std::sprintf(var, "w_%d", di); @@ -1851,11 +1849,11 @@ void ConsLaw::load_parameters(const char *infile){ // The boundary info for (unsigned int b = 0; b < MAX_BD; b++) { - std::vector flags(n_components(), OUTFLOW_BC); + std::vector flags(EulerEquations::n_components, OUTFLOW_BC); // Define a parser for every boundary, though it may be // unused. - SideCondition *sd = new SideCondition(n_components()); + SideCondition *sd = new SideCondition(EulerEquations::n_components); char bd[512]; std::sprintf(bd, "boundary_%d", b); prm.enter_subsection(bd); @@ -1863,7 +1861,7 @@ void ConsLaw::load_parameters(const char *infile){ const std::string &nopen = prm.get("no penetration"); // Determine how each component is handled. - for (unsigned int di = 0; di < n_components(); di++) { + for (unsigned int di = 0; di < EulerEquations::n_components; di++) { char var[512]; std::sprintf(var, "w_%d", di); std::string btype = prm.get(var); @@ -1889,7 +1887,7 @@ void ConsLaw::load_parameters(const char *infile){ // Initial conditions. prm.enter_subsection("initial condition"); - for (unsigned int di = 0; di < n_components(); di++) { + for (unsigned int di = 0; di < EulerEquations::n_components; di++) { char var[512]; std::sprintf(var, "w_%d value", di); -- 2.39.5