From 6fb2a209ca2b055f475df0f54d7a9d9621dfa994 Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Wed, 14 Mar 2018 00:50:27 +0100 Subject: [PATCH] Fix two_phase_flow --- two_phase_flow/LevelSetSolver.cc | 72 ++++++++++++++-------------- two_phase_flow/MultiPhase.cc | 24 +++++----- two_phase_flow/NavierStokesSolver.cc | 33 ++++++------- 3 files changed, 65 insertions(+), 64 deletions(-) diff --git a/two_phase_flow/LevelSetSolver.cc b/two_phase_flow/LevelSetSolver.cc index bf6ceb2..7a63cb3 100644 --- a/two_phase_flow/LevelSetSolver.cc +++ b/two_phase_flow/LevelSetSolver.cc @@ -76,7 +76,7 @@ class LevelSetSolver ///////////////////////// // BOUNDARY CONDITIONS // ///////////////////////// - void set_boundary_conditions(std::vector boundary_values_id_u, + void set_boundary_conditions(std::vector &boundary_values_id_u, std::vector boundary_values_u); ////////////////// // SET VELOCITY // @@ -169,10 +169,10 @@ class LevelSetSolver // MY PETSC WRAPPERS // /////////////////////// void get_vector_values(PETScWrappers::VectorBase &vector, - const std::vector &indices, + const std::vector &indices, std::vector &values); void get_vector_values(PETScWrappers::VectorBase &vector, - const std::vector &indices, + const std::vector &indices, std::map &map_from_Q1_to_Q2, std::vector &values); @@ -202,7 +202,7 @@ class LevelSetSolver std_cxx1x::shared_ptr MC_preconditioner; // BOUNDARIES - std::vector boundary_values_id_u; + std::vector boundary_values_id_u; std::vector boundary_values_u; ////////////// @@ -343,7 +343,7 @@ void LevelSetSolver::initial_condition (PETScWrappers::MPI::Vector un, ////////// BOUNDARY CONDITIONS ////////// ///////////////////////////////////////// template -void LevelSetSolver::set_boundary_conditions(std::vector boundary_values_id_u, +void LevelSetSolver::set_boundary_conditions(std::vector &boundary_values_id_u, std::vector boundary_values_u) { this->boundary_values_id_u = boundary_values_id_u; @@ -597,7 +597,7 @@ void LevelSetSolver::assemble_ML() const unsigned int n_q_points = quadrature_formula.size(); Vector cell_ML (dofs_per_cell); - std::vector local_dof_indices (dofs_per_cell); + std::vector local_dof_indices (dofs_per_cell); typename DoFHandler::active_cell_iterator cell_LS = dof_handler_LS.begin_active(), @@ -650,7 +650,7 @@ void LevelSetSolver::assemble_MC() const unsigned int n_q_points = quadrature_formula.size(); FullMatrix cell_MC (dofs_per_cell, dofs_per_cell); - std::vector local_dof_indices (dofs_per_cell); + std::vector local_dof_indices (dofs_per_cell); std::vector shape_values(dofs_per_cell); typename DoFHandler::active_cell_iterator @@ -810,7 +810,7 @@ void LevelSetSolver::assemble_K_times_vector(PETScWrappers::MPI::Vector &so Vector un_dofs(dofs_per_cell); - std::vector indices_LS (dofs_per_cell); + std::vector indices_LS (dofs_per_cell); // loop on cells typename DoFHandler::active_cell_iterator @@ -860,10 +860,10 @@ void LevelSetSolver::assemble_K_DL_DH_times_vector dLij_matrix = 0; dCij_matrix = 0; - int ncolumns; - const int *gj; - const double *Cxi, *Cyi, *Czi, *CTxi, *CTyi, *CTzi; - const double *EntResi, *SuppSizei, *MCi; + PetscInt ncolumns; + const PetscInt *gj; + const PetscScalar *Cxi, *Cyi, *Czi, *CTxi, *CTyi, *CTzi; + const PetscScalar *EntResi, *SuppSizei, *MCi; double solni; Tensor<1,dim> vi,vj; @@ -873,7 +873,7 @@ void LevelSetSolver::assemble_K_DL_DH_times_vector for (;idofs_iter!=locally_owned_dofs_LS.end(); idofs_iter++) { - int gi = *idofs_iter; + PetscInt gi = *idofs_iter; //double ith_K_times_solution = 0; // read velocity of i-th DOF @@ -897,7 +897,7 @@ void LevelSetSolver::assemble_K_DL_DH_times_vector MatGetRow(MC_matrix,gi,&ncolumns,&gj,&MCi); // get vector values for column indices - const std::vector gj_indices (gj,gj+ncolumns); + const std::vector gj_indices (gj,gj+ncolumns); std::vector soln(ncolumns); std::vector vx(ncolumns); std::vector vy(ncolumns); @@ -1115,7 +1115,7 @@ void LevelSetSolver::compute_bounds(PETScWrappers::MPI::Vector &un_solution int gi = *idofs_iter; // get solution at DOFs on the sparsity pattern of i-th DOF - std::vector gj_indices = sparsity_pattern[gi]; + std::vector gj_indices = sparsity_pattern[gi]; std::vector soln(gj_indices.size()); get_vector_values(un_solution,gj_indices,soln); // compute bounds, ith row of flux matrix, P vectors @@ -1138,7 +1138,7 @@ void LevelSetSolver::check_max_principle(PETScWrappers::MPI::Vector &unp1_s { // compute min and max vectors const unsigned int dofs_per_cell = fe_LS.dofs_per_cell; - std::vector local_dof_indices (dofs_per_cell); + std::vector local_dof_indices (dofs_per_cell); double tol=1e-10; typename DoFHandler::active_cell_iterator @@ -1204,9 +1204,9 @@ void LevelSetSolver::compute_MPP_uH // loop on locally owned i-DOFs (rows) IndexSet::ElementIterator idofs_iter = locally_owned_dofs_LS.begin(); - int ncolumns; - const int *gj; - const double *MCi, *dLi, *dCi; + PetscInt ncolumns; + const PetscInt *gj; + const PetscScalar *MCi, *dLi, *dCi; double solni, mi, solLi, solHi; for (;idofs_iter!=locally_owned_dofs_LS.end(); idofs_iter++) @@ -1224,7 +1224,7 @@ void LevelSetSolver::compute_MPP_uH MatGetRow(dCij_matrix,gi,&ncolumns,&gj,&dCi); // get vector values for support of i-th DOF - const std::vector gj_indices (gj,gj+ncolumns); + const std::vector gj_indices (gj,gj+ncolumns); std::vector soln(ncolumns); std::vector solH(ncolumns); get_vector_values(solution,gj_indices,soln); @@ -1288,7 +1288,7 @@ void LevelSetSolver::compute_MPP_uH MatGetRow(A_matrix,gi,&ncolumns,&gj,&Ai); // get vector values for column indices - const std::vector gj_indices (gj,gj+ncolumns); + const std::vector gj_indices (gj,gj+ncolumns); std::vector Rpos(ncolumns); std::vector Rneg(ncolumns); get_vector_values(R_pos_vector,gj_indices,Rpos); @@ -1327,9 +1327,9 @@ void LevelSetSolver::compute_MPP_uH_with_iterated_FCT LxAkp1_matrix.copy_from(LxA_matrix); // loop in num of FCT iterations - int ncolumns; - const int *gj; - const double *Akp1i; + PetscInt ncolumns; + const PetscInt *gj; + const PetscScalar *Akp1i; double mi; for (int iter=0; iter::compute_MPP_uH_with_iterated_FCT // get i-th row of matrices MatGetRow(Akp1_matrix,gi,&ncolumns,&gj,&Akp1i); // get vector values for support of i-th DOF - const std::vector gj_indices (gj,gj+ncolumns); + const std::vector gj_indices (gj,gj+ncolumns); std::vector soln(ncolumns); get_vector_values(un_solution,gj_indices,soln); @@ -1397,7 +1397,7 @@ void LevelSetSolver::compute_MPP_uH_with_iterated_FCT MatGetRow(Akp1_matrix,gi,&ncolumns,&gj,&Akp1i); // get vector values for column indices - const std::vector gj_indices(gj,gj+ncolumns); + const std::vector gj_indices(gj,gj+ncolumns); std::vector Rpos(ncolumns); std::vector Rneg(ncolumns); get_vector_values(R_pos_vector,gj_indices,Rpos); @@ -1489,13 +1489,13 @@ void LevelSetSolver::get_sparsity_pattern() { // loop on DOFs IndexSet::ElementIterator idofs_iter = locally_owned_dofs_LS.begin(); - int ncolumns; - const int *gj; - const double *MCi; + PetscInt ncolumns; + const PetscInt *gj; + const PetscScalar *MCi; for (;idofs_iter!=locally_owned_dofs_LS.end(); idofs_iter++) { - int gi = *idofs_iter; + PetscInt gi = *idofs_iter; // get i-th row of mass matrix (dummy, I just need the indices gj) MatGetRow(MC_matrix,gi,&ncolumns,&gj,&MCi); sparsity_pattern[gi] = std::vector(gj,gj+ncolumns); @@ -1508,9 +1508,9 @@ void LevelSetSolver::get_map_from_Q1_to_Q2() { map_from_Q1_to_Q2.clear(); const unsigned int dofs_per_cell_LS = fe_LS.dofs_per_cell; - std::vector local_dof_indices_LS (dofs_per_cell_LS); + std::vector local_dof_indices_LS (dofs_per_cell_LS); const unsigned int dofs_per_cell_U = fe_U.dofs_per_cell; - std::vector local_dof_indices_U (dofs_per_cell_U); + std::vector local_dof_indices_U (dofs_per_cell_U); typename DoFHandler::active_cell_iterator cell_LS = dof_handler_LS.begin_active(), @@ -1565,7 +1565,7 @@ void LevelSetSolver::save_old_vel_solution() // ------------------------------------------------------------------------------- // template void LevelSetSolver::get_vector_values (PETScWrappers::VectorBase &vector, - const std::vector &indices, + const std::vector &indices, std::vector &values) { // PETSc wrapper to get sets of values from a petsc vector. @@ -1581,7 +1581,7 @@ void LevelSetSolver::get_vector_values (PETScWrappers::VectorBase &vector, IndexSet ghost_indices = locally_relevant_dofs_LS; ghost_indices.subtract_set(locally_owned_dofs_LS); - int n_idx, begin, end, i; + PetscInt n_idx, begin, end, i; n_idx = indices.size(); VecGetOwnershipRange (vector, &begin, &end); @@ -1609,7 +1609,7 @@ void LevelSetSolver::get_vector_values (PETScWrappers::VectorBase &vector, template void LevelSetSolver::get_vector_values (PETScWrappers::VectorBase &vector, - const std::vector &indices, + const std::vector &indices, std::map &map_from_Q1_to_Q2, std::vector &values) { @@ -1627,7 +1627,7 @@ void LevelSetSolver::get_vector_values (PETScWrappers::VectorBase &vector, IndexSet ghost_indices = locally_relevant_dofs_U; ghost_indices.subtract_set(locally_owned_dofs_U); - int n_idx, begin, end, i; + PetscInt n_idx, begin, end, i; n_idx = indices.size(); VecGetOwnershipRange (vector, &begin, &end); diff --git a/two_phase_flow/MultiPhase.cc b/two_phase_flow/MultiPhase.cc index a64e800..dd525c6 100644 --- a/two_phase_flow/MultiPhase.cc +++ b/two_phase_flow/MultiPhase.cc @@ -67,7 +67,7 @@ public: private: void set_boundary_inlet(); void get_boundary_values_U(); - void get_boundary_values_phi(std::vector &boundary_values_id_phi, + void get_boundary_values_phi(std::vector &boundary_values_id_phi, std::vector &boundary_values_phi); void output_results(); void output_vectors(); @@ -109,9 +109,9 @@ private: PETScWrappers::MPI::Vector completely_distributed_solution_v; PETScWrappers::MPI::Vector completely_distributed_solution_p; // BOUNDARY VECTORS - std::vector boundary_values_id_u; - std::vector boundary_values_id_v; - std::vector boundary_values_id_phi; + std::vector boundary_values_id_u; + std::vector boundary_values_id_v; + std::vector boundary_values_id_phi; std::vector boundary_values_u; std::vector boundary_values_v; std::vector boundary_values_phi; @@ -268,9 +268,9 @@ void MultiPhase::init_constraints() template void MultiPhase::get_boundary_values_U() { - std::map map_boundary_values_u; - std::map map_boundary_values_v; - std::map map_boundary_values_w; + std::map map_boundary_values_u; + std::map map_boundary_values_v; + std::map map_boundary_values_w; // NO-SLIP CONDITION if (PROBLEM==BREAKING_DAM || PROBLEM==FALLING_DROP) @@ -325,8 +325,8 @@ void MultiPhase::get_boundary_values_U() boundary_values_id_v.resize(map_boundary_values_v.size()); boundary_values_u.resize(map_boundary_values_u.size()); boundary_values_v.resize(map_boundary_values_v.size()); - std::map::const_iterator boundary_value_u =map_boundary_values_u.begin(); - std::map::const_iterator boundary_value_v =map_boundary_values_v.begin(); + std::map::const_iterator boundary_value_u =map_boundary_values_u.begin(); + std::map::const_iterator boundary_value_v =map_boundary_values_v.begin(); for (int i=0; boundary_value_u !=map_boundary_values_u.end(); ++boundary_value_u, ++i) { @@ -372,10 +372,10 @@ void MultiPhase::set_boundary_inlet() } template -void MultiPhase::get_boundary_values_phi(std::vector &boundary_values_id_phi, +void MultiPhase::get_boundary_values_phi(std::vector &boundary_values_id_phi, std::vector &boundary_values_phi) { - std::map map_boundary_values_phi; + std::map map_boundary_values_phi; unsigned int boundary_id=0; set_boundary_inlet(); @@ -383,7 +383,7 @@ void MultiPhase::get_boundary_values_phi(std::vector &boundar VectorTools::interpolate_boundary_values (dof_handler_LS,boundary_id,BoundaryPhi(1.0),map_boundary_values_phi); boundary_values_id_phi.resize(map_boundary_values_phi.size()); boundary_values_phi.resize(map_boundary_values_phi.size()); - std::map::const_iterator boundary_value_phi = map_boundary_values_phi.begin(); + std::map::const_iterator boundary_value_phi = map_boundary_values_phi.begin(); for (int i=0; boundary_value_phi !=map_boundary_values_phi.end(); ++boundary_value_phi, ++i) { boundary_values_id_phi[i]=boundary_value_phi->first; diff --git a/two_phase_flow/NavierStokesSolver.cc b/two_phase_flow/NavierStokesSolver.cc index d7d5c4c..761d28e 100644 --- a/two_phase_flow/NavierStokesSolver.cc +++ b/two_phase_flow/NavierStokesSolver.cc @@ -84,12 +84,12 @@ public: PETScWrappers::MPI::Vector locally_relevant_solution_w, PETScWrappers::MPI::Vector locally_relevant_solution_p); //boundary conditions - void set_boundary_conditions(std::vector boundary_values_id_u, - std::vector boundary_values_id_v, std::vector boundary_values_u, + void set_boundary_conditions(std::vector boundary_values_id_u, + std::vector boundary_values_id_v, std::vector boundary_values_u, std::vector boundary_values_v); - void set_boundary_conditions(std::vector boundary_values_id_u, - std::vector boundary_values_id_v, - std::vector boundary_values_id_w, std::vector boundary_values_u, + void set_boundary_conditions(std::vector boundary_values_id_u, + std::vector boundary_values_id_v, + std::vector boundary_values_id_w, std::vector boundary_values_u, std::vector boundary_values_v, std::vector boundary_values_w); void set_velocity(PETScWrappers::MPI::Vector locally_relevant_solution_u, PETScWrappers::MPI::Vector locally_relevant_solution_v); @@ -184,9 +184,9 @@ private: ConstraintMatrix constraints; ConstraintMatrix constraints_psi; - std::vector boundary_values_id_u; - std::vector boundary_values_id_v; - std::vector boundary_values_id_w; + std::vector boundary_values_id_u; + std::vector boundary_values_id_v; + std::vector boundary_values_id_w; std::vector boundary_values_u; std::vector boundary_values_v; std::vector boundary_values_w; @@ -357,8 +357,8 @@ void NavierStokesSolver::initial_condition(PETScWrappers::MPI::Vector local } template -void NavierStokesSolver::set_boundary_conditions(std::vector boundary_values_id_u, - std::vector boundary_values_id_v, +void NavierStokesSolver::set_boundary_conditions(std::vector boundary_values_id_u, + std::vector boundary_values_id_v, std::vector boundary_values_u, std::vector boundary_values_v) { @@ -369,9 +369,9 @@ void NavierStokesSolver::set_boundary_conditions(std::vector } template -void NavierStokesSolver::set_boundary_conditions(std::vector boundary_values_id_u, - std::vector boundary_values_id_v, - std::vector boundary_values_id_w, +void NavierStokesSolver::set_boundary_conditions(std::vector boundary_values_id_u, + std::vector boundary_values_id_v, + std::vector boundary_values_id_w, std::vector boundary_values_u, std::vector boundary_values_v, std::vector boundary_values_w) @@ -641,7 +641,7 @@ void NavierStokesSolver::assemble_system_U() std::vector > grad_psiqn(n_q_points); std::vector > grad_psiqnm1(n_q_points); - std::vector local_dof_indices(dofs_per_cell); + std::vector local_dof_indices(dofs_per_cell); std::vector > shape_grad(dofs_per_cell); std::vector shape_value(dofs_per_cell); @@ -868,7 +868,7 @@ void NavierStokesSolver::assemble_system_dpsi_q() { std::vector > gvnp1(n_q_points); std::vector > gwnp1(n_q_points); - std::vector local_dof_indices(dofs_per_cell); + std::vector local_dof_indices(dofs_per_cell); std::vector shape_value(dofs_per_cell); std::vector > shape_grad(dofs_per_cell); @@ -913,7 +913,7 @@ void NavierStokesSolver::assemble_system_dpsi_q() { for (unsigned int i=0; i::assemble_system_dpsi_q() { cell_S(i,j)+=shape_grad[i]*shape_grad[j]*JxW; cell_M(i,j)+=shape_value[i]*shape_value[j]*JxW; } + } } } cell_P->get_dof_indices(local_dof_indices); -- 2.39.5