Add Examples

Add Examples/ADVDIFF example for advection-diffusion

Add exotic grid generators to Examples/GridTools

Add Examples/TestSuits for tests of diffusion and advection-diffusion
problems

Examples/ADFVDiscr/main.cpp

@@ -145,7 +145,7 @@ int main(int argc,char ** argv)
 					tag_BC[*face][1] = 0; //neumann
 					tag_BC[*face][2] = func(x,0);//face->Mean(func, 0);
 				}
-		}`
+		}
 		
 		if(m->HaveTag("REFERENCE_SOLUTION") )
 			phi_ref = m->GetTag("REFERENCE_SOLUTION");
@@ -304,44 +304,47 @@ int main(int argc,char ** argv)
 
 			ttt = Timer();
 
-			Tag error = m->CreateTag("error",DATA_REAL,CELL,NONE,1);
-
-			double err_C = 0.0, err_L2 = 0.0;
+			
+			if( phi_ref.isValid() )
+			{
+				Tag error = m->CreateTag("error",DATA_REAL,CELL,NONE,1);
+				double err_C = 0.0, err_L2 = 0.0;
 #if defined(USE_OMP)
 #pragma omp parallel
 #endif
-			{
-				double local_err_C = 0;
+				{
+					double local_err_C = 0;
 #if defined(USE_OMP)
 #pragma omp for reduction(+:err_L2)
 #endif
-				for( Storage::integer icell = 0; icell < m->CellLastLocalID(); ++icell )
-				{
-					Cell cell = Cell(m,ComposeCellHandle(icell));
-					if( cell->GetStatus() != Element::Ghost )
+					for( Storage::integer icell = 0; icell < m->CellLastLocalID(); ++icell )
 					{
-						double old = phi[cell];
-						double exact = phi_ref[cell];
-						double res = Update[Phi.Index(cell)];
-						double sol = old-res;
-						double err = fabs (sol - exact);
-						if (err > local_err_C) local_err_C = err;
-						err_L2 += err * err * cell->Volume();
-						cell->Real(error) = err;
-						phi[cell] = sol;
+						Cell cell = Cell(m,ComposeCellHandle(icell));
+						if( cell->GetStatus() != Element::Ghost )
+						{
+							double old = phi[cell];
+							double exact = phi_ref[cell];
+							double res = Update[Phi.Index(cell)];
+							double sol = old-res;
+							double err = fabs (sol - exact);
+							if (err > local_err_C) local_err_C = err;
+							err_L2 += err * err * cell->Volume();
+							cell->Real(error) = err;
+							phi[cell] = sol;
+						}
 					}
-				}
 #if defined(USE_OMP)
 #pragma omp critical
 #endif
-				{
-					if( local_err_C > err_C ) err_C = local_err_C;
+					{
+						if( local_err_C > err_C ) err_C = local_err_C;
+					}
 				}
+				err_C = m->AggregateMax(err_C); // Compute the maximal C norm for the error
+				err_L2 = sqrt(m->Integrate(err_L2)); // Compute the global L2 norm for the error
+				if( m->GetProcessorRank() == 0 ) std::cout << "err_C  = " << err_C << std::endl;
+				if( m->GetProcessorRank() == 0 ) std::cout << "err_L2 = " << err_L2 << std::endl;
 			}
-			err_C = m->AggregateMax(err_C); // Compute the maximal C norm for the error
-			err_L2 = sqrt(m->Integrate(err_L2)); // Compute the global L2 norm for the error
-			if( m->GetProcessorRank() == 0 ) std::cout << "err_C  = " << err_C << std::endl;
-			if( m->GetProcessorRank() == 0 ) std::cout << "err_L2 = " << err_L2 << std::endl;
 		}
 		BARRIER;
 		if( m->GetProcessorRank() == 0 ) std::cout << "Compute true residual: " << Timer()-ttt << std::endl;

Examples/ADMFD/main.cpp

@@ -165,7 +165,7 @@ int main(int argc,char ** argv)
 #endif
 			{
 				rMatrix NK, R, Areas;
-				rMatrix x(1,3), xf(1,3), n(1,3)
+				rMatrix x(1,3), xf(1,3), n(1,3);
 				double area; //area of the face
 				double volume; //volume of the cell
 				Areas.Zero();
@@ -193,14 +193,14 @@ int main(int argc,char ** argv)
 						 faces[k].Centroid(xf.data());
 						 faces[k].OrientedUnitNormal(cell->self(),n.data());
 						 // assemble matrix of directions
-						 R(k,k+1,0,3) = (yf-x)*area;
+						 R(k,k+1,0,3) = (xf-x)*area;
 						 // assemble matrix of co-normals
 						 NK(k,k+1,0,3) = n*K;
 						 Areas(k,k) = area;
 					 } //end of loop over faces
 					 W = NK*(NK.Transpose()*R).Invert(true).first*NK.Transpose(); //stability part
 					 W+=(rMatrix::Unit(NF) - R*(R.Transpose()*R).Invert(true).first*R.Transpose())*
-						(2.0/(static_cast<real>(NF)*vP)*(NK*K.Invert(true).first*NK.Transpose()).Trace());
+						(2.0/(static_cast<real>(NF)*volume)*(NK*K.Invert(true).first*NK.Transpose()).Trace());
 			 		 W = Areas*W*Areas;
 					 //access data structure for gradient matrix in mesh
 					 real_array store_W = cell->RealArrayDV(tag_W);

Examples/ADVDIFF/CMakeLists.txt

+project(ADVDIFF)
+set(SOURCE main.cpp stencil.cpp checks.cpp limited_average.cpp save_mesh.cpp conv_diff.cpp)
+set(HEADER stencil.h checks.h limited_average.h save_mesh.h conv_diff.h)
+
+add_executable(ADVDIFF ${SOURCE} ${HEADER})
+
+target_link_libraries(ADVDIFF inmost)
+
+if(USE_SOLVER)
+  if(USE_SOLVER_ANI)
+    message("linking ADVDIFF with ani3d and BLAS")
+    target_link_libraries(ADVDIFF ani3d ${BLAS_LIBRARIES})
+    if(BLAS_LINKER_FLAGS)
+      set_target_properties(ADVDIFF PROPERTIES LINK_FLAGS "${BLAS_LINKER_FLAGS}")
+    endif()
+  endif()
+  if(USE_SOLVER_PETSC)
+    message("linking ADVDIFF with PETSc")
+    target_link_libraries(ADVDIFF ${PETSC_LIBRARIES})
+  endif()
+  if(USE_SOLVER_TRILINOS)
+    message("linking ADVDIFF with Trilinos")
+    target_link_libraries(ADVDIFF ${Trilinos_LIBRARIES} ${Trilinos_TPL_LIBRARIES})
+  endif()
+  if(USE_SOLVER_METIS)
+    message("linking ADVDIFF with Metis")
+    target_link_libraries(ADVDIFF ${METIS_LIBRARIES})
+  endif()
+  if(USE_SOLVER_MONDRIAAN)
+    message("linking ADVDIFF with Mondriaan")
+    target_link_libraries(ADVDIFF ${MONDRIAAN_LIBRARIES})
+  endif()
+  if(USE_SOLVER_SUPERLU)
+    message("linking ADVDIFF with SuperLU")
+    target_link_libraries(ADVDIFF ${SUPERLU_LIBRARIES})
+  endif()
+endif()
+
+
+if(USE_PARTITIONER)
+  if(USE_PARTITIONER_ZOLTAN)
+    message("linking ADVDIFF with Zoltan")
+    target_link_libraries(ADVDIFF ${ZOLTAN_LIBRARIES})
+  endif()
+  if(USE_PARTITIONER_PARMETIS)
+    message("linking ADVDIFF with ParMETIS")
+    target_link_libraries(ADVDIFF ${PARMETIS_LIBRARIES})
+  endif()
+endif()
+
+if(USE_MPI)
+  message("linking ADVDIFF with MPI")
+  target_link_libraries(ADVDIFF ${MPI_LIBRARIES}) 
+  if(MPI_LINK_FLAGS)
+    set_target_properties(ADVDIFF PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
+  endif() 
+endif(USE_MPI)
+
+install(TARGETS ADVDIFF EXPORT inmost-targets RUNTIME DESTINATION bin)
\ No newline at end of file

Examples/ADVDIFF/checks.cpp

+#include "checks.h"
+
+using namespace INMOST;
+
+
+//shortcuts
+typedef Storage::bulk       bulk;
+typedef Storage::enumerator enumerator;
+typedef Storage::real       real;
+
+bool check_flux_properties(int i, const variable & flux, const std::string & name, bool print)
+{
+	bool test = true;
+	const Sparse::Row & r = flux.GetRow();
+	for(int k = 0; k < (int)r.Size(); ++k)
+	{
+		if( r.GetIndex(k) == i ) //diagonal element
+		{
+			if( r.GetValue(k) < 0.0 )
+			{
+				if( print ) std::cout << "flux " << name << " gives negative diagonal contribution to " << i << "-th variable" << std::endl;
+				test = false;
+			}
+		}
+		else //off-diagonal element
+		{
+			if( r.GetValue(k) > 0.0 )
+			{
+				if( print ) std::cout << "flux " << name << " gives positive off-diagonal contribution to " << i << "-th variable" << std::endl;
+				test = false;
+			}
+		}
+	}
+	if( !test && print ) 
+	{
+		r.Print();
+		//scanf("%*c");
+	}
+	return test;
+}
+
+
+bool check_matrix_properties(const Sparse::Matrix & A, bool print)
+{
+	bool is_m_matrix = true;
+	enumerator mbeg,mend;
+	A.GetInterval(mbeg,mend);
+	for(enumerator i = mbeg; i < mend; ++i)
+	{
+		const Sparse::Row & r = A[i];
+		real rowsum = 0.0;
+		for(enumerator k = 0; k < r.Size(); ++k)
+		{
+			rowsum += r.GetValue(k);
+			if( r.GetIndex(k) != i ) 
+			{
+				if( r.GetValue(k) > 0.0 )
+				{
+					if( print ) std::cout << "Off-diagonal element (" << i << "," << r.GetIndex(k) << ") is positive: " << r.GetValue(k) << std::endl;
+					is_m_matrix = false;
+				}
+			}
+			else 
+			{
+				if( r.GetValue(k) < 0.0 )
+				{
+					if( print ) std::cout << "Diagonal element (" << i << "," << r.GetIndex(k) << ") is negative: " << r.GetValue(k) << std::endl;
+					is_m_matrix = false;
+				}
+			}
+		}
+		if( rowsum < -1.0e-9 )
+		{
+			std::cout << "Row-sum is " << rowsum << std::endl;
+			is_m_matrix = false;
+		}
+	}
+	if( !is_m_matrix ) std::cout << "Not an M-matrix" << std::endl;
+	return is_m_matrix;
+}
\ No newline at end of file

Examples/ADVDIFF/checks.h

+#ifndef __CHECKS_H
+#define __CHECKS_H
+
+#include "inmost.h"
+
+/// Check that this will give M-matrix for i-th unknown.
+/// @param i Position of the diagonal unknown.
+/// @param flux Computed flux expression with derivatives.
+/// @param name Name of the flux to be printed out.
+/// @param print Print out errors.
+bool check_flux_properties(int i, const INMOST::variable & flux, const std::string & name, bool print = true);
+
+
+///Check monotonicity of the matrix.
+/// @param A Input matrix to be checked.
+/// @param print Print out errors.
+bool check_matrix_properties(const INMOST::Sparse::Matrix & A, bool print = true);
+
+
+#endif
\ No newline at end of file

Examples/ADVDIFF/conv_diff.h

+#ifndef __CONV_DIFF_H
+#define __CONV_DIFF_H
+
+#include "inmost.h"
+#include "limited_average.h"
+
+///Comment me
+class ConvectionDiffusion
+{
+	INMOST::Mesh * m;
+	//tags for convection
+	INMOST::Tag tag_CONV_CB;	// Coefficients for upwind corrector at back cell
+	INMOST::Tag tag_CONV_EB;	// Elements for upwind corrector at back cell
+	INMOST::Tag tag_CONV_RB;	// Right hand side for upwind corrector at back cell
+	INMOST::Tag tag_CONV_CF;	// Coefficients for upwind corrector at front cell
+	INMOST::Tag tag_CONV_EF;	// Elements for upwind corrector at front cell
+	INMOST::Tag tag_CONV_RF;	// Right hand side for upwind corrector at front cell
+	INMOST::Tag tag_CONV_CU;	// Upwind coefficient
+	INMOST::Tag tag_CONV_RU;	// Right hand side for upwind
+	INMOST::Tag tag_CONV_EU;	// Upstream cell
+	INMOST::Tag tag_CONV_F;		// A flag indicating is there
+								// a full nonlinear part  (internal     = 2),
+								// a one-sided correction (on outlet BC = 1),
+								// no correction          (on inlet BC  = 0)
+	INMOST::Tag tag_CONV_VU;	// Two vectors for upwind correctors at back and front cells
+	//tags for diffusion
+	INMOST::Tag tag_DIFF_CB;	// Coefficients for transversal correction at back cell
+	INMOST::Tag tag_DIFF_EB;	// Elements for transversal correction at back cell
+	INMOST::Tag tag_DIFF_RB;	// Right hand side for two transversal correction at back cell
+	INMOST::Tag tag_DIFF_CF;	// Coefficients for transversal correction at front cell
+	INMOST::Tag tag_DIFF_EF;	// Elements for transversal correction at front cell
+	INMOST::Tag tag_DIFF_RF;	// Right hand side for two transversal correction at back cell
+	INMOST::Tag tag_DIFF_CT;	// Two-point transmissibility
+	INMOST::Tag tag_DIFF_RT;	// Two-point right hand side
+	INMOST::Tag tag_DIFF_F;		// A flag indicating is there
+								// a full nonlinear part   (internal = 3),
+								// no nonlinear part       (internal = 2),
+								// a one-sided correction  (on BC    = 1),
+								// no one-sided correction (on BC    = 0)
+	INMOST::Tag tag_DIFF_VT;	// Vector for transversal correction to two-point part
+
+
+	INMOST::Tag tag_U;  // Normal velocity vector on faces
+	INMOST::Tag tag_K;  // Diffusion tensor
+	INMOST::Tag tag_BC; // Boundary conditions
+	
+
+	INMOST::MarkerType build_flux; // defines wheather to build flux approximation on interface in parallel
+	INMOST::MarkerType boundary_face; //defines boundary faces in parallel
+
+	bool initialization_failure; //there is a failure during intialization
+
+	bool perform_correction_diffusion; //add nonlinear correction to two-point flux
+	bool perform_correction_convection; //add nonlinear correction to single point upstream 
+public:
+	ConvectionDiffusion(INMOST::Mesh * _m, INMOST::Tag _tag_U, INMOST::Tag _tag_K, INMOST::Tag _tag_BC, INMOST::MarkerType _boundary_face, bool correct_convection, bool correct_diffusion);
+	bool Failure() const;
+	~ConvectionDiffusion();
+	void DiffusionFluxes(INMOST::abstract_variable & param, INMOST::Face fKL, INMOST::variable & fluxT, INMOST::variable & corrK, INMOST::variable & corrL, INMOST::variable & corrK_cK, INMOST::variable & corrL_cL) const;
+	void AdvectionFluxes(INMOST::abstract_variable & param, INMOST::Face fKL, INMOST::variable & fluxT, INMOST::variable & corrK, INMOST::variable & corrL, INMOST::variable & corrK_cK, INMOST::variable & corrL_cL) const;
+	void Averaging(SchemeType scheme_type, INMOST_DATA_REAL_TYPE regularization, const INMOST::variable & fluxT, const INMOST::variable & corrK, const INMOST::variable & corrL, const INMOST::variable & corrK_cK, const INMOST::variable & corrL_cL, INMOST::variable & fluxK, INMOST::variable & fluxL, bool boundary) const;
+	
+	bool BuildFlux(INMOST::Face f) const;
+};
+
+#endif //__CONV_DIFF_H
\ No newline at end of file

Examples/ADVDIFF/limited_average.cpp

+#include "limited_average.h"
+
+using namespace INMOST;
+
+//shortcuts
+typedef Storage::real       real;
+typedef Storage::integer    integer;
+typedef Storage::enumerator enumerator;
+
+/// Function that performs non-linear weighted convex combination of two fluxes in
+/// a style of non-linear two-point flux approximation.
+/// Non-linear two-point flux approximation has positivity-preserving property.
+/// Depending on the template type it will account or not for derivatives in
+/// convex combination, which may result in faster convergence but loss of 
+/// positivity-preserving property. 
+/// @param fluxK Flux for back cell.
+/// @param fluxL Flux for front cell.
+/// @param fluxK_cK Flux for back cell only part with back cell.
+/// @param fluxK_cL Flux for back cell only part with front cell.
+/// @param fluxL_cK Flux for front cell only part with back cell.
+/// @param fluxL_cL Flux for front cell only part with front cell.
+/// @param outK Flux for the back cell.
+/// @param outL FLux for the front cell. May be different from outK.
+/// @param regularization Regularization parameter in computation of expression.
+template<typename VarType> 
+void LimitedAverageTwoPoint(const INMOST::variable & fluxK,
+	                        const INMOST::variable & fluxL,
+							const INMOST::variable & fluxK_cK,
+							const INMOST::variable & fluxL_cL,
+							INMOST::variable & outK, 
+							INMOST::variable & outL,
+							real regularization)
+{
+	VarType muK, muL, div; //convex combination that elemenates additional entries
+	assign(muL, -(fluxK - fluxK_cK));
+	assign(muK, +(fluxL - fluxL_cL));
+	// diffusion and advection fluxes come with a different sign,
+	// this trick should regularize both. Regularization is needed
+	// in order to avoid tottaly zero flux with zero derivatives.
+	if( muK < 0.0 || muL < 0.0 )
+	{
+		muK -= regularization;
+		muL -= regularization;
+	}
+	else
+	{
+		muK += regularization;
+		muL += regularization;
+	}
+	div = muK+muL;
+	muK /= div;
+	muL /= div;
+	// two-point transmissibilities for each cell
+	outK = outL = (fluxK_cK*muK) + (fluxL_cL*muL);
+}
+
+/// Function performs non-linear weighted convex combination for two fluxes with Van-Albada limited averages.
+/// Van-Albada limited averages function is supposed to enforce TVD property. 
+/// It does not lead to always positive approximation which is required to satisfy discrete maximum principle on each iteration.
+/// This funciton returns single flux definition, thus it is conservative on each iteration.
+/// Depending on the template type it will account or not for derivatives in convex combination.
+/// @param inK Flux at one side of the face.
+/// @param inL Flux at another side of the face.
+/// @param outK Limited conservative flux for one side of the face.
+/// @param outL Limited conservative flux for another side of the face.
+/// @param regularization Regularization parameter in computation of expression.
+template<typename VarType> 
+void LimitedAverageVanAlbada(const variable & inK, const variable & inL, variable & outK, variable & outL, real regularization)
+{
+	VarType vK, vL; //values on each side
+	VarType muK, muL; //multipliers for each side
+	VarType denomenator;
+	assign(vK,inK);
+	assign(vL,inL);
+	muK = vL*vL + regularization;
+	muL = vK*vK + regularization;
+	//compute nonlinear coefficients
+	denomenator = muK + muL;
+	//multipliers for each side
+	muK /= denomenator;
+	muL /= denomenator;
+	//output
+	outK = outL = muK*inK + muL*inL;
+}
+
+
+/// Function performs non-linear weighted convex combination for two fluxes with coefficients in power of 4.
+/// This limited averages function is supposed to enforce TVD property. 
+/// It does not lead to always positive approximation which is required to satisfy discrete maximum principle on each iteration.
+/// This funciton returns single flux definition, thus it is conservative on each iteration.
+/// Depending on the template type it will account or not for derivatives in convex combination.
+/// @param inK Flux at one side of the face.
+/// @param inL Flux at another side of the face.
+/// @param outK Limited conservative flux for one side of the face.
+/// @param outL Limited conservative flux for another side of the face.
+/// @param regularization Regularization parameter in computation of expression.
+template<typename VarType> 
+void LimitedAverageQuadratic(const variable & inK, const variable & inL, variable & outK, variable & outL, real regularization)
+{
+	VarType vK, vL; //values on each side
+	VarType muK, muL; //multipliers for each side
+	VarType denomenator;
+	assign(vK,inK);
+	assign(vL,inL);
+	muK = vL*vL*vL*vL + regularization;
+	muL = vK*vK*vK*vK + regularization;
+	//compute nonlinear coefficients
+	denomenator = muK + muL;
+	//multipliers for each side
+	muK /= denomenator;
+	muL /= denomenator;
+	//output
+	outK = outL = muK*inK + muL*inL;
+}
+
+/// Function performs non-linear weighted convex combination for two fluxes with Van-Albada limited averages.
+/// Van-Albada limited averages function is supposed to enforce TVD property. 
+/// It does not lead to always positive approximation which is required to satisfy discrete maximum principle on each iteration.
+/// This funciton returns dual flux definition, thus it is not conservative on each iteration.
+/// No large difference from LimitedAverageVanAlbada if we account for all derivative.
+/// @param inK Flux at one side of the face.
+/// @param inL Flux at another side of the face.
+/// @param outK Limited conservative flux for one side of the face.
+/// @param outL Limited conservative flux for another side of the face.
+/// @param regularization Regularization parameter in computation of expression.
+void LimitedAverageVanAlbadaDual(const variable & inK, const variable & inL, variable & outK, variable & outL, real regularization)
+{
+	real vK, vL; //values on each side
+	real cK, cL; //multipliers for each side
+	real denomenator;
+	assign(vK,inK);
+	assign(vL,inL);
+	//compute nonlinear coefficients
+	denomenator = vK*vK + vL*vL + 2*regularization;
+	//multipliers for each side
+	cK = (vL*(vK + vL) + regularization*2)/denomenator;
+	cL = (vK*(vK + vL) + regularization*2)/denomenator;
+	//output
+	outK = cK*inK;
+	outL = cL*inL;
+}
+
+/// Function performs non-linear weighted convex combination for two fluxes with Van-Albada limited averages with correction.
+/// Corrected Van-Albada limited averages function is supposed to enforce discrete maximum principle property on convergence. 
+/// This funciton returns single flux definition, thus it is conservative on each iteration.
+/// Depending on the template type it will account or not for derivatives in convex combination.
+/// @param inK Flux at one side of the face.
+/// @param inL Flux at another side of the face.
+/// @param outK Limited conservative flux for one side of the face.
+/// @param outL Limited conservative flux for another side of the face.
+/// @param regularization Regularization parameter in computation of expression.
+template<typename VarType> 
+void LimitedAverageVanAlbadaCorrected(const variable & inK, const variable & inL, variable & outK, variable & outL, real regularization)
+{
+	VarType vK, vL; //values on each side
+	VarType muK, muL; //multipliers for each side
+	VarType denomenator;
+	VarType sign; //sign of the product of fluxes
+	VarType corr; //correction
+	assign(vK,inK);
+	assign(vL,inL);
+	sign = soft_sign(vK*vL,regularization);
+	//sign = 2.0/(1+exp(-vK*vL))-1.0;
+	corr = 0.5*(1.0-sign);
+	muK = vL*vL - corr*vK*vL + regularization;
+	muL = vK*vK - corr*vK*vL + regularization;
+	//compute nonlinear coefficients
+	denomenator = muK + muL;
+	//multipliers for each side
+	muK /= denomenator;
+	muL /= denomenator;
+	//output
+	outK = outL = muK*inK + muL*inL;
+}
+
+
+/// Function performs non-linear weighted convex combination for two fluxes with coefficients in power of 4and correction.
+/// Corrected quadratic limited averages function is supposed to enforce discrete maximum principle property on convergence. 
+/// This funciton returns single flux definition, thus it is conservative on each iteration.
+/// Depending on the template type it will account or not for derivatives in convex combination.
+/// @param inK Flux at one side of the face.
+/// @param inL Flux at another side of the face.
+/// @param outK Limited conservative flux for one side of the face.
+/// @param outL Limited conservative flux for another side of the face.
+/// @param regularization Regularization parameter in computation of expression.
+template<typename VarType> 
+void LimitedAverageQuadraticCorrected(const variable & inK, const variable & inL, variable & outK, variable & outL, real regularization)
+{
+	VarType vK, vL; //values on each side
+	VarType muK, muL; //multipliers for each side
+	VarType denomenator;
+	VarType sign; //sign of the product of fluxes
+	VarType corr; //correction
+	assign(vK,inK);
+	assign(vL,inL);
+	sign = soft_sign(vK*vL,regularization);
+	//sign = 2.0/(1+exp(-vK*vL))-1.0;
+	corr = 0.5*(1.0-sign);
+	muK = vL*vL*vL*vL - corr*vL*vK*vK*vK + regularization;
+	muL = vK*vK*vK*vK - corr*vK*vL*vL*vL + regularization;
+	//compute nonlinear coefficients
+	denomenator = muK + muL;
+	//multipliers for each side
+	muK /= denomenator;
+	muL /= denomenator;
+	//output
+	outK = outL = muK*inK + muL*inL;
+}
+
+/// Function performs non-linear weighted convex combination for two fluxes with Van-Albada limited averages 
+/// with correction and dual flux definition.
+/// Not conservative on each iteration but should satisfy discrete maximum principle on each iteration (up to regularization).
+/// @param inL Flux at another side of the face.
+/// @param outK Limited conservative flux for one side of the face.
+/// @param outL Limited conservative flux for another side of the face.
+/// @param regularization Regularization parameter in computation of expression.
+void LimitedAverageVanAlbadaCorrectedDual(const variable & inK, const variable & inL, variable & outK, variable & outL, real regularization)
+{
+	real vK, vL; //values on each side
+	real cK, cL; //multipliers for each side
+	real sign, corr, denomenator;
+	vK = get_value(inK);
+	vL = get_value(inL);
+	//positivity corrector
+	sign = soft_sign(vK*vL,regularization);
+	corr = 0.5*(1.0 - sign);
+	//compute nonlinear coefficients
+	denomenator = vK*vK + vL*vL - 2*corr*vK*vL + 2*regularization;
+	//multipliers for each side
+	cK = (vL*(vK + vL)*(1-corr) + regularization*2)/denomenator;
+	cL = (vK*(vK + vL)*(1-corr) + regularization*2)/denomenator;
+	//output
+	outK = cK*inK;
+	outL = cL*inL;
+}
+
+
+/// Function performs non-linear weighted convex combination for two fluxes with Van-Leer limited averages.
+/// Van-Leer limited averages function is supposed to enforce discrete maximum principle property. 
+/// This funciton returns single flux definition, thus it is conservative on each iteration.
+/// Depending on the template type it will account or not for derivatives in
+/// convex combination.
+/// This method will satsify discrete maximum principle only on convergence, but not on each
+/// nonlinear iteration. To satisfy discrete maximum principle on each iteration see dual flux 
+/// Van-Leer limited averages function.
+/// @param inK Flux at one side of the face.
+/// @param inL Flux at another side of the face.
+/// @param outK Limited conservative flux for one side of the face.
+/// @param outL Limited conservative flux for another side of the face.
+/// @param regularization Regularization parameter in computation of expression.
+template<typename VarType> 
+void LimitedAverageVanLeer(const variable & inK, const variable & inL, variable & outK, variable & outL, real regularization)
+{