From 76eb1972fbe91fe428cd028a896f628bdef16c41 Mon Sep 17 00:00:00 2001 From: Roger Fu Date: Thu, 23 Feb 2017 15:01:07 -0500 Subject: [PATCH] updated ceres.cc to include output of core-mantle boundary --- CeresFE/src/ceres.cc | 194 ++++++++++++++++++++++++++----------------- 1 file changed, 118 insertions(+), 76 deletions(-) diff --git a/CeresFE/src/ceres.cc b/CeresFE/src/ceres.cc index 3c36cbd..945aac6 100644 --- a/CeresFE/src/ceres.cc +++ b/CeresFE/src/ceres.cc @@ -170,6 +170,7 @@ private: DoFHandler dof_handler; unsigned int n_u = 0, n_p = 0; unsigned int plastic_iteration = 0; + unsigned int last_max_plasticity = 0; QGauss quadrature_formula; std::vector< std::vector > > quad_viscosities; // Indices for this object are [cell][q][q coords, eta] @@ -684,6 +685,14 @@ void StokesProblem::assemble_system() { Assert( local_quadrature_points_history < &quadrature_point_history.back(), ExcInternalError()); + + double cell_area = cell->measure(); + if(cell_area<0) + append_physical_times(-1); + AssertThrow(cell_area > 0 + , + ExcInternalError()); + unsigned int m_id = cell->material_id(); @@ -726,12 +735,18 @@ void StokesProblem::assemble_system() { } for (unsigned int q = 0; q < n_q_points; ++q) { - const SymmetricTensor<2, dim> &old_stress = + SymmetricTensor<2, dim> &old_stress = local_quadrature_points_history[q].old_stress; double &local_old_phiphi_stress = local_quadrature_points_history[q].old_phiphi_stress; double r_value = fe_values.quadrature_point(q)[0]; double z_value = fe_values.quadrature_point(q)[1]; + + // if(system_parameters::present_timestep == system_parameters::initial_elastic_iterations) + // { + // old_stress *= 0; + // local_old_phiphi_stress = 0; + // } // get local density based on mat id double local_density = system_parameters::rho[m_id]; @@ -1099,17 +1114,28 @@ void StokesProblem::solution_stesses() { std::vector > vector_values(0); std::vector < std::vector > > gradient_values(0); std::vector failing_cells; + // Write the stresses from the previous step into vectors + std::vector> old_stress; + std::vector old_phiphi_stress; + std::vector cell_Gs; for (typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(); cell != dof_handler.end(); ++cell) { + // Makes pointer to data in quadrature_point_history + PointHistory *local_quadrature_points_history = + reinterpret_cast *>(cell->user_pointer()); + fe_values.reinit(cell); fe_values.get_function_gradients(solution, velocity_grads); fe_values.get_function_values(solution, velocities); Vector current_cell_velocity(dim+1); std::vector> current_cell_grads(dim+1); + SymmetricTensor<2, dim> current_cell_old_stress; + current_cell_old_stress = 0; + double current_cell_old_phiphi_stress = 0; double cell_area = 0; - // Averages across each cell to find mean velocities and gradients + // Averages across each cell to find mean velocities, gradients, and old stresses for (unsigned int q = 0; q < quadrature_formula.size(); ++q) { cell_area += fe_values.JxW(q); @@ -1120,13 +1146,24 @@ void StokesProblem::solution_stesses() { velocity_grads[q][i] *= fe_values.JxW(q); current_cell_grads[i] += velocity_grads[q][i]; } + current_cell_old_stress += local_quadrature_points_history[q].old_stress * fe_values.JxW(q); + current_cell_old_phiphi_stress += local_quadrature_points_history[q].old_phiphi_stress * fe_values.JxW(q); } current_cell_velocity /= cell_area; for (unsigned int i = 0; i < (dim+1); i++) current_cell_grads[i] /= cell_area; + current_cell_old_stress /= cell_area; + current_cell_old_phiphi_stress /= cell_area; vector_values.push_back(current_cell_velocity); gradient_values.push_back(current_cell_grads); + old_stress.push_back(current_cell_old_stress); + old_phiphi_stress.push_back(current_cell_old_phiphi_stress); + + // Get cell shear modulus: assumes it's constant for the cell + unsigned int mat_id = cell->material_id(); + double local_G = rheology.get_G(mat_id); + cell_Gs.push_back(local_G); } //tracks where failure occurred @@ -1152,6 +1189,17 @@ void StokesProblem::solution_stesses() { current_cell_viscosity = cell_viscosities[i]; } + + double cell_eta_ve = 2 + / ((1 / current_cell_viscosity) + + (1 / cell_Gs[i] + / system_parameters::current_time_interval)); + double cell_chi_ve = 1 + / (1 + + (cell_Gs[i] + * system_parameters::current_time_interval + / current_cell_viscosity)); + //find local pressure double cell_p = vector_values[i].operator()(2); //find stresses tensor @@ -1162,11 +1210,15 @@ void StokesProblem::solution_stesses() { A << gradient_values[i][0][0] << 0 << sigma13 << endr << 0 << vector_values[i].operator()(0) / points_list[i].operator()(0)<< 0 << endr << sigma13 << 0 << gradient_values[i][1][1] << endr; + mat olddevstress; + olddevstress << old_stress[i][0][0] << 0 << old_stress[i][0][1] << endr + << 0 << old_phiphi_stress[i] << 0 << endr + << old_stress[i][0][1] << 0 << old_stress[i][1][1] << endr; vec P; P << cell_p << cell_p << cell_p; mat Pmat = diagmat(P); mat B; - B = (2 * current_cell_viscosity * A) - Pmat; + B = (cell_eta_ve * A + cell_chi_ve * olddevstress) - Pmat; //finds principal stresses vec eigval; @@ -1204,7 +1256,7 @@ void StokesProblem::solution_stesses() { if (sigma3 < 0) temp_reductionfactor = 100; else - temp_reductionfactor = 1.5 * sigma1 / 5 / sigma3; + temp_reductionfactor = 1.9 * sigma1 / 5 / sigma3; reduction_factor.push_back(temp_reductionfactor); total_fails++; @@ -1275,7 +1327,13 @@ void StokesProblem::solution_stesses() { // If there are enough failed cells, update eta at all quadrature points and perform smoothing std::cout << " Number of failing cells: " << total_fails << "\n"; - if (total_fails <= 120) + double last_max_plasticity_double = last_max_plasticity; + double total_fails_double = total_fails; + double decrease_in_plasticity = ((last_max_plasticity_double - total_fails_double) / last_max_plasticity_double); + if(plastic_iteration == 0) + decrease_in_plasticity = 1; + last_max_plasticity = total_fails; + if (total_fails <= 100 || decrease_in_plasticity <= 0.2) { system_parameters::continue_plastic_iterations = false; for(unsigned int j=0; j < triangulation.n_active_cells(); j++) @@ -1529,10 +1587,14 @@ void StokesProblem::update_quadrature_point_history() { template void StokesProblem::update_time_interval() { - double move_goal_per_step = system_parameters::initial_disp_target - + double move_goal_per_step = system_parameters::initial_disp_target; + if(system_parameters::present_timestep > system_parameters::initial_elastic_iterations) + { + move_goal_per_step = system_parameters::initial_disp_target - ((system_parameters::initial_disp_target - system_parameters::final_disp_target) / - system_parameters::total_viscous_steps * - (system_parameters::present_timestep - system_parameters::initial_elastic_iterations)); + system_parameters::total_viscous_steps * + (system_parameters::present_timestep - system_parameters::initial_elastic_iterations)); + } double zero_tolerance = 1e-3; double max_velocity = 0; @@ -1574,11 +1636,12 @@ void StokesProblem::update_time_interval() } } } - std:: cout << "Breaker " << max_velocity; // NOTE: It is possible for this time interval to be very different from that used in the viscoelasticity calculation. system_parameters::current_time_interval = move_goal_per_step / max_velocity; double step_time_yr = system_parameters::current_time_interval / SECSINYEAR; - std::cout << "\n Viscous time for moving mesh: " << step_time_yr << " yr"; + std::cout << "Timestep interval changed to: " + << step_time_yr + << " years\n"; } //====================== MOVE MESH ====================== @@ -1598,7 +1661,6 @@ void StokesProblem::move_mesh() { for (unsigned int d = 0; d < dim; ++d) vertex_displacement[d] = solution( cell->vertex_dof_index(v, d)); - cell->vertex(v) += vertex_displacement * system_parameters::current_time_interval; } @@ -1675,8 +1737,11 @@ void StokesProblem::write_vertices(unsigned char boundary_that_we_need) { std::ofstream fout_final_vertices(vertices_output.str().c_str()); fout_final_vertices.close(); - // Figure out if the vertex is on the boundary of the domain std::vector vertex_touched(triangulation.n_vertices(), false); + + if (boundary_that_we_need == 0) + { + // Figure out if the vertex is on the boundary of the domain for (typename Triangulation::active_cell_iterator cell = triangulation.begin_active(); cell != triangulation.end(); ++cell) for (unsigned int f = 0; f < GeometryInfo::faces_per_cell; ++f) @@ -1694,6 +1759,37 @@ void StokesProblem::write_vertices(unsigned char boundary_that_we_need) { } } } + } + else + { + // Figure out if the vertex is on an internal boundary + for (typename Triangulation::active_cell_iterator cell = + triangulation.begin_active(); cell != triangulation.end(); ++cell) + for (unsigned int f = 0; f < GeometryInfo::faces_per_cell; ++f) + { + if (cell->neighbor(f) != triangulation.end()) { + if (cell->material_id() != cell->neighbor(f)->material_id()) //finds face is at internal boundary + { + int high_mat_id = std::max(cell->material_id(), + cell->neighbor(f)->material_id()); + if (high_mat_id == boundary_that_we_need) //finds faces at the correct internal boundary + { + for (unsigned int v = 0; + v < GeometryInfo::vertices_per_face; + ++v) + if (vertex_touched[cell->face(f)->vertex_index( + v)] == false) { + vertex_touched[cell->face(f)->vertex_index( + v)] = true; + std::ofstream fout_final_vertices(vertices_output.str().c_str(), std::ios::app); + fout_final_vertices << cell->face(f)->vertex(v) << "\n"; + fout_final_vertices.close(); + } + } + } + } + } + } } //====================== SETUP INITIAL MESH ====================== @@ -1777,32 +1873,6 @@ void StokesProblem::setup_initial_mesh() { triangulation.refine_global(system_parameters::global_refinement); -// -////refines region near r=0 -// if (system_parameters::small_r_refinement != 0) { -// for (unsigned int step = 0; -// step < system_parameters::small_r_refinement; ++step) { -// // std::cout << "iteration " << step + 1 << "\n"; -// typename dealii::Triangulation::active_cell_iterator cell = -// triangulation.begin_active(), endc = triangulation.end(); -// -// for (; cell != endc; ++cell) -// for (unsigned int v = 0; -// v < GeometryInfo::vertices_per_cell; ++v) { -// Point current_vertex = cell->vertex(v); -// -// const double x_coord = current_vertex.operator()(0); -// -// if (std::fabs(x_coord) < 1e-10) { -// cell->set_refine_flag(); -// break; -// } -// -// } -// triangulation.execute_coarsening_and_refinement(); -// } -// } -// //refines crustal region if (system_parameters::crustal_refinement != 0) { @@ -1836,30 +1906,6 @@ void StokesProblem::setup_initial_mesh() { } } -// -// //refines surface region -// if (system_parameters::surface_refinement != 0) { -// for (unsigned int step = 0; -// step < system_parameters::surface_refinement; ++step) { -// // std::cout << "iteration " << step + 1 << "\n"; -// typename dealii::Triangulation::active_cell_iterator cell = -// triangulation.begin_active(), endc = triangulation.end(); -// for (; cell != endc; ++cell) -// for (unsigned int f = 0; f < GeometryInfo::faces_per_cell; -// ++f) { -// if (cell->face(f)->at_boundary()) { -// if (cell->face(f)->center()[0] != 0) { -// if (cell->face(f)->center()[1] != 0) -// if (cell->face(f)->center()[0] > 1) { -// cell->set_refine_flag(); -// break; -// } -// } -// } -// } -// triangulation.execute_coarsening_and_refinement(); -// } -// } // output initial mesh in eps std::ostringstream refined_mesh_file; @@ -1927,9 +1973,10 @@ void StokesProblem::do_elastic_steps() if (system_parameters::present_timestep == 0) initialize_eta_and_G(); - - system_parameters::current_time_interval = - system_parameters::elastic_time; //This is the time interval needed in assembling the problem + + if(elastic_iteration == 0) + system_parameters::current_time_interval = + system_parameters::viscous_time; //This is the time interval needed in assembling the problem std::cout << " Assembling..." << std::endl << std::flush; assemble_system(); @@ -1940,15 +1987,14 @@ void StokesProblem::do_elastic_steps() output_results(); update_quadrature_point_history(); - // std::cout << std::endl << "\a"; + append_physical_times(0); elastic_iteration++; system_parameters::present_timestep++; - system_parameters::present_time = system_parameters::present_time + system_parameters::current_time_interval; - move_mesh(); do_ellipse_fits(); write_vertices(0); - write_vertices(1); + write_vertices(1); write_mesh(); + update_time_interval(); } } @@ -1992,13 +2038,10 @@ void StokesProblem::run() times_filename << system_parameters::output_folder << "/physical_times.txt"; std::ofstream fout_times(times_filename.str().c_str()); fout_times.close(); - append_physical_times(0); // Computes elastic timesteps do_elastic_steps(); - // Computes viscous timesteps - system_parameters::current_time_interval = system_parameters::viscous_time; unsigned int VEPstep = 0; while (system_parameters::present_timestep < (system_parameters::initial_elastic_iterations @@ -2010,18 +2053,17 @@ void StokesProblem::run() // Computes plasticity do_flow_step(); update_quadrature_point_history(); - update_time_interval(); + move_mesh(); append_physical_times(plastic_iteration); system_parameters::present_timestep++; system_parameters::present_time = system_parameters::present_time + system_parameters::current_time_interval; - move_mesh(); do_ellipse_fits(); write_vertices(0); write_vertices(1); write_mesh(); VEPstep++; } - append_physical_times(0); + append_physical_times(-1); } } @@ -2083,4 +2125,4 @@ int main(int argc, char* argv[]) { } return 0; -} \ No newline at end of file +} -- 2.39.5