Merge pull request #46 from nagelt/Minkley_Burgers_Fixes

Minkley burgers fixes

FEM/CMakeLists.txt

@@ -1,5 +1,6 @@
 set( HEADERS
 	BoundaryCondition.h
+	burgers.h
 	PhysicalConstant.h
 	Constrained.h
 	conversion_rate.h
@@ -16,9 +17,11 @@ set( HEADERS
 	files0.h
 	GeoInfo.h
 	InitialCondition.h
+	invariants.h
 	LinearFunctionData.h
 	mathlib.h
 	matrix_class.h
+	minkley.h
 	Output.h
 	pcs_dm.h
 	problem.h
@@ -54,6 +57,7 @@ set( HEADERS
 
 set( SOURCES
 	BoundaryCondition.cpp
+	burgers.cpp
 	CAP_IO.cpp
 	conversion_rate.cpp
 	DistributionInfo.cpp
@@ -72,9 +76,11 @@ set( SOURCES
 	files0.cpp
 	GeoInfo.cpp
 	InitialCondition.cpp
+	invariants.cpp
 	LinearFunctionData.cpp
 	mathlib.cpp
 	matrix_class.cpp
+	minkley.cpp
 	Output.cpp
 	pcs_dm.cpp
 	problem.cpp

FEM/burgers.cpp

@@ -0,0 +1,156 @@
+/**
+ * \copyright
+ * Copyright (c) 2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "burgers.h"
+#include "PhysicalConstant.h"
+namespace Burgers
+{
+SolidBurgers::SolidBurgers(const Math_Group::Matrix& data)
+{
+	GK0 = data(0); // Kelvin shear modulus
+	mK = data(1); // dependency parameter for "
+	etaK0 = data(2); // Kelvin viscosity
+	mvK = data(3); // dependency parameter for "
+	GM0 = data(4); // Maxwell shear modulus
+	KM0 = data(5); // Maxwell bulk modulus
+	etaM0 = data(6); // Maxwell viscosity
+	mvM = data(7); // dependency parameter for "
+	m_GM = data(8); // slope of elesticity temperature dependence
+	m_KM = data(9); // slope of elesticity temperature dependence
+	T_ref = data(10); // reference temperature dependency parameter for "
+	B = data(11); // constant factor for Arrhenius term
+	Q = data(12); // activation energy in Arrhenius term
+
+	GM = GM0;
+	KM = KM0;
+	GK = GK0;
+	etaK = etaK0;
+	etaM = etaM0;
+
+	if (!T_Process)
+	{
+		m_GM = 0.;
+		m_KM = 0.;
+		B = 1.;
+		Q = 0.; // for cutting off Arrhenius term
+		T_ref = 273.15;
+	}
+}
+
+/**************************************************************************
+   FEMLib-Method: Burgers::UpdateBurgersProperties()
+   Task: Updates BURGERS material parameters in LUBBY2 fashion
+   Programing:
+   07/2014 TN Implementation
+**************************************************************************/
+void SolidBurgers::UpdateBurgersProperties(double s_eff, const double Temperature)
+{
+	const double dT(Temperature - T_ref);
+	GM = GM0 + m_GM * dT;
+	KM = KM0 + m_KM * dT;
+
+	s_eff *= GM;
+
+	GK = GK0 * std::exp(mK * s_eff);
+	etaK = etaK0 * std::exp(mvK * s_eff);
+	etaM = etaM0 * std::exp(mvM * s_eff) * B
+	       * std::exp(Q * (-dT) / (PhysicalConstant::IdealGasConstant * Temperature * T_ref));
+	//	if (etaM / etaM0 < 1.e-2)
+	//		std::cout << "WARNING: Maxwell viscosity sank to 100th of original value." << std::endl;
+}
+
+/**************************************************************************
+   FEMLib-Method: Burgers::CalResidualBurgers()
+   Task: Calculates the 12x1 residual vector. Implementation fully implicit only.
+   Programing:
+   06/2014 TN Implementation
+**************************************************************************/
+void SolidBurgers::CalResidualBurgers(const double dt, const KVec& strain_curr, const KVec& stress_curr,
+                                      KVec& strain_Kel_curr, const KVec& strain_Kel_t, KVec& strain_Max_curr,
+                                      const KVec& strain_Max_t, Eigen::Matrix<double, 18, 1>& res)
+{
+	// calculate stress residual
+	res.block<6, 1>(0, 0) = stress_curr - 2. * (strain_curr - strain_Kel_curr - strain_Max_curr);
+
+	// calculate Kelvin strain residual
+	res.block<6, 1>(6, 0) = 1. / dt * (strain_Kel_curr - strain_Kel_t)
+	                        - 1. / (2. * etaK) * (GM * stress_curr - 2. * GK * strain_Kel_curr);
+
+	// calculate Maxwell strain residual
+	res.block<6, 1>(12, 0) = 1. / dt * (strain_Max_curr - strain_Max_t) - GM / (2. * etaM) * stress_curr;
+}
+
+/**************************************************************************
+   FEMLib-Method: Burgers::CalJacobianBurgers()
+   Task: Calculates the 12x12 Jacobian. Implementation fully implicit only.
+   Programing:
+   06/2014 TN Implementation
+**************************************************************************/
+void SolidBurgers::CalJacobianBurgers(const double dt, Eigen::Matrix<double, 18, 18>& Jac, const double s_eff,
+                                      const KVec& sig_i, const KVec& eps_K_i)
+{
+	// Check Dimension of Jacobian
+	// assert(Jac.cols() == 18 && Jac.rows() && 18);
+	Jac.setZero(18, 18);
+
+	// build G_11
+	Jac.block<6, 6>(0, 0) = SolidMath::ident;
+
+	// build G_12
+	Jac.block<6, 6>(0, 6) = 2. * SolidMath::ident;
+
+	// build G_13
+	Jac.block<6, 6>(0, 12) = 2. * SolidMath::ident;
+
+	// build G_21
+	Jac.block<6, 6>(6, 0) = -GM / (2. * etaK) * SolidMath::ident;
+
+	// build G_22
+	Jac.block<6, 6>(6, 6) = (1. / dt + GK / etaK) * SolidMath::ident;
+
+	// nothing to do for G_23
+
+	// build G_31
+	Jac.block<6, 6>(12, 0) = -GM / (2. * etaM) * SolidMath::ident;
+
+	// nothing to do for G_32
+
+	// build G_33
+	Jac.block<6, 6>(12, 12) = 1. / dt * SolidMath::ident;
+
+	if (s_eff > 0.)
+	{
+		const KVec dG_K = mK * 3. * GK * GM / (2. * s_eff) * sig_i;
+		const KVec dmu_vK = mvK * 3. * GM * etaK / (2. * s_eff) * sig_i;
+		const KVec dmu_vM = mvM * 3. * GM * etaM / (2. * s_eff) * sig_i;
+		const KVec eps_K_aid = 1. / (etaK * etaK) * (GM * sig_i - 2. * GK * eps_K_i);
+
+		// build G_21
+		Jac.block<6, 6>(6, 0) += 0.5 * eps_K_aid * dmu_vK.transpose() + 1. / etaK * eps_K_i * dG_K.transpose();
+
+		// build G_31
+		Jac.block<6, 6>(12, 0) += GM / (2. * etaM * etaM) * sig_i * dmu_vM.transpose();
+	}
+}
+
+/**************************************************************************
+   FEMLib-Method: Burgers::CaldGdEBurgers()
+   Task: Calculates the 12x6 derivative of the residuals with respect to total strain. Implementation fully implicit
+only.
+   Programing:
+   06/2014 TN Implementation
+**************************************************************************/
+void SolidBurgers::CaldGdEBurgers(Eigen::Matrix<double, 18, 6>& dGdE)
+{
+	// Check Dimension of dGdE
+	// assert(dGdE.cols() == 6 && dGdE.rows() == 18);
+	dGdE.setZero(18, 6);
+	dGdE.block<6, 6>(0, 0) = -2. * SolidMath::P_dev;
+}
+}

FEM/burgers.h

@@ -0,0 +1,44 @@
+/**
+ * \copyright
+ * Copyright (c) 2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef BURGERS_H
+#define BURGERS_H
+
+#include <cfloat>
+#include "Eigen/Dense"
+// FEM-Makros
+#include "makros.h"
+#include "tools.h"
+#include "invariants.h"
+
+namespace Burgers
+{
+typedef Eigen::Matrix<double, 6, 1> KVec;
+typedef Eigen::Matrix<double, 6, 6> KMat;
+
+class SolidBurgers
+{
+public:
+	explicit SolidBurgers(const Math_Group::Matrix& data);
+	~SolidBurgers() {}
+	// basic material parameters
+	double GK0, GM0, KM0, etaK0, etaM0, mK, mvK, mvM, l0, T_ref, m_GM, m_KM, B, Q;
+	// solution dependent values
+	double GM, KM, GK, etaK, etaM;
+
+	void UpdateBurgersProperties(double s_eff, const double Delta_T);
+	void CalResidualBurgers(const double dt, const KVec& strain_curr, const KVec& stress_curr, KVec& strain_Kel_curr,
+	                        const KVec& strain_Kel_t, KVec& strain_Max_curr, const KVec& strain_Max_t,
+	                        Eigen::Matrix<double, 18, 1>& res);
+	void CalJacobianBurgers(const double dt, Eigen::Matrix<double, 18, 18>& Jac, const double s_eff, const KVec& sig_i,
+	                        const KVec& eps_K_i);
+	void CaldGdEBurgers(Eigen::Matrix<double, 18, 6>& dGdE);
+};
+}
+#endif // BURGERS_H

FEM/equation_class.h

@@ -120,7 +120,7 @@ class Linear_EQS
 	double X(const long i) const { return x[i]; }
 	double RHS(const long i) const { return b[i]; }
 	double NormX();
-	double NormRHS() { return bNorm; }
+	double NormRHS() { return Norm(b); }
 #if defined(USE_MPI)
 	int DOF() { return A->Dof(); }
 	long Size() { return A->Size(); }