try damping in leaf elements for shock

JuntaoHuang · Feb 14, 2024 · 705e53c · 705e53c
1 parent c2c1ab4
commit 705e53c
Show file tree

Hide file tree

Showing 4 changed files with 82 additions and 0 deletions.
diff --git a/include/DGAdapt.h b/include/DGAdapt.h
@@ -67,6 +67,9 @@ class DGAdapt :
 	// refine based on reference solution generated by Euler forward
 	void refine();
 
+	// add damping for the leaf elements
+	void damping_leaf(const double damp_coef);
+
 	// refine to some max mesh level in given dimension
 	// NOT TESTED
 	void refine(const int max_mesh, const std::vector<int> dim);

diff --git a/include/DGSolution.h b/include/DGSolution.h
@@ -249,6 +249,8 @@ friend class IO;
 	 */
 	void source_elem_separable(Element & elem, const VecMultiD<double> & coefficient_1D, const int index_var = 0);
 
+	void print_order_all_basis_in_dg();
+
 protected:
 
 	const bool sparse;		///< control sparse grid (true) or full grid (false)

diff --git a/source/DGAdapt.cpp b/source/DGAdapt.cpp
@@ -403,6 +403,61 @@ void DGAdapt::refine()
     update_order_all_basis_in_dgmap();
 }
 
+void DGAdapt::damping_leaf(const double damp_coef)
+{
+	for (auto & iter : this->dg)
+	{	
+		if ((iter.second.level[0] == NMAX))
+		{
+			iter.second.ucoe_alpt[0] *= damp_coef;
+		}
+	}
+
+	for (auto & i : viscosity_element)	
+	{
+		// auto iter = *i;
+		// print(iter.level[0]);
+		// iter.ucoe_alpt[0] *= damp_coef;
+		// for (auto & iter_par : iter.hash_ptr_par)
+		// {
+		// 	iter_par.second->ucoe_alpt[0] *= damp_coef;
+		// }
+
+		// if ((iter.second.level[0] == NMAX))
+		// {
+		// 	iter.second.ucoe_alpt[0] *= damp_coef;
+		// }
+        // // max mesh level of this element in all dimensions
+        // const int max_mesh_level = *(std::max_element(std::begin(iter.second->level), std::end(iter.second->level)));
+
+        // // only consider element in finest mesh level
+        // if (max_mesh_level == NMAX)
+		// {	
+		// 	// iter.second->ucoe_alpt[0].print();
+		// 	iter.second->ucoe_alpt[0] *= damp_coef;
+
+		// 	// for (auto & iter_par : iter.second->hash_ptr_par)
+		// 	// {
+		// 	// 	iter_par.second->ucoe_alpt[0] *= damp_coef;
+		// 	// }
+		// }
+
+		// // loop over all its parent elements
+		// const std::set<std::array<std::vector<int>,2>> index_chd_elem = index_all_par(iter.second->level, iter.second->suppt);
+		// for (auto const & index : index_chd_elem)
+		// {
+		// 	int hash_key = hash.hash_key(index);
+
+		// 	// if index is not in current child index set, then add it to dg solution
+		// 	if ( iter.second->hash_ptr_chd.find(hash_key) == iter.second->hash_ptr_chd.end())
+		// 	{
+		// 		Element elem(index[0], index[1], all_bas, hash);
+		// 		add_elem(elem);
+		// 	}                
+		// }			
+	}
+}
+
 // NOT TESTED
 void DGAdapt::refine(const int max_mesh, const std::vector<int> dim)
 {

diff --git a/source/DGSolution.cpp b/source/DGSolution.cpp
@@ -906,7 +906,29 @@ void DGSolution::update_order_all_basis_in_dgmap()
 	}
 }
 
+void DGSolution::print_order_all_basis_in_dg()
+{
+	int order_alpt_basis_in_dgmap = 0;
 
+	for (auto & iter : dg)
+	{
+		for (size_t num_vec = 0; num_vec < ind_var_vec.size(); num_vec++)
+		{
+			// update order of alpert basis
+			for (size_t num_basis = 0; num_basis < iter.second.size_alpt(); num_basis++)
+			{
+				const std::vector<int> & order_local_basis = iter.second.order_local_alpt[num_basis];
+				iter.second.order_alpt_basis_in_dg[num_vec].at(order_local_basis) = order_alpt_basis_in_dgmap;
+
+				std::cout << "order: " << order_alpt_basis_in_dgmap
+						<< " level: " << iter.second.level[0]
+						<< " suppt: " << iter.second.suppt[0] <<  std::endl;
+
+				order_alpt_basis_in_dgmap++;
+			}
+		}
+	}
+}
 
 
 std::vector<int> DGSolution::max_mesh_level_vec() const