From 4a61e9a1859f21d84d391d8ab87dc6381ca17b67 Mon Sep 17 00:00:00 2001
From: Klaus <klaus.boehnlein@tu-dresden.de>
Date: Wed, 27 Jul 2022 17:15:37 +0200
Subject: [PATCH] Check Solvers.. bad results with CG detected! use QR-solver
instead
---
dune/microstructure/matrix_operations.hh | 26 ++++++
inputs/cellsolver.parset | 25 +++---
src/Cell-Problem.cc | 103 ++++++++++++++---------
3 files changed, 105 insertions(+), 49 deletions(-)
diff --git a/dune/microstructure/matrix_operations.hh b/dune/microstructure/matrix_operations.hh
index 35b80833..7efd263b 100644
--- a/dune/microstructure/matrix_operations.hh
+++ b/dune/microstructure/matrix_operations.hh
@@ -105,6 +105,24 @@ namespace MatrixOperations {
auto t1 = 2.0 * mu * sym(E1) + lambda * trace(sym(E1)) * Id();
return scalarProduct(t1,sym(E2));
}
+//
+// static double linearizedStVenantKirchhoffDensity(double mu, double lambda, MatrixRT E1, MatrixRT E2) // CHANGED
+// {
+// auto t1 = 2.0 * mu * sym(E1) + lambda * trace(sym(E1)) * Id();
+// auto tmp1 = scalarProduct(t1,sym(E2));
+//
+// auto t2 = 2.0 * mu * sym(E2) + lambda * trace(sym(E2)) * Id();
+// auto tmp2 = scalarProduct(t2,sym(E1));
+//
+// if(abs(tmp1-tmp2) > 1e-8 )
+// {
+// std::cout << "linearizedStVenantKirchhoffDensity NOT Symmetric!" << std::endl;
+// std::cout << "tmp1: " << tmp1 << std::endl;
+// std::cout << "tmp2: " << tmp2 << std::endl;
+// }
+// return tmp1;
+//
+// }
// --- Generalization: Define Quadratic QuadraticForm
@@ -117,6 +135,14 @@ namespace MatrixOperations {
}
+// static double linearizedStVenantKirchhoffDensity(const double mu, const double lambda, MatrixRT F, MatrixRT G){
+// /// Write this whole File as a Class that uses lambda,mu as members ?
+//
+// // Define L via Polarization-Identity from QuadratifForm
+// // <LF,G> := (1/2)*(Q(F+G) - Q(F) - Q(G) )
+// return (1.0/2.0)*(QuadraticForm(mu,lambda,F+G) - QuadraticForm(mu,lambda,F) - QuadraticForm(mu,lambda,G) );
+// }
+
static double generalizedDensity(const double mu, const double lambda, MatrixRT F, MatrixRT G){
/// Write this whole File as a Class that uses lambda,mu as members ?
diff --git a/inputs/cellsolver.parset b/inputs/cellsolver.parset
index f2bb27d6..0ff8aa2a 100644
--- a/inputs/cellsolver.parset
+++ b/inputs/cellsolver.parset
@@ -9,8 +9,8 @@
# Choose Output-path for logfile
#############################################
### Remove/Comment this when running via Python-Script:
-outputPath=/home/stefan/DUNE/dune-microstructure/outputs
-#outputPath=/home/klaus/Desktop/DUNE/dune-microstructure/outputs
+#outputPath=/home/stefan/DUNE/dune-microstructure/outputs
+outputPath=/home/klaus/Desktop/Dune-Testing/dune-microstructure/outputs
#############################################
@@ -27,7 +27,7 @@ cellDomain=1
## {start,finish} computes on all grid from 2^(start) to 2^finish refinement
#----------------------------------------------------
-numLevels=2 5
+numLevels= 2 4
#numLevels = 1 1 # computes all levels from first to second entry
#numLevels = 2 2 # computes all levels from first to second entry
#numLevels = 1 3 # computes all levels from first to second entry
@@ -41,7 +41,7 @@ numLevels=2 5
########################################################################################
# --- Choose scale ratio gamma:
-gamma=1.0
+gamma=0.1
#############################################
@@ -59,19 +59,22 @@ alpha=8.0
#--- Lame-Parameters
-mu1=1.0
-lambda1=1.0
+mu1=8.0
+lambda1=5.0
# better: pass material parameters as a vector
# Poisson ratio nu = 1/4 => lambda = 0.4*mu
mu=1.2 1.2 1
-lambda=0.48 0.48 0.4
+#lambda=0.48 0.48 0.4
#rho=1.0 0
+#mu=80 80 60
+lambda=80 80 25
+
# ---volume fraction (default value = 1.0/4.0)
-theta=0.5
+theta=0.25
#theta = 0.3
#theta = 0.75
#theta = 0.125
@@ -80,8 +83,8 @@ theta=0.5
#--- choose composite-Implementation:
-geometryFunctionPath = /home/stefan/DUNE/dune-microstructure/geometries
-#geometryFunctionPath = /home/klaus/Desktop/DUNE/dune-microstructure/geometries
+#geometryFunctionPath = /home/stefan/DUNE/dune-microstructure/geometries
+geometryFunctionPath = /home/klaus/Desktop/Dune-Testing/dune-microstructure/geometries
material_prestrain_imp= "material_neukamm" #(Python-indicator-function with same name determines material phases)
#material_prestrain_imp= "two_phase_material_2" #(Python-indicator-function with same name determines material phases)
#material_prestrain_imp= "two_phase_material_3" #(Python-indicator-function with same name determines material phases)
@@ -99,7 +102,7 @@ material_prestrain_imp= "material_neukamm" #(Python-indicator-function with sa
# --- (Optional output) write Material / prestrain / Corrector functions to .vtk-Files (default=false):
write_materialFunctions = true
write_prestrainFunctions = true # VTK norm of B ,
-#write_VTK = true
+write_VTK = true
# --- (Optional output) L2Error, compute integral mean:
diff --git a/src/Cell-Problem.cc b/src/Cell-Problem.cc
index 82b4a5ed..c8c4561a 100644
--- a/src/Cell-Problem.cc
+++ b/src/Cell-Problem.cc
@@ -424,6 +424,8 @@ void computeElementLoadVector( const LocalView& localView,
std::vector< FieldVector< double, dim> > gradients(referenceGradients.size());
for (size_t i=0; i< gradients.size(); i++)
jacobian.mv(referenceGradients[i][0], gradients[i]);
+
+
// "f*phi"-part
for (size_t i=0; i < nSf; i++)
@@ -439,6 +441,8 @@ void computeElementLoadVector( const LocalView& localView,
defGradientV = crossSectionDirectionScaling((1/gamma),defGradientV);
double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos),forceTerm(quadPos), defGradientV );
+// double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos),forceTerm(element.geometry().global(quadPos)), defGradientV ); //TEST
+
size_t row = localView.tree().child(k).localIndex(i);
elementRhs[row] += energyDensity * quadPoint.weight() * integrationElement;
}
@@ -447,6 +451,7 @@ void computeElementLoadVector( const LocalView& localView,
for (size_t m=0; m<3; m++)
{
double energyDensityfG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), forceTerm(quadPos),basisContainer[m] );
+// double energyDensityfG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), forceTerm(element.geometry().global(quadPos)),basisContainer[m] );//TEST
elementRhs[localPhiOffset+m] += energyDensityfG * quadPoint.weight() * integrationElement;
}
}
@@ -964,8 +969,8 @@ auto test_derivative(const Basis& basis,
auto geometry = e.geometry();
const auto& localFiniteElement = localView.tree().child(0).finiteElement();
- int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1 + 5 ); // TEST
-// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
+// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1 + 5 ); // TEST
+ int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(e.type(), orderQR);
for (const auto& quadPoint : quad)
@@ -1141,8 +1146,8 @@ auto check_Orthogonality(const Basis& basis,
auto geometry = e.geometry();
const auto& localFiniteElement = localView.tree().child(0).finiteElement();
- int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1 + 5 ); // TEST
-// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
+// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1 + 5 ); // TEST
+ int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(e.type(), orderQR);
for (const auto& quadPoint : quad)
@@ -1169,6 +1174,8 @@ auto check_Orthogonality(const Basis& basis,
// strain1[m][n] += M[m][n];
auto G = matrixFieldG(quadPos);
+// auto G = matrixFieldG(e.geometry().global(quadPos)); //TEST
+
auto tmp = G + *M1 + strain1;
@@ -1250,8 +1257,8 @@ auto computeFullQ(const Basis& basis,
auto geometry = e.geometry();
const auto& localFiniteElement = localView.tree().child(0).finiteElement();
- int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1 + 5 ); // TEST
-// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
+// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1 + 5 ); // TEST
+ int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(e.type(), orderQR);
for (const auto& quadPoint : quad)
@@ -1279,11 +1286,13 @@ auto computeFullQ(const Basis& basis,
auto G1 = matrixFieldG1(quadPos);
auto G2 = matrixFieldG2(quadPos);
+// auto G1 = matrixFieldG1(e.geometry().global(quadPos)); //TEST
+// auto G2 = matrixFieldG2(e.geometry().global(quadPos)); //TEST
auto X1 = G1 + Chi1;
auto X2 = G2 + Chi2;
- double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), tmp1, tmp2);
+ double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), X1, X2);
elementEnergy += energyDensity * quadPoint.weight() * integrationElement;
@@ -1386,6 +1395,14 @@ int main(int argc, char *argv[])
// double len = parameterSet.get<double>("len", 10.0); //length
// double height = parameterSet.get<double>("h", 0.1); //rod thickness
// double eps = parameterSet.get<double>("eps", 0.1); //size of perioticity cell
+
+ if(imp == "material_neukamm")
+ {
+ std::cout <<"mu: " <<parameterSet.get<std::array<double,3>>("mu", {1.0,3.0,2.0}) << std::endl;
+ std::cout <<"lambda: " << parameterSet.get<std::array<double,3>>("lambda", {1.0,3.0,2.0}) << std::endl;
+ }
+
+
///////////////////////////////////
// Get Prestrain/Parameters
@@ -1604,13 +1621,13 @@ int main(int argc, char *argv[])
//TEST
- std::cout << "Test crossSectionDirectionScaling:" << std::endl;
-
+// std::cout << "Test crossSectionDirectionScaling:" << std::endl;
+/*
MatrixRT T {{1.0, 1.0, 1.0}, {1.0, 1.0, 1.0}, {1.0, 1.0, 1.0}};
printmatrix(std::cout, T, "Matrix T", "--");
auto ST = crossSectionDirectionScaling((1.0/5.0),T);
- printmatrix(std::cout, ST, "scaled Matrix T", "--");
+ printmatrix(std::cout, ST, "scaled Matrix T", "--");*/
//TEST
// auto QuadraticForm = [] (const double mu, const double lambda, const MatrixRT& M) {
@@ -1661,9 +1678,9 @@ int main(int argc, char *argv[])
// printmatrix(std::cout, stiffnessMatrix_CE, "StiffnessMatrix", "--");
// printvector(std::cout, load_alpha1, "load_alpha1", "--");
-
+ //TODO Add Option for this
// CHECK SYMMETRY:
- checkSymmetry(Basis_CE,stiffnessMatrix_CE);
+// checkSymmetry(Basis_CE,stiffnessMatrix_CE);
// set one basis-function to zero
@@ -1700,7 +1717,9 @@ int main(int argc, char *argv[])
ilu0, //NULL,
1e-8, // desired residual reduction factorlack
iter, // maximum number of iterations
- Solver_verbosity); // verbosity of the solver
+ Solver_verbosity,
+ true // Try to estimate condition_number
+ ); // verbosity of the solver
InverseOperatorResult statistics;
std::cout << "solve linear system for x_1.\n";
solver.apply(x_1, load_alpha1, statistics);
@@ -1709,6 +1728,11 @@ int main(int argc, char *argv[])
std::cout << "solve linear system for x_3.\n";
solver.apply(x_3, load_alpha3, statistics);
log << "Solver-type used: " <<" CG-Solver" << std::endl;
+
+
+ std::cout << "statistics.converged " << statistics.converged << std::endl;
+ std::cout << "statistics.condition_estimate: " << statistics.condition_estimate << std::endl;
+ std::cout << "statistics.iterations: " << statistics.iterations << std::endl;
}
////////////////////////////////////////////////////////////////////////////////////
@@ -1736,7 +1760,7 @@ int main(int argc, char *argv[])
InverseOperatorResult statistics;
// solve for different Correctors (alpha = 1,2,3)
- solver.apply(x_1, load_alpha1, statistics);
+ solver.apply(x_1, load_alpha1, statistics); //load_alpha1 now contains the corresponding residual!!
solver.apply(x_2, load_alpha2, statistics);
solver.apply(x_3, load_alpha3, statistics);
log << "Solver-type used: " <<" GMRES-Solver" << std::endl;
@@ -1993,40 +2017,43 @@ int main(int argc, char *argv[])
- //VARIANT 2
+ //---VARIANT 2
//Compute effective elastic law Q
MatrixRT Q_2(0);
for(size_t a = 0; a < 3; a++)
for(size_t b=0; b < 3; b++)
{
+ std::cout << "check_Orthogonality:" << check_Orthogonality(Basis_CE, muLocal, lambdaLocal,gamma,mContainer[a],mContainer[b],*phiDerContainer[a],*phiDerContainer[b],x3MatrixBasis[a]) << std::endl;
Q_2[a][b] = computeFullQ(Basis_CE, muLocal, lambdaLocal,gamma,mContainer[a],mContainer[b],*phiDerContainer[a],*phiDerContainer[b],x3MatrixBasis[a],x3MatrixBasis[b]);
}
printmatrix(std::cout, Q_2, "Matrix Q_2", "--");
- // VARIANT 3
+ Q = Q_2;
+
+ //--- VARIANT 3
// Compute effective elastic law Q
- MatrixRT Q_3(0);
- for(size_t a = 0; a < 3; a++)
- for(size_t b=0; b < 3; b++)
- {
- assembleCellLoad(Basis_CE, muLocal, lambdaLocal, gamma, tmp1 ,x3MatrixBasis[b]); // <L i(M_alpha) + sym(grad phi_alpha), i(x3G_beta) >
-
- double GGterm = 0.0;
- GGterm = energy(Basis_CE, muLocal, lambdaLocal, x3MatrixBasis[a] , x3MatrixBasis[b] ); // <L i(x3G_alpha) , i(x3G_beta) >
-
- // TEST
- // std::setprecision(std::numeric_limits<float>::digits10);
-
- Q_3[a][b] = (coeffContainer[a]*tmp1) + GGterm; // seems symmetric...check positiv definitness?
-
- if (print_debug)
- {
- std::cout << "analyticGGTERM:" << (mu1*(1-theta)+mu2*theta)/6.0 << std::endl;
- std::cout << "GGTerm:" << GGterm << std::endl;
- std::cout << "coeff*tmp: " << coeffContainer[a]*tmp1 << std::endl;
- }
- }
- printmatrix(std::cout, Q_3, "Matrix Q_3", "--");
+// MatrixRT Q_3(0);
+// for(size_t a = 0; a < 3; a++)
+// for(size_t b=0; b < 3; b++)
+// {
+// assembleCellLoad(Basis_CE, muLocal, lambdaLocal, gamma, tmp1 ,x3MatrixBasis[b]); // <L i(M_alpha) + sym(grad phi_alpha), i(x3G_beta) >
+//
+// double GGterm = 0.0;
+// GGterm = energy(Basis_CE, muLocal, lambdaLocal, x3MatrixBasis[a] , x3MatrixBasis[b] ); // <L i(x3G_alpha) , i(x3G_beta) >
+//
+// // TEST
+// // std::setprecision(std::numeric_limits<float>::digits10);
+//
+// Q_3[a][b] = (coeffContainer[a]*tmp1) + GGterm; // seems symmetric...check positiv definitness?
+//
+// if (print_debug)
+// {
+// std::cout << "analyticGGTERM:" << (mu1*(1-theta)+mu2*theta)/6.0 << std::endl;
+// std::cout << "GGTerm:" << GGterm << std::endl;
+// std::cout << "coeff*tmp: " << coeffContainer[a]*tmp1 << std::endl;
+// }
+// }
+// printmatrix(std::cout, Q_3, "Matrix Q_3", "--");
--
GitLab