diff --git a/inputFiles/constitutiveDriver/friction/frictionDriver_Coulomb.xml b/inputFiles/constitutiveDriver/friction/frictionDriver_Coulomb.xml new file mode 100644 index 00000000000..a24ffec9f26 --- /dev/null +++ b/inputFiles/constitutiveDriver/friction/frictionDriver_Coulomb.xml @@ -0,0 +1,32 @@ + + + + + + + + + + + + diff --git a/inputFiles/constitutiveDriver/friction/frictionDriver_base.xml b/inputFiles/constitutiveDriver/friction/frictionDriver_base.xml new file mode 100644 index 00000000000..e9fe6c2fe4f --- /dev/null +++ b/inputFiles/constitutiveDriver/friction/frictionDriver_base.xml @@ -0,0 +1,63 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/constitutiveDriver/friction/tables/constant.geos b/inputFiles/constitutiveDriver/friction/tables/constant.geos new file mode 100644 index 00000000000..f53496892fa --- /dev/null +++ b/inputFiles/constitutiveDriver/friction/tables/constant.geos @@ -0,0 +1,2 @@ +1e2 +1e2 diff --git a/inputFiles/constitutiveDriver/friction/tables/djumps.geos b/inputFiles/constitutiveDriver/friction/tables/djumps.geos new file mode 100644 index 00000000000..9aea9e0ce5f --- /dev/null +++ b/inputFiles/constitutiveDriver/friction/tables/djumps.geos @@ -0,0 +1,2 @@ +0 +0 \ No newline at end of file diff --git a/inputFiles/constitutiveDriver/friction/tables/jumps.geos b/inputFiles/constitutiveDriver/friction/tables/jumps.geos new file mode 100644 index 00000000000..d9a5218f002 --- /dev/null +++ b/inputFiles/constitutiveDriver/friction/tables/jumps.geos @@ -0,0 +1,2 @@ +0 +-5e-3 diff --git a/inputFiles/constitutiveDriver/friction/tables/time.geos b/inputFiles/constitutiveDriver/friction/tables/time.geos new file mode 100644 index 00000000000..fd3306644d5 --- /dev/null +++ b/inputFiles/constitutiveDriver/friction/tables/time.geos @@ -0,0 +1,2 @@ +0 +5 diff --git a/inputFiles/constitutiveDriver/friction/tables/tractions.geos b/inputFiles/constitutiveDriver/friction/tables/tractions.geos new file mode 100644 index 00000000000..8fd577fe1f0 --- /dev/null +++ b/inputFiles/constitutiveDriver/friction/tables/tractions.geos @@ -0,0 +1,2 @@ +0 +-5e2 diff --git a/src/coreComponents/constitutive/contact/CoulombFriction.hpp b/src/coreComponents/constitutive/contact/CoulombFriction.hpp index 31fd39cadc4..071cbdeb679 100644 --- a/src/coreComponents/constitutive/contact/CoulombFriction.hpp +++ b/src/coreComponents/constitutive/contact/CoulombFriction.hpp @@ -189,6 +189,15 @@ class CoulombFriction : public FrictionBase */ KernelWrapper createKernelUpdates() const; + // TODO (jafranc) make sure it is only helpers for FrictionDriver + /// getting cohesion value + real64 getCohesion() const + { return m_cohesion[0]; } + + /// getting friction coeff + real64 getFrictionCoeff() const + { return m_frictionCoefficient[0]; } + /** * @struct Set of "char const *" and keys for data specified in this class. */ diff --git a/src/coreComponents/constitutiveDrivers/CMakeLists.txt b/src/coreComponents/constitutiveDrivers/CMakeLists.txt index e3b7a38416c..3b6268fe2e7 100644 --- a/src/coreComponents/constitutiveDrivers/CMakeLists.txt +++ b/src/coreComponents/constitutiveDrivers/CMakeLists.txt @@ -29,12 +29,16 @@ set( constitutiveDrivers_headers relativePermeability/RelpermDriver.hpp relativePermeability/RelpermDriverRunTest.hpp solid/TriaxialDriver.hpp + contact/FrictionDriver.hpp + contact/FrictionDriverRunTest.hpp ) # # Specify all sources # set( constitutiveDrivers_sources ConstitutiveDriver.cpp + contact/FrictionDriver.cpp + contact/FrictionDriverRunTest.cpp fluid/multiFluid/PVTDriver.cpp fluid/multiFluid/constant/PVTDriverRunTestInvariantImmiscibleFluid.cpp fluid/multiFluid/blackOil/PVTDriverRunTestDeadOilFluid.cpp diff --git a/src/coreComponents/constitutiveDrivers/ConstitutiveDriver.cpp b/src/coreComponents/constitutiveDrivers/ConstitutiveDriver.cpp index 2126e23da78..8080e7c197a 100644 --- a/src/coreComponents/constitutiveDrivers/ConstitutiveDriver.cpp +++ b/src/coreComponents/constitutiveDrivers/ConstitutiveDriver.cpp @@ -145,9 +145,9 @@ void ConstitutiveDriver::outputToFile() const integer const precision = LvArray::math::max( LvArray::math::min( m_precision, maxPrecision ), minPrecision ); file << std::scientific << std::setprecision( precision ); - integer const width = precision+7; + integer const width = precision+7+1;//spacing for minus sign - for( integer step = 0; step <= m_numSteps; ++step ) + for( integer step = 0; step < m_table.size(0); ++step ) { file << std::setw( width ) << m_table( step, 0 ); for( integer col = 1; col < numColumns; ++col ) @@ -213,7 +213,7 @@ void ConstitutiveDriver::outputToConsole() const TableData tableData; stdVector< TableData::CellData > tableRow( numColumns ); - for( integer step = 0; step <= m_numSteps; ++step ) + for( integer step = 0; step < m_table.size(0); ++step ) { for( integer col = 0; col < numColumns; ++col ) { diff --git a/src/coreComponents/constitutiveDrivers/contact/FrictionDriver.cpp b/src/coreComponents/constitutiveDrivers/contact/FrictionDriver.cpp new file mode 100644 index 00000000000..fae495bbe2a --- /dev/null +++ b/src/coreComponents/constitutiveDrivers/contact/FrictionDriver.cpp @@ -0,0 +1,972 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +#include "FrictionDriver.hpp" + +#include "constitutive/ConstitutiveManager.hpp" +#include "constitutiveDrivers/LogLevelsInfo.hpp" +#include "constitutive/contact/FrictionBase.hpp" +#include "constitutive/contact/FrictionSelector.hpp" + +#include "functions/FunctionManager.hpp" +#include "functions/TableFunction.hpp" + +#define TEST_PBQ 0 + +namespace geos +{ + + namespace LUsolver{ + +bool luDecompose(std::vector& A, int n, std::vector& piv) +{ + piv.resize(n); + for (int i = 0; i < n; ++i) piv[i] = i; + + for (int k = 0; k < n; ++k) + { + // partial pivot: largest |entry| in column k, rows k..n-1 + int p = k; + double maxVal = std::abs(A[k*n + k]); + for (int i = k+1; i < n; ++i) + { + double v = std::abs(A[i*n + k]); + if (v > maxVal) { maxVal = v; p = i; } + } + if (maxVal < 1e-13) return false; // singular / near-singular + + if (p != k) + { + for (int j = 0; j < n; ++j) std::swap(A[k*n+j], A[p*n+j]); + std::swap(piv[k], piv[p]); + } + + for (int i = k+1; i < n; ++i) + { + double factor = A[i*n + k] / A[k*n + k]; + A[i*n + k] = factor; // store multiplier (L entry) + for (int j = k+1; j < n; ++j) + A[i*n + j] -= factor * A[k*n + j]; + } + } + return true; +} + + +// Solve A x = b given the LU-factored A (from luDecompose) and its pivot. +std::vector luSolve(const std::vector& LU, int n, + const std::vector& piv, + const std::vector& b) +{ + std::vector x(n); + for (int i = 0; i < n; ++i) x[i] = b[piv[i]]; // apply permutation + + // forward substitution, L has implicit unit diagonal + for (int i = 0; i < n; ++i) + for (int j = 0; j < i; ++j) + x[i] -= LU[i*n + j] * x[j]; + + // backward substitution + for (int i = n-1; i >= 0; --i) + { + for (int j = i+1; j < n; ++j) + x[i] -= LU[i*n + j] * x[j]; + x[i] /= LU[i*n + i]; + } + return x; +} + +// Convenience: solve a SMALL fixed-size system (used for the 3x3 bubble +// condensation) via the same LU machinery, working on a flat copy. +template +std::array,N> invertSmall(const std::array,N>& M) +{ + std::vector A(N*N); + for (int i=0;i piv; + bool ok = luDecompose(A, N, piv); + if (!ok) std::cerr << "WARNING: singular bubble block!\n"; + + std::array,N> Inv{}; + for (int col=0; col e(N, 0.0); e[col] = 1.0; + std::vector x = luSolve(A, N, piv, e); + for (int row=0; row buildD(double K, double G) + { + array2d D{6,6}; + + D[0][0] = D[1][1] = D[2][2] = K + 4/3*G; + D[0][1] = D[1][0] = D[0][2] = D[2][0] = D[1][2] = D[2][1] = K - 2./3 * G; + D[3][3] = D[4][4] = D[5][5] = G; + + + return D; + } + +void shapeAndGrad(double xi, double eta, double zeta, + std::array& N, + std::array,8>& dNdXi) +{ + for (int a = 0; a < 8; ++a) + { + double xa = refCorners[a][0], ya = refCorners[a][1], za = refCorners[a][2]; + N[a] = 0.125 * (1+xa*xi) * (1+ya*eta) * (1+za*zeta); + dNdXi[a][0] = 0.125 * xa * (1+ya*eta) * (1+za*zeta); + dNdXi[a][1] = 0.125 * (1+xa*xi) * ya * (1+za*zeta); + dNdXi[a][2] = 0.125 * (1+xa*xi) * (1+ya*eta) * za; + } +} + +// Face bubble (Frigo Eq. 20a) for the face at xi_j = sign, on an axis-aligned +// hex. faceAxis in {0,1,2} = which parent coordinate is pinned to the face; +// faceSign = +1 or -1. +void bubbleAndGrad(double xi, double eta, double zeta, + int faceAxis, double faceSign, + double& b, array1d& dbdXi) +{ + std::array c{xi, eta, zeta}; + double linear = 0.5 * (1.0 + faceSign * c[faceAxis]); + double dlinear = 0.5 * faceSign; // wrt c[faceAxis] + + // product of (1 - c_i^2) over the two axes orthogonal to faceAxis + std::array other; + { int k=0; for (int i=0;i<3;++i) if (i!=faceAxis) other[k++]=i; } + + double q0 = 1.0 - c[other[0]]*c[other[0]]; + double q1 = 1.0 - c[other[1]]*c[other[1]]; + + b = linear * q0 * q1; + + array1d grad{3}; + grad[faceAxis] = dlinear * q0 * q1; + grad[other[0]] = linear * (-2.0 * c[other[0]]) * q1; + grad[other[1]] = linear * q0 * (-2.0 * c[other[1]]); + dbdXi = grad;//[1x3] +} + +// Voigt B-column for a generic scalar shape function's physical gradient +array2d bColumn(double dNdx, double dNdy, double dNdz) +{ + array2d Bc{6,3}; Bc.zero(); + // Bc.resize(6,3); + + Bc(0,0)=dNdx; Bc(1,1)=dNdy; Bc(2,2)=dNdz; + Bc(3,1)=dNdz; Bc(3,2)=dNdy; + Bc(4,0)=dNdz; Bc(4,2)=dNdx; + Bc(5,0)=dNdy; Bc(5,1)=dNdx; + return Bc; +} + +array2d computeJacobian(const std::array,8>& coords, + const std::array,8>& dNdXi) +{ + array2d J{3,3}; J.zero(); + for (int a=0;a<8;++a) + { + array1d dN{3}; dN[0] = dNdXi[a][0]; dN[1] = dNdXi[a][1]; dN[2] = dNdXi[a][2]; + LvArray::tensorOps::Rij_add_AiBj< 3, 3 >( J, coords[a], dN ); // J += coords[a] (x) dNdXi[a] + } + return J; +} + + +LocalBlocks integrateElement(const array2d& D, const HexElem& e) +{ + LocalBlocks L; + L.Kuu.zero(); L.Kub.zero(); L.Kbb.zero(); + + double detJ = 0.; + const double w = 1.0; + + for (const auto& gpt : gaussPts) + { + double xi=gpt[0], eta=gpt[1], zeta=gpt[2]; + std::array N; + std::array,8> dNdXi; + shapeAndGrad(xi, eta, zeta, N, dNdXi);//TODO use GEOS's function as in H1_TriangleFace_Lagrange1_Gauss.hpp + + array2d mat(3,3); mat.zero(); + array2d mat36(6,3); mat36.zero(); + array2d Ba[8]; + // L.Kuu.resize(L.nnodes,L.nnodes); + // L.Kbb.resize(L.ncenters,L.ncenters); + // L.Kub.resize(L.nnodes,L.ncenters); + + array2d J = computeJacobian(e.coords, dNdXi), invJ{3,3}; + detJ = LvArray::tensorOps::invert<3>(invJ,J); + + for (int a=0;a<8;++a){ + Ba[a].resize(6,3); + array1d g{3}; + array1d dN{3}; + dN[0] = dNdXi[a][0]; dN[1] = dNdXi[a][1]; dN[2] = dNdXi[a][2]; + LvArray::tensorOps::Ri_eq_AijBj<3,3>(g, invJ, dN); + Ba[a] = bColumn(g[0],g[1],g[2]); + } + + for (int a=0;a<8;++a) + for (int b=0;b<8;++b){ + // L.Kuu.block<3,3>(3*a,3*b) += gaussW*detJ * (Ba[a].transpose()*D*Ba[b]); + mat.zero(); + #if TEST_PBQ + LvArray::tensorOps::Rij_eq_PkiBklQlj<6,3>(mat,Ba[a],D,Ba[b]); + #else + mat36.zero(); + LvArray::tensorOps::Rij_eq_AkiBkj<3,6,6>(mat36, Ba[a],D); + LvArray::tensorOps::Rij_eq_AikBkj<3,3,6>(mat,mat36,Ba[b]); + #endif + LvArray::tensorOps::scale<3,3>(mat, w*gaussW*detJ); + for (int i=0; i<3; ++i) + for (int j=0; j<3; ++j) + L.Kuu[3*a + i][3*b + j ] += mat[i][j]; + } + + if (e.hasBubble) + { + double bub; array1d dbdXi; + bubbleAndGrad(xi,eta,zeta, e.bubbleFaceAxis, e.bubbleFaceSign, bub, dbdXi); + array1d gb{3}; + array1d db{3}; + db[0] = dbdXi[0]; db[1] = dbdXi[1]; db[2] = dbdXi[2]; + LvArray::tensorOps::Ri_eq_AjiBj<3,3>(gb, invJ, db); + auto Bb = bColumn(gb[0],gb[1],gb[2]); + + for (int a=0;a<8;++a){ + // L.Kub.block<3,3>(3*a,0) += gaussW*detJ * (Ba[a].transpose()*D*Bb); + mat.zero(); + #if TEST_PBQ + LvArray::tensorOps::Rij_eq_PkiBklQlj<6,3>(mat,Ba[a],D,Bb); + #else + mat36.zero(); + LvArray::tensorOps::Rij_eq_AkiBkj<3,6,6>(mat36, Ba[a],D); + LvArray::tensorOps::Rij_eq_AikBkj<3,3,6>(mat,mat36,Bb); + #endif + LvArray::tensorOps::scale<3,3>(mat, w*gaussW*detJ); + for (int i=0; i<3; ++i) + for (int j=0; j<3; ++j) + L.Kuu[3*a + i][ j ] += mat[i][j]; + } + + + + // L.Kbb += gaussW*detJ * (Bb.transpose()*D*Bb); + mat.zero(); + #if TEST_PBQ + LvArray::tensorOps::Rij_eq_PkiBklPlj<6,3>(mat,Bb,D); + #else + mat36.zero(); + LvArray::tensorOps::Rij_eq_AkiBkj<3,6,6>(mat36, Bb,D); + LvArray::tensorOps::Rij_eq_AikBkj<3,3,6>(mat,mat36,Bb); + #endif + LvArray::tensorOps::scale<3,3>(mat, w*gaussW*detJ); + LvArray::tensorOps::add<3,3>(L.Kbb, mat); + } + } + return L; +} + +// Kcond = Kuu - Kub * Kbb^-1 * Kub^T (3x3 inverse via our own LU, small enough +// that we go through the same generic solver rather than a closed form) +array2d condense(const LocalBlocks& L) +{ + // std::vector Kbb3(9); + // for (int i=0;i<3;++i) for (int j=0;j<3;++j) Kbb3[i*3+j] = L.Kbb(i,j); + // std::vector piv; + // luDecompose(Kbb3, 3, piv); + + array2d KbbInv(3,3); + LvArray::tensorOps::invert(KbbInv, L.Kbb); + // for (int col=0; col<3; ++col) + // { + // std::vector e(3,0.0); e[col]=1.0; + // auto x = luSolve(Kbb3, 3, piv, e); + // for (int row=0; row<3; ++row) KbbInv(row,col) = x[row]; + // } + + array2d Kcond(24,24), mat(24,24), mat243(24,3); + LvArray::tensorOps::copy(Kcond,L.Kuu); + // Eigen::Matrix Kcond = L.Kuu - L.Kub * KbbInv * L.Kub.transpose(); + //TODO (recheck formula) + #if TEST_PBQ + LvArray::tensorOps::Rij_eq_PikBklQjl(mat, L.Kub, KbbInv, L.Kub); + LvArray::tensorOps::scaledAdd(Kcond, mat, -1); + #else// main road + LvArray::tensorOps::Rij_eq_AikBkj(mat243, L.Kub, KbbInv); + LvArray::tensorOps::Rij_eq_AikBjk(mat, mat243, L.Kub); + LvArray::tensorOps::scaledAdd(Kcond, mat, -1); + #endif + + // for (int r=0;r<24;++r) + // for (int c=0;c<24;++c) + // { + // double s=0.0; + // for (int k1=0;k1<3;++k1) + // for (int k2=0;k2<3;++k2) + // s += L.Kub(r,k1)*KbbInv(k1,k2)*L.Kub(c,k2); + // Kcond(r,c) = L.Kuu(r,c) - s; + // } + return Kcond; +} + +CRSMat assemble(const array2d& D, const ContactState& contact) +{ + array2d R{3,3}; + HexElem elemA, elemB; + setElement(contact, elemA,elemB,R); + + const int nNodes = 16, nDof = nNodes*3; + + // ---------------- assemble & condense each element --------------- + // Eigen::Matrix KcondA = assembleCondensed(elemA); + // Eigen::Matrix KcondB = assembleCondensed(elemB); + array2d KcondA = mimicAssembly::condense(integrateElement(D,elemA)); + array2d KcondB = mimicAssembly::condense(integrateElement(D,elemB)); + + // ---------------- GEOS-style two-phase CRSMatrix assembly ---------------- + // Phase 1 (symbolic): reserve every (row,col) pair that will ever be + // written -- elastic pairs from both elements, PLUS the 4x4 cross-face + // contact coupling pairs (needed even though epsN only activates them + // conditionally, since CRSMatrix pre-allocates its sparsity). + CRSMat Kelastic; + Kelastic.resize(nDof, nDof, 24); // generous initial row capacity + + auto elemDofs = [&](const HexElem& e) + { + std::array dofs; + for (int a=0;a<8;++a) for (int c=0;c<3;++c) dofs[3*a+c] = 3*e.node[a]+c; + return dofs; + }; + auto dofsA = elemDofs(elemA); + auto dofsB = elemDofs(elemB); + + for (auto row : dofsA) for (auto col : dofsA) Kelastic.insertNonZero(row, col, 0.0); + for (auto row : dofsB) for (auto col : dofsB) Kelastic.insertNonZero(row, col, 0.0); + + + // NOTE: the contact tangent (state-dependent, rebuilt every Newton + // iteration) is deliberately NOT part of Kelastic's sparsity -- it's a + // tiny 8x8 dense coupling added directly to the dense Jacobian buffer + // each iteration (see solveStep). Kelastic only ever holds the constant + // elastic operator. + + // Phase 2 helper (numeric): accumulate the CONSTANT elastic contribution + // via addToRow -- exactly the GEOS kernel pattern (columns must be sorted + // per row, which our simple 8-node connectivity already is after a sort). + auto scatterElastic = [&](const std::array& dofs, const array2d& Kc) + { + for (int a=0;a<24;++a) + { + std::vector cols(dofs.begin(), dofs.end()); + std::vector vals(24); + for (int b=0;b<24;++b) vals[b] = Kc(a,b); + // sort cols and vals together (addToRow requires sorted columns) + std::vector order(24); for (int i=0;i<24;++i) order[i]=i; + std::sort(order.begin(), order.end(), [&](int i,int j){return cols[i] colsSorted(24); std::vector valsSorted(24); + for (int i=0;i<24;++i){ colsSorted[i]=cols[order[i]]; valsSorted[i]=vals[order[i]]; } + Kelastic.addToRow< RAJA::seq_atomic >(dofs[a], colsSorted.data(), valsSorted.data(), 24); + } + }; + scatterElastic(dofsA, KcondA); + scatterElastic(dofsB, KcondB); + + // Kelastic now holds the constant, condensed elastic operator. Each + // Newton iteration below extracts it to a dense buffer and adds the + // state-dependent contact tangent on top -- Kelastic itself is never + // touched again. + + return Kelastic; +}; + +void setElement(const ContactState& contact, HexElem& elemA, HexElem& elemB, array2d& RR ) +{ + // elemA.x0=0; elemA.y0=0; elemA.z0=0; + elemA.node = {0,1,2,3,4,5,6,7}; + elemA.hasBubble = true; elemA.bubbleId = 0; + elemA.bubbleFaceAxis = 0; elemA.bubbleFaceSign = +1.0; // xi=+1 face -> x=1 + + // elemB.x0=1; elemB.y0=0; elemB.z0=0; + elemB.node = {8,9,10,11,12,13,14,15}; + elemB.hasBubble = true; elemB.bubbleId = 1; + elemB.bubbleFaceAxis = 0; elemB.bubbleFaceSign = -1.0; // xi=-1 face -> x=1 + + + //rotations and geometry + auto rotZ = [](double a){ array2d R{3,3}; R.zero(); R[0][0]=cos(a);R[0][1]=-sin(a);R[0][2]=0; + R[1][0]=sin(a);R[1][1]=cos(a); R[1][2]=0; + R[2][0]=0; R[2][1]=0; R[2][2]=1; return R; }; + auto rotY = [](double a){ array2d R{3,3}; R.zero(); R[0][0]=cos(a); R[0][1]=0;R[0][2]=sin(a); + R[1][0]=0; R[1][1]=1;R[1][2]=0; + R[2][0]=-sin(a);R[2][1]=0;R[2][2]=cos(a); return R; }; + + + //setters + LvArray::tensorOps::Rij_eq_AikBkj<3,3,3>(RR,rotY(contact.theta),rotZ(contact.phi)); + auto unitCubeCorner = [&](array1d origin, int localIndex) + { + array1d c{3}; c.zero(); + c[0] = origin[0] + 0.5*(1+refCorners[localIndex][0]); + c[1] = origin[1] + 0.5*(1+refCorners[localIndex][1]); + c[2] = origin[2] + 0.5*(1+refCorners[localIndex][2]); + array1d res{3}; + LvArray::tensorOps::Ri_eq_AijBj<3,3>(res, RR, c); + return res; + }; + + + array1d e0{3}, e1{3}; + e0.zero(); e1.zero(); e1[0] = 1.; + for (int a=0;a<8;++a) elemA.coords[a] = unitCubeCorner(e0, a); + for (int a=0;a<8;++a) elemB.coords[a] = unitCubeCorner(e1, a); + +} + +std::vector solveStep(double imposedUx, ContactState& contact, CRSMat const & Kelastic, double epsN, int kMaxNewton, std::function& updateNormalTraction ) + { + std::vector results; + //prep -- duplicate (pass element A and B ?) + HexElem elemA, elemB; + array2d R{3,3}; + setElement(contact, elemA,elemB, R); + + // ---------------- fault normal DERIVED from geometry -------------------- + // Two edge vectors of the (planar, quad) fault face, physical corners + // taken directly from elemA's rotated coordinates -- no assumption + // about which global axis the fault points along. + array1d nf{3}; nf[0] = R[0][0]; nf[1] = R[1][0]; nf[2] = R[2][0]; + + // Orient outward from "-" element into "+" element: check against the + // vector between element centroids and flip if pointing the wrong way. + array1d centroidA{3}, centroidB{3}; + centroidA.zero(); centroidB.zero(); + for (int a=0;a<8;++a) + for (int c=0;c<3;++c) + { centroidA[c]+=elemA.coords[a][c]/8.0; centroidB[c]+=elemB.coords[a][c]/8.0; } + + array1d aToB{3}; aToB.zero(); + aToB[0] = centroidB[0]-centroidA[0]; aToB[1] = centroidB[1]-centroidA[1]; aToB[2] = centroidB[2]-centroidA[2]; + + if (LvArray::tensorOps::AiBi<3>(nf, aToB) < 0.0) { LvArray::tensorOps::scale<3>(nf,-1); } + + std::cout << std::fixed << std::setprecision(6); + std::cout << "Derived fault normal: (" << nf[0] << ", " << nf[1] << ", " << nf[2] << ")\n" + << "Expected (rotated x-axis): (" << R[0][0] << ", " << R[1][0] << ", " << R[2][0] << ")\n\n"; + + + + // const int nNodes = 16, nDof = nNodes*3; + const auto nDof = Kelastic.numRows(); + // elemA reference corners with xi=+1 are local nodes {1,2,5,6} -> fault face "-" + // elemB reference corners with xi=-1 are local nodes {0,3,4,7} -> fault face "+" + std::array faceLocalMinus = {1,2,5,6}; + std::array faceLocalPlus = {0,3,4,7}; + std::array faceGlobalMinus, faceGlobalPlus; + for (int i=0;i<4;++i) + { + // std::cout << " elemA.node : " << elemA.node[faceLocalMinus[i]] << "\n elemB.node : " << elemB.node[faceLocalPlus[i]] << std::endl; + faceGlobalMinus[i] = elemA.node[faceLocalMinus[i]]; + faceGlobalPlus[i] = elemB.node[faceLocalPlus[i]]; + } + std::array,4> faceDofMinus, faceDofPlus; + for (int i=0;i<4;++i) + for (int c=0;c<3;++c) + { + // std::cout << " dofMinus: " << 3*faceGlobalMinus[i] + c << "\n dofPlus : " << 3*faceGlobalPlus[i] + c << std::endl; + faceDofMinus[i][c] = 3*faceGlobalMinus[i]+c; + faceDofPlus[i][c] = 3*faceGlobalPlus[i]+c; + } + + std::vector u(nDof, 0.0); + double tNold = contact.tN; + const double faceWeight = 1.0/4.0; + + auto denseFromCRS = [&](const CRSMat& M) + { + std::vector A(nDof*nDof, 0.0); + auto view = M.toViewConst(); + for (IndexType row=0; row(jumpAvg,nf); + trial = updateNormalTraction(tNold, gN, epsN); + std::cout << "\n[ " << newtonIt << " ] trial : " << contact.tN << std::endl; + } + + + // residual r = K_elastic*u + contact virtual work (tangent is + // NOT part of this -- residual uses ONLY the elastic operator, + // matching Frigo's r_u = (elastic term) + (traction term); the + // tangent stiffness kc belongs to the JACOBIAN only, never to + // the residual's K*u term) + std::vector Kdense = denseFromCRS(Kelastic); + std::vector r(nDof, 0.0); + for (int i=0;i(std::cout,",")); + + for (int i=0;i<4;++i) + for (int c=0;c<3;++c) + { + r[faceDofPlus[i][c]] += faceWeight*trial.tN*nf[c]; + r[faceDofMinus[i][c]] -= faceWeight*trial.tN*nf[c]; + } + + // Jacobian = Kglobal + contact tangent (if closed) + std::vector J = Kdense; + if (!trial.open) + { + double kc = faceWeight*epsN/4.0; + for (int i=0;i<4;++i) + for (int j=0;j<4;++j) + for (int a=0;a<3;++a) + for (int b=0;b<3;++b) + { + int pi=faceDofPlus[i][a], pj=faceDofPlus[j][b]; + int mi=faceDofMinus[i][a], mj=faceDofMinus[j][b]; + J[pi*nDof+pj] += kc * nf[a] * nf[b]; + J[pi*nDof+mj] -= kc* nf[a] * nf[b]; + J[mi*nDof+pj] -= kc* nf[a] * nf[b]; + J[mi*nDof+mj] += kc* nf[a] * nf[b]; + } + } + + // BCs + std::array fixedLocal = {0,3,4,7}; + std::array drivenLocal = {1,2,5,6}; + std::vector fixedDofs, drivenDofsX, drivenDofsYZ; + for (int ln : fixedLocal) + for(int c=0;c<3;++c) + fixedDofs.push_back(3*elemA.node[ln]+c); + + for (int ln : drivenLocal) + { + drivenDofsX.push_back(3*elemB.node[ln]+0); + drivenDofsYZ.push_back(3*elemB.node[ln]+1); + drivenDofsYZ.push_back(3*elemB.node[ln]+2); + } + std::vector isFixed(nDof,false); + for (auto d: fixedDofs) isFixed[d]=true; + for (auto d: drivenDofsYZ) isFixed[d]=true; + for (auto d: drivenDofsX) isFixed[d]=true; + + // apply BCs by row/col elimination + std::vector rhs(nDof); + for (int i=0;i piv; + bool ok = LUsolver::luDecompose(J, (int)nDof, piv); + if (!ok) { std::cerr << "singular Jacobian!\n"; break; } + std::vector du = LUsolver::luSolve(J, (int)nDof, piv, rhs); + + double duNorm = 0.0; + for (int i=0;i(jumpAvg,nf); + results.emplace_back(u, gN, newtonIt, (std::sqrt(duNorm) < 1e-12)); + break; + } + else{ + results.emplace_back(u, gN, newtonIt, (std::sqrt(duNorm) < 1e-12)); + } + } + + return results; + }; + +};//end mimicAssembly + + + +using namespace dataRepository; +using namespace constitutive; + +FrictionDriver::FrictionDriver( const string & name, Group * const parent ) + : ConstitutiveDriver( name, parent ) +{ + registerWrapper( viewKeyStruct::frictionNameString(), &m_frictionName ). + setRTTypeName( rtTypes::CustomTypes::groupNameRef ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "Friction model to test" ); + + registerWrapper( viewKeyStruct::numStepsString(), &m_numSteps ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "Number of sample step to take in both jumps and traction increments" ); + + registerWrapper( viewKeyStruct::displacementFunctionString(), &m_dispFunctionName ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "Name of the input function representing displacement function along world x-axis" ); + + // registerWrapper( viewKeyStruct::dJumpFunctionString(), &m_dJumpFunctionName ).//should be derive from convergence history + // setInputFlag( InputFlags::REQUIRED ). + // setDescription( "Name of the input function representing deltaDisplacementJump function along world x-axis" ); + + registerWrapper( viewKeyStruct::stressFunctionLString(), &m_stressFunctionsNamesL ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "Name of the input functions representing stresses in the left cell"); + + registerWrapper( viewKeyStruct::stressFunctionRString(), &m_stressFunctionsNamesR ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "Name of the input functions representing stresses in the right cell"); + + registerWrapper( viewKeyStruct::thetaString(), &m_theta ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "y-Tilt angle in degree" ); + + registerWrapper( viewKeyStruct::phiString(), &m_phi ). + setInputFlag( InputFlags::OPTIONAL ). + setDefaultValue(0.). + setDescription( "z-Tilt angle in degree" ); + + //first batch of parameters + registerWrapper( viewKeyStruct::normalDispTolFac(), &m_normalDispTolFac ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "normal Displacement Tolerance (scale as inverse of average Young modulus)." ); + + registerWrapper( viewKeyStruct::normalTractionTolFac(), &m_normalTracTolFac ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "normal Traction Tolerance" ); + + registerWrapper( viewKeyStruct::slidingTolFac(), &m_slidingTolFac ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "tangential Displacement Tolerance" ); + + registerWrapper( viewKeyStruct::iterPenNFac(), &m_iterPenNFac ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "normal Penalty Factor" ); + + registerWrapper( viewKeyStruct::iterPenTFac(), &m_iterPenTFac ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "tangential Penatly Factor" ); + + + //geometry + registerWrapper( viewKeyStruct::area(), &m_area ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "Face area" ); + + registerWrapper( viewKeyStruct::volume(), &m_volume ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "Neighboring cells' volume" ); + + registerWrapper( viewKeyStruct::bulk(), &m_bulkModulus ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "Neighboring cells' bulk Modulus" ); + + registerWrapper( viewKeyStruct::shear(), &m_shearModulus ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "Neighboring cells' shear Modulus" ); + + //algo tune + registerWrapper( viewKeyStruct::simultaneous(), &m_isSimultaneous ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "isSimultaneous" ); + + registerWrapper( viewKeyStruct::maxNewtonIterString(), &m_maxNewtonIter ). + setInputFlag( InputFlags::OPTIONAL ). + setApplyDefaultValue( 20 ). + setDescription( "Maximum Newton iterations" ); + + addLogLevel< logInfo::LogOutput >(); +} + +void FrictionDriver::postInputInitialization() +{ + ConstitutiveDriver::postInputInitialization(); + + // Check that the functions exist + FunctionManager & functionManager = FunctionManager::getInstance(); + GEOS_ERROR_IF( !functionManager.hasGroup< TableFunction >( m_dispFunctionName ), + GEOS_FMT( "Jump function with name '{}' not found", m_dispFunctionName ), + getWrapperDataContext( viewKeyStruct::displacementFunctionString() ) ); + + // GEOS_ERROR_IF( !functionManager.hasGroup< TableFunction >( m_dJumpFunctionName ), + // GEOS_FMT( "dJump function with name '{}' not found", m_dJumpFunctionName ), + // getWrapperDataContext( viewKeyStruct::dJumpFunctionString() ) ); + + GEOS_ERROR_IF( !functionManager.hasGroup< TableFunction >( m_stressFunctionsNamesR ), + GEOS_FMT( "Stress functions with name '{}' not found", m_stressFunctionsNamesR ), + getWrapperDataContext( viewKeyStruct::stressFunctionRString() ) ); + + GEOS_ERROR_IF( !functionManager.hasGroup< TableFunction >( m_stressFunctionsNamesL ), + GEOS_FMT( "Stress functions with name '{}' not found", m_stressFunctionsNamesL ), + getWrapperDataContext( viewKeyStruct::stressFunctionLString() ) ); + + string_array columnNames; + getColumnNames( columnNames ); + integer const numCols = static_cast< integer >(columnNames.size()); + + // initialize functions + TableFunction & dispFunction = functionManager.getGroup< TableFunction >( m_dispFunctionName ); + // TableFunction & dJumpFunction = functionManager.getGroup< TableFunction >( m_dJumpFunctionName ); + TableFunction & tractionFunctionR00 = functionManager.getGroup< TableFunction >( m_stressFunctionsNamesR ); + TableFunction & tractionFunctionL00 = functionManager.getGroup< TableFunction >( m_stressFunctionsNamesL ); + // TableFunction & tractionFunction = functionManager.getGroup< TableFunction >( m_tractionFunctionName ); + + dispFunction.initializeFunction(); + // dJumpFunction.initializeFunction(); + tractionFunctionR00.initializeFunction(); + tractionFunctionL00.initializeFunction(); + // tractionFunction.initializeFunction(); + + // // TODO: Maybe we should take the maximum extent of jumpFunction and tractionFunction + ArrayOfArraysView< real64 > coordinates = dispFunction.getCoordinates(); + real64 const minTime = coordinates[0][0]; + real64 const maxTime = coordinates[0][coordinates.sizeOfArray( 0 )-1] * 20; + + // Allocate the data + allocateTable( numCols, minTime, maxTime ); + m_table.resize( (m_numSteps+1)*m_maxNewtonIter, numCols ); + real64 const dt = (maxTime-minTime) / m_numSteps; + + // set time columns + for( integer step = 0; step < m_table.size(0) ; ++step ) + { + m_table( step, TIME ) = minTime + (step/m_maxNewtonIter)*dt; + } + // set input columns + initializeTable(); +} + +bool FrictionDriver::execute() +{ + FrictionBase & baseFriction = getFriction(); + + GEOS_LOG_LEVEL_RANK_0( logInfo::LogOutput, "Launching Friction Driver" ); + GEOS_LOG_LEVEL_RANK_0( logInfo::LogOutput, " Friction ............... " << m_frictionName ); + GEOS_LOG_LEVEL_RANK_0( logInfo::LogOutput, " Type ................... " << baseFriction.getCatalogName() ); + GEOS_LOG_LEVEL_RANK_0( logInfo::LogOutput, " Steps .................. " << m_numSteps ); + GEOS_LOG_LEVEL_RANK_0( logInfo::LogOutput, " Output ................. " << m_outputFile ); + + // create a dummy discretization with one quadrature point for + // storing constitutive data + conduit::Node node; + dataRepository::Group rootGroup( "root", node ); + dataRepository::Group discretization( "discretization", &rootGroup ); + + integer const numRows = m_table.size( 0 ); + discretization.resize( numRows ); // numRows elements + baseFriction.allocateConstitutiveData( discretization, 1 ); // one quadrature point + + constitutiveUpdatePassThru( baseFriction, [&]( auto & selectedFrictionModel ) + { + using FRICTION_TYPE = TYPEOFREF( selectedFrictionModel ); + runTest< FRICTION_TYPE >( selectedFrictionModel, m_table ); + } ); + + return false; +} + +void FrictionDriver::getColumnNames( string_array & columnNames ) const +{ + columnNames.emplace_back( "time" ); + columnNames.emplace_back( "traction,normal" ); + columnNames.emplace_back( "traction,tangent1" ); + columnNames.emplace_back( "traction,tangent2" ); + columnNames.emplace_back( "delta displacement jump,normal" ); + columnNames.emplace_back( "delta displacement jump,tangent1" ); + columnNames.emplace_back( "delta displacement jump,tangent2" ); + columnNames.emplace_back( "encoded constaint (0:converged, 1:stick & gn>0 (opening), 2: interpenetration, 3: stick & gt>lim (disp-sliding), 4: tau>taulim (trac-sliding) )" ); + columnNames.emplace_back( "fracture state (0:stick, 1:slip , 2: new slip, 3: open)" ); + columnNames.emplace_back( "newtraction,normal" ); + columnNames.emplace_back( "newtraction,tangent1" ); + columnNames.emplace_back( "newtraction,tangent2" ); + columnNames.emplace_back( "displacement jump,normal" ); + columnNames.emplace_back( "displacement jump,tangent1" ); + columnNames.emplace_back( "displacement jump,tangent2" ); + columnNames.emplace_back( "derived disp tol, normal" ); + columnNames.emplace_back( "derived disp tol, tangent" ); + columnNames.emplace_back( "derived traction tol, normal" ); + columnNames.emplace_back( "iterative penalty, normal" ); + columnNames.emplace_back( "iterative penalty, tangent" ); + columnNames.emplace_back( "mimicAssembly FEM solve, converged normal traction" ); + columnNames.emplace_back( "mimicAssembly FEM solve, converged normal displacement jump" ); + columnNames.emplace_back( "mimicAssembly FEM solve, Newton iteration count" ); + columnNames.emplace_back( "mimicAssembly FEM solve, converged flag (1=converged, 0=not converged)" ); + columnNames.emplace_back( "iteration" ); + + + if( dynamic_cast< CoulombFriction const * >(&getFriction()) != nullptr ) + { + columnNames.emplace_back( "tau limit" ); + } +} + +void FrictionDriver::initializeTable() +{ + integer const numRows = m_table.size( 0 ); + + FunctionManager & functionManager = FunctionManager::getInstance(); + TableFunction const & dispFunction = functionManager.getGroup< TableFunction >( m_dispFunctionName ); + // TableFunction const & dJumpFunction = functionManager.getGroup< TableFunction >( m_dJumpFunctionName ); + + TableFunction const & tractionFunctionR00 = functionManager.getGroup< TableFunction >( m_stressFunctionsNamesR );//should be Voigt full from stress + // TableFunction const & tractionFunctionR11 = functionManager.getGroup< TableFunction >( m_stressFunctionsNamesR );//1 + // TableFunction const & tractionFunctionR22 = functionManager.getGroup< TableFunction >( m_stressFunctionsNamesR );//2 + // TableFunction const & tractionFunctionR21 = functionManager.getGroup< TableFunction >( m_stressFunctionsNamesR );//3 + // TableFunction const & tractionFunctionR02 = functionManager.getGroup< TableFunction >( m_stressFunctionsNamesR );//4 + // TableFunction const & tractionFunctionR01 = functionManager.getGroup< TableFunction >( m_stressFunctionsNamesR );//5 + // + TableFunction const & tractionFunctionL00 = functionManager.getGroup< TableFunction >( m_stressFunctionsNamesL ); + // TableFunction const & tractionFunctionL11 = functionManager.getGroup< TableFunction >( m_stressFunctionsNamesL ); + // TableFunction const & tractionFunctionL22 = functionManager.getGroup< TableFunction >( m_stressFunctionsNamesL ); + // TableFunction const & tractionFunctionL21 = functionManager.getGroup< TableFunction >( m_stressFunctionsNamesL ); + // TableFunction const & tractionFunctionL02 = functionManager.getGroup< TableFunction >( m_stressFunctionsNamesL ); + // TableFunction const & tractionFunctionL01 = functionManager.getGroup< TableFunction >( m_stressFunctionsNamesL ); + + real64 const cos_theta = cos( m_theta * M_PI/180.0 ); + real64 const sin_theta = sin( m_theta * M_PI/180.0 ); + + real64 const cos_phi = cos( m_phi * M_PI/180.0 ); + real64 const sin_phi = sin( m_phi * M_PI/180.0 ); + + for( integer index = 0; index < numRows; index+=m_maxNewtonIter ) + { + std::cerr << " index: " << index << "/" << numRows << std::endl; + real64 const time = m_table( index, TIME )/20; + std::cerr << " time: " << time << std::endl; + + + real64 const leftStress[6] = { tractionFunctionR00.evaluate(&time), + tractionFunctionR00.evaluate(&time), + tractionFunctionR00.evaluate(&time), + 0.,// tractionFunctionR21.evaluate(&time), + 0.,// tractionFunctionR02.evaluate(&time), + 0.// tractionFunctionR01.evaluate(&time) + };//should be full 6-Voigt + + + real64 const rightStress[6] = { tractionFunctionL00.evaluate(&time), + tractionFunctionL00.evaluate(&time), + tractionFunctionL00.evaluate(&time), + 0.,// tractionFunctionL21.evaluate(&time), + 0.,// tractionFunctionL02.evaluate(&time), + 0.// tractionFunctionL01.evaluate(&time) + };//should be full 6-Voigt + + real64 const n[3] = {cos_theta,sin_theta*sin_phi,-cos_theta*sin_phi}; + + m_table( index, DISP ) = dispFunction.evaluate(&time); + //rotate and project stress + for(const auto& sigma : {leftStress, rightStress }){ + + constexpr real64 NCELLS = 2.0; + + const real64 sigman[3] = { + sigma[0]*n[0]+sigma[5]*n[1]+sigma[4]*n[2], + sigma[5]*n[0]+sigma[1]*n[1]+sigma[3]*n[2], + sigma[4]*n[0]+sigma[3]*n[1]+sigma[2]*n[2], + }; + + m_table( index, NTRAC ) += 1./NCELLS*(sigman[0]*cos_phi + sigman[1]*sin_phi*sin_theta + sigman[2]*cos_theta*sin_phi); + m_table( index, STRAC0 ) += 1./NCELLS*(sigman[1]*cos_theta - sigman[2]*sin_theta); + m_table( index, STRAC1 ) += 1./NCELLS*(-sigman[0]*sin_phi + sigman[1]*sin_theta*cos_phi + sigman[2]*cos_theta*cos_phi); + + } + + //Project and average Left and Right stress + //jump = tLocal/(kn kt) + // real64 const jump = jumpFunction.evaluate( &time ); + // m_table( index, NJUMP ) = m_table( index, NTRAC )/m_iterPenNFac;//might need compute tol --> move to evaluation part + // m_table( index, SLIP0 ) = m_table( index, STRAC0 )/m_iterPenTFac; + // m_table( index, SLIP1 ) = m_table( index, STRAC1 )/m_iterPenTFac; + + // real64 const dJump = dJumpFunction.evaluate( &time ); + // m_table( index, NDJUMP ) = dJump*cos_phi*cos_theta; + // m_table( index, DSLIP0 ) = dJump*cos_theta*sin_phi; + // m_table( index, DSLIP1 ) = dJump*sin_theta; + m_table( index, NDJUMP ) = 0.; + m_table( index, DSLIP0 ) = 0.; + m_table( index, DSLIP1 ) = 0.; + + m_table( index, CC ) = 0; + m_table( index, FS ) = fields::contact::FractureState::Stick; + } + + // if( CoulombFriction const * coulombFriction = dynamic_cast< CoulombFriction const * >(&getFriction()) ) + if( dynamic_cast< CoulombFriction const * >(&getFriction()) ) + { + real64 const cohesion = 0e6;//coulombFriction->getCohesion(); //sized-0 mesh crash --> getDefault + real64 const frictionCoeff = 0.01;//coulombFriction->getFrictionCoeff(); + for( integer index = 0; index < numRows; index+=m_maxNewtonIter ) + { + + real64 const normal_traction = m_table( index, NTRAC ); + m_table( index, TLIM ) = cohesion - normal_traction * frictionCoeff; + } + } +} + +FrictionBase & FrictionDriver::getFriction() +{ + return getConstitutiveManager().getGroup< FrictionBase >( m_frictionName ); +} + +FrictionBase const & FrictionDriver::getFriction() const +{ + return getConstitutiveManager().getGroup< FrictionBase >( m_frictionName ); +} + + +REGISTER_CATALOG_ENTRY( TaskBase, + FrictionDriver, + string const &, dataRepository::Group * const ) + +} diff --git a/src/coreComponents/constitutiveDrivers/contact/FrictionDriver.hpp b/src/coreComponents/constitutiveDrivers/contact/FrictionDriver.hpp new file mode 100644 index 00000000000..77f6eb3d222 --- /dev/null +++ b/src/coreComponents/constitutiveDrivers/contact/FrictionDriver.hpp @@ -0,0 +1,292 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +#ifndef GEOS_CONSTITUTIVEDRIVERS_CONTACT_FRICTIONDRIVER_HPP +#define GEOS_CONSTITUTIVEDRIVERS_CONTACT_FRICTIONDRIVER_HPP + + +#include "constitutiveDrivers/ConstitutiveDriver.hpp" +#include "physicsSolvers/solidMechanics/contact/ContactSolverBase.hpp" +// #include "physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp" + +namespace geos +{ +namespace constitutive +{ +class FrictionBase; +} + +namespace LUsolver{ + +// ============================================================================ +// Minimal dense linear algebra: LU with partial pivoting +// ============================================================================ + +// In-place LU decomposition of an n x n matrix stored row-major in A +// (A is overwritten: U in upper incl. diag, L multipliers below diag, +// unit diagonal for L implied). piv[i] = original row that ended up at row i. +// Returns false if the matrix is (numerically) singular. +bool luDecompose(std::vector& A, int n, std::vector& piv); + + +// Solve A x = b given the LU-factored A (from luDecompose) and its pivot. +std::vector luSolve(const std::vector& LU, int n, + const std::vector& piv, + const std::vector& b); + +// Convenience: solve a SMALL fixed-size system (used for the 3x3 bubble +// condensation) via the same LU machinery, working on a flat copy. +template +std::array,N> invertSmall(const std::array,N>& M); + +} + +namespace mimicAssembly{ + +using IndexType = std::ptrdiff_t; +using CRSMat = CRSMatrix< real64, IndexType, IndexType >; + +array2d buildD(double K, double G); + +// ---------------------------------------------------------------------------- +// Reference hex: trilinear shape functions and their parent-coordinate +// derivatives at (xi, eta, zeta) in [-1,1]^3. Standard node ordering. +// ---------------------------------------------------------------------------- +constexpr std::array,8> refCorners = {{ + {-1,-1,-1},{ 1,-1,-1},{ 1, 1,-1},{-1, 1,-1}, + {-1,-1, 1},{ 1,-1, 1},{ 1, 1, 1},{-1, 1, 1} +}}; + +void shapeAndGrad(double xi, double eta, double zeta, + std::array& N, + std::array,8>& dNdXi); + +// Face bubble (Frigo Eq. 20a) for the face at xi_j = sign, on an axis-aligned +// hex. faceAxis in {0,1,2} = which parent coordinate is pinned to the face; +// faceSign = +1 or -1. +void bubbleAndGrad(double xi, double eta, double zeta, + int faceAxis, double faceSign, + double& b, array1d& dbdXi); + +// 2x2x2 Gauss rule +const double gp = 1.0 / std::sqrt(3.0); +const std::array,8> gaussPts = {{ + {-gp,-gp,-gp},{ gp,-gp,-gp},{ gp, gp,-gp},{-gp, gp,-gp}, + {-gp,-gp, gp},{ gp,-gp, gp},{ gp, gp, gp},{-gp, gp, gp} +}}; +const double gaussW = 1.0; // weight 1 each for 2-pt rule per axis + +// Voigt B-column for a generic scalar shape function's physical gradient +array2d bColumn(double dNdx, double dNdy, double dNdz); +array2d computeJacobian(const std::array,8>& coords, + const std::array,8>& dNdXi); +// ---------------------------------------------------------------------------- +// Element geometry: axis-aligned unit cube, origin = (x0,y0,z0). +// Physical gradient = parent gradient * 2/L (L=1 here) -> factor 2. +// detJ = (L/2)^3 = 1/8. +// ---------------------------------------------------------------------------- +struct HexElem +{ + // double x0, y0, z0; // corner with min x,y,z + std::array,8> coords; // PHYSICAL corner coordinates (arbitrary, + // not assumed axis-aligned or unit-length) + + std::array node; // global node ids + bool hasBubble = false; + int bubbleId = -1; // global bubble block index + int bubbleFaceAxis = 0; // which parent axis the fault face sits on + double bubbleFaceSign = 1.0; +}; + +// Local nodal-nodal (24x24), nodal-bubble (24x3), bubble-bubble (3x3) +struct LocalBlocks +{ + static constexpr int nnodes = 24; + static constexpr int ncenters = 3; + + array2d Kuu{nnodes,nnodes}; + array2d Kub{nnodes,ncenters}; + array2d Kbb{ncenters,ncenters}; +}; + +LocalBlocks integrateElement(const array2d& D, const HexElem& e); +struct ContactState { double tN = 0.0; bool open = true; double theta = 0.; double phi = 0.; }; +void setElement(const ContactState& contact, HexElem& elemA, HexElem& elemB, array2d& RR ); +// Kcond = Kuu - Kub * Kbb^-1 * Kub^T (3x3 inverse via our own LU, small enough +// that we go through the same generic solver rather than a closed form) +array2d condense(const LocalBlocks& L); + +CRSMat assemble(const array2d& D, const ContactState& contact); + +struct SolverStepResult +{ + SolverStepResult(std::vector& u_,double gN_, int newtonIterations_, bool converged_): + u(u_),gN(gN_),newtonIterations(newtonIterations_),converged(converged_){}; + std::vector u; + double gN = 0.; + int newtonIterations = 0; + bool converged = false; +}; + +std::vector solveStep(double imposedUx, ContactState& contact, CRSMat const & Kelastic, + double epsN, int kMaxNewton, std::function& updateNormalTraction); + + + +} + + + + +class FrictionDriver : public ConstitutiveDriver +{ +public: + FrictionDriver( const string & name, + Group * const parent ); + + static string catalogName() + { return "FrictionDriver"; } + + void postInputInitialization() override; + + bool execute() override; + + void getColumnNames( string_array & columnNames ) const override; + + template< typename FRICTION_TYPE > + void + runTest( FRICTION_TYPE & friction, + const arrayView2d< real64, 1 > & table ); + +private: + /** + * @brief Get the friction model from the catalog + */ + constitutive::FrictionBase & getFriction(); + constitutive::FrictionBase const & getFriction() const; + + void initializeTable(); + + /** + * @struct viewKeyStruct holds char strings and viewKeys for fast lookup + */ + struct viewKeyStruct : ConstitutiveDriver::viewKeyStruct + { + constexpr static char const * frictionNameString() + { return "friction"; } + + constexpr static char const * contactNameString() + { return "contact"; } + + constexpr static char const * displacementFunctionString() + { return "dispControl"; } + + // constexpr static char const * dJumpFunctionString() + // { return "dJumpControl"; } + + constexpr static char const * stressFunctionRString() + { return "stressControlsR"; } + + constexpr static char const * stressFunctionLString() + { return "stressControlsL"; } + + constexpr static char const * thetaString() + { return "yTiltAngle";} + + constexpr static char const * phiString() + { return "zTiltAngle";} + + constexpr static char const * normalDispTolFac() + { return "tolJumpN"; } + + constexpr static char const * normalTractionTolFac() + { return "tolNormalTrac"; } + + constexpr static char const * slidingTolFac() + { return "tolJumpT"; } + + constexpr static char const * iterPenNFac() + { return "iterPenNFac"; } + + constexpr static char const * iterPenTFac() + { return "iterPenTFac"; } + + //geometry + constexpr static char const * area() + { return "faceArea"; } + + constexpr static char const * volume() + { return "neighborsVolume"; } + + constexpr static char const * shear() + { return "neighborsShear"; } + + constexpr static char const * bulk() + { return "neighborsBulk"; } + + constexpr static char const * simultaneous() + { return "simultaneous"; } + + constexpr static char const * maxNewtonIterString() + { return "maxNewtonIter"; } + + }; + + // Time is defined in base class +// enum columnKeys { NTRAC=1, STRAC0, STRAC1, NDJUMP, DSLIP0, DSLIP1, CC, FS, +// NEWTRAC, SNEWTRAC0, SNEWTRAC1, +// NJUMP, SLIP0, SLIP1, +// NTOL, TTOL, NTRACTOL, +// ITERPEN0, ITERPEN1, TLIM, +// ITER }; +enum columnKeys { NTRAC=1, STRAC0, STRAC1, DISP, NDJUMP, DSLIP0, DSLIP1, CC, FS, + NEWTRAC, SNEWTRAC0, SNEWTRAC1, + NJUMP, SLIP0, SLIP1, + NTOL, TTOL, NTRACTOL, + ITERPEN0, ITERPEN1, + FEMNTRAC, FEMGN, FEMNEWTONITER, FEMCONVERGED, // NEW + TLIM, + ITER }; + + string m_dispFunctionName; ///< + // string m_dJumpFunctionName; ///< + string m_stressFunctionsNamesR; ///< + string m_stressFunctionsNamesL; ///< + + real64 m_theta{0.0}; ///< x-tilt of fault + real64 m_phi{0.0}; ///< y-tilt of fault + + real64 m_normalDispTolFac{1.e-8}; + real64 m_normalTracTolFac{100.}; + real64 m_slidingTolFac{1e-5}; + + real64 m_iterPenNFac{100}; + real64 m_iterPenTFac{1}; + + //geometry for iterPen + real64 m_area{0.0}; + array1d< real64 > m_volume{}; + array1d< real64 > m_shearModulus{}; + array1d< real64 > m_bulkModulus{}; + + integer m_isSimultaneous{1}; + integer m_maxNewtonIter{20}; ///< Maximum Newton iterations + + string m_frictionName; ///< frictionType identifier +}; + +} + +#endif //GEOS_CONSTITUTIVEDRIVERS_CONTACT_FRICTIONDRIVER_HPP diff --git a/src/coreComponents/constitutiveDrivers/contact/FrictionDriverRunTest.cpp b/src/coreComponents/constitutiveDrivers/contact/FrictionDriverRunTest.cpp new file mode 100644 index 00000000000..f06fe66436d --- /dev/null +++ b/src/coreComponents/constitutiveDrivers/contact/FrictionDriverRunTest.cpp @@ -0,0 +1,27 @@ +#include "FrictionDriverRunTest.hpp" +#include "constitutive/contact/CoulombFriction.hpp" +#include "constitutive/contact/FrictionlessContact.hpp" +#include "constitutive/contact/RateAndStateFriction.hpp" +#include + + +namespace geos +{ + +template +void +FrictionDriver::runTest( constitutive::CoulombFriction &, const arrayView2d< real64 > & ); + +template +void +FrictionDriver::runTest( constitutive::FrictionlessContact &, const arrayView2d< real64 > & ); + +template +void +FrictionDriver::runTest( constitutive::RateAndStateFriction< std::integral_constant< bool, true > > &, const arrayView2d< real64 > & ); + +template +void +FrictionDriver::runTest( constitutive::RateAndStateFriction< std::integral_constant< bool, false > > &, const arrayView2d< real64 > & ); + +} diff --git a/src/coreComponents/constitutiveDrivers/contact/FrictionDriverRunTest.hpp b/src/coreComponents/constitutiveDrivers/contact/FrictionDriverRunTest.hpp new file mode 100644 index 00000000000..be8237c4b62 --- /dev/null +++ b/src/coreComponents/constitutiveDrivers/contact/FrictionDriverRunTest.hpp @@ -0,0 +1,243 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC* + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +#ifndef GEOS_FRICTIONDRIVERRUNTEST_HPP_ +#define GEOS_FRICTIONDRIVERRUNTEST_HPP_ + +#include "constitutiveDrivers/contact/FrictionDriver.hpp" +#include "physicsSolvers/solidMechanics/contact/FractureState.hpp" +#include "physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp" +#include "constitutive/solid/SolidFields.hpp" + +#include +#include + +namespace geos +{ + +template< typename FRICTION_TYPE > +void +FrictionDriver::runTest( FRICTION_TYPE & friction, + const arrayView2d< real64 > & table ) +{ + + array1d< integer > const ghostRank( 1 ); ghostRank[0] = -1; + + array1d< integer > fractureState( 1 ); + + fractureState[0] = fields::contact::FractureState::Stick; + array2d< real64 > traction( 1, 3 ); + array2d< real64 > jump( 1, 3 ); + array2d< real64 > djump( 1, 3 ); + + bool isSimultaneous = m_isSimultaneous; + real64 const slidingCheckTol = .05; //default + real64 cos_phi = cos( m_phi * M_PI/180 ); + real64 cos_theta = cos( m_theta * M_PI/180 ); + real64 sin_phi = sin( m_phi * M_PI/180 ); + real64 sin_theta = sin( m_theta * M_PI/180 ); + + real64 normalDispTolFac = m_normalDispTolFac; + real64 normalTracTolFac = m_normalTracTolFac; + real64 slidingTolFac = m_slidingTolFac; + real64 iterPenNFac = m_iterPenNFac; + real64 iterPenTFac = m_iterPenTFac; + + real64 area = m_area; + real64 volumes[2]{m_volume[0], m_volume[1]}; + real64 shearModuli[2]{m_shearModulus[0], m_shearModulus[1]}; + real64 bulkModuli[2]{m_bulkModulus[0], m_bulkModulus[1]}; + + // mimicAssembly only models a single homogeneous material for its two-hex + // mock geometry -- average the two neighboring cells' moduli. + array2d const D = mimicAssembly::buildD( + 0.5*(bulkModuli[0]+bulkModuli[1]), + 0.5*(shearModuli[0]+shearModuli[1]) ); + + // Persistent contact state carried across rows. theta/phi are fixed + // geometry inputs, so Kelastic only needs to be built once. + mimicAssembly::ContactState contact{ + table( 0, NTRAC ), // initial normal traction guess + true, // start unconstrained + m_theta * M_PI/180.0, + m_phi * M_PI/180.0 }; + + mimicAssembly::CRSMat const Kelastic = mimicAssembly::assemble( D, contact ); + + real64 cumulativeImposedUx = 0.0; // running sum of table(ei,NDJUMP) + + //TODO computeTolerance eleme to Elem + typename FRICTION_TYPE::KernelWrapper const kernelWrapper = friction.createKernelUpdates(); + + integer const numRows = m_table.size( 0 ); + // forAll< parallelDevicePolicy<> >( numRows, + // [&friction, &table, &kernelWrapper, + // &ghostRank, + // // &normalDisplacementTol, &normalTractionTol, &slidingTol, + // &normalDispTolFac, &normalTracTolFac, &slidingTolFac, + // &iterPenNFac, &iterPenTFac, + // &area, &volumes, &bulkModuli, &shearModuli, + // &cos_phi, &sin_phi, &cos_theta, &sin_theta, + // &slidingCheckTol, &isSimultaneous, + // &jump, &djump, + // &fractureState, &traction ] + // GEOS_HOST_DEVICE ( integer const ei ) + + for(integer ei=0; ei rotationMatrix( 3, 3 ); + rotationMatrix[0][0] = cos_phi*cos_theta; rotationMatrix[0][1] = -sin_phi; rotationMatrix[0][2] = cos_phi*sin_theta; + rotationMatrix[1][0] = sin_phi*sin_theta; rotationMatrix[1][1] = cos_phi; rotationMatrix[1][2] = sin_theta*sin_phi; + rotationMatrix[2][0] = -sin_theta; rotationMatrix[2][1] = 0.; rotationMatrix[2][2] = cos_theta; + + SolidMechanicsAugmentedLagrangianContact::computeTolerancePerFace( area, + volumes, + bulkModuli, + shearModuli, + rotationMatrix, + normalDispTolFac, + normalTracTolFac, + slidingTolFac, + iterPenNFac, + iterPenTFac, + normalDispTol, + slidingTol, + normalTractionTol, + iterativePen ); + + //small adapter + array2d< real64 > const iterPen( 1, 5 ); + iterPen[0][0] = iterativePen[0]; + iterPen[0][1] = iterativePen[1]; + iterPen[0][2] = iterativePen[1]; iterPen[0][3] = iterativePen[1]; iterPen[0][4] = 0.; + + array1d< real64 > const a_normalDisplacementTol( 1 ); a_normalDisplacementTol[0]=normalDispTol;//normalDispTol should scale as 1/E + array1d< real64 > const a_normalTractionTol( 1 ); a_normalTractionTol[0]=normalTractionTol; + array1d< real64 > const a_slidingTol( 1 ); a_slidingTol[0]=slidingTol; + + + + + auto [newTraction, condCov] = SolidMechanicsAugmentedLagrangianContact::updateTractionAndConstraintCheck( 1, + friction, + isSimultaneous, + slidingCheckTol, + a_normalDisplacementTol, + a_normalTractionTol, + a_slidingTol, + iterPen, + jump, + djump, + ghostRank, + fractureState.toView(), + traction.toView() + ); + + kernelWrapper.updateFractureState( 0, + jump[0], + newTraction[0], + fractureState[0] ); + + + + + cumulativeImposedUx += table( ei, DISP ); + + // Bridges solveStep's Newton iterate (tNold, gN) to the actual configured + // friction model, reusing the same call as the primary update above. Built + // fresh each row: captures this row's tolerance/iterPen locals by reference. + std::function< mimicAssembly::ContactState(double,double,double) > updateNormalTraction = + [&]( double tNold, double gN, double /*epsNArg*/ ) -> mimicAssembly::ContactState + { + jump[0][0] = gN; jump[0][1] = 0.0; jump[0][2] = 0.0; + traction[0][0] = tNold; traction[0][1] = table( ei, STRAC0 ); traction[0][2] = table( ei, STRAC1 ); + + kernelWrapper.updateFractureState( 0, jump[0], traction[0], fractureState[0] ); + + std::tie(newTraction, condCov) = SolidMechanicsAugmentedLagrangianContact::updateTractionAndConstraintCheck( + 1, friction, isSimultaneous, slidingCheckTol, + a_normalDisplacementTol, a_normalTractionTol, a_slidingTol, + iterPen, jump, djump, ghostRank, fractureState.toView(), traction.toView() ); + + kernelWrapper.updateFractureState( 0, jump[0], newTraction[0], fractureState[0] ); + + return mimicAssembly::ContactState{ + newTraction[0][0], + fractureState[0] == fields::contact::FractureState::Open, + contact.theta, contact.phi }; + }; + + // epsN: reuse this row's normal iterative penalty (already folds in + // area/volume/moduli scaling via computeTolerancePerFace) so the FEM + // Newton loop's contact tangent stays consistent with the table-driven path. + const auto stepResult = + mimicAssembly::solveStep( cumulativeImposedUx, contact, Kelastic, iterativePen[0], m_maxNewtonIter, updateNormalTraction ); + + int iter = 0; + for(const auto& step : stepResult ){ + table( ei*m_maxNewtonIter + iter, FEMNTRAC ) = contact.tN; + table( ei*m_maxNewtonIter + iter, FEMGN ) = step.gN; + table( ei*m_maxNewtonIter + iter, FEMNEWTONITER ) = step.newtonIterations; + table( ei*m_maxNewtonIter + iter, FEMCONVERGED ) = step.converged ? 1.0 : 0.0; + + //old cols + table( ei*m_maxNewtonIter + iter, CC ) = condCov[0]; + table( ei*m_maxNewtonIter + iter, FS ) = fractureState[0]; + table( ei*m_maxNewtonIter + iter, NEWTRAC ) = newTraction[0][0]; + table( ei*m_maxNewtonIter + iter, SNEWTRAC0 ) = newTraction[0][1]; + table( ei*m_maxNewtonIter + iter, SNEWTRAC1 ) = newTraction[0][2]; + + + table( ei*m_maxNewtonIter + iter, NTOL ) = normalDispTol; + table( ei*m_maxNewtonIter + iter, TTOL ) = slidingTol; + table( ei*m_maxNewtonIter + iter, NTRACTOL ) = normalTractionTol; + + table( ei*m_maxNewtonIter + iter, ITERPEN0 ) = iterPen[0][0]; + table( ei*m_maxNewtonIter + iter, ITERPEN1 ) = iterPen[0][1]; + ++iter; + } + // } ); //TODO restore once forAll<> +} +}//function +}//namespace + + +#endif //GEOS_FRICTIONDRIVERRUNTEST_HPP_ diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp index a928891a726..a7f648e114e 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp @@ -41,7 +41,9 @@ #include "finiteElement/FiniteElementDiscretization.hpp" #include "mesh/DomainPartition.hpp" +#include #include +#include #if defined( GEOS_USE_CUDA ) #include @@ -1067,6 +1069,84 @@ void SolidMechanicsAugmentedLagrangianContact::updateState( DomainPartition & do GEOS_UNUSED_VAR( domain ); } + +std::tuple< array2d< real64 >, array1d< int > > SolidMechanicsAugmentedLagrangianContact::updateTractionAndConstraintCheck( + std::ptrdiff_t const rsize, + FrictionBase const & frictionLaw, + bool isSimultaneous, + real64 const slidingCheckTolerance, + arrayView1d< real64 const > const & normalDisplacementTolerance, + arrayView1d< real64 const > const & normalTractionTolerance, + arrayView1d< real64 const > const & slidingTolerance, + arrayView2d< real64 const > const & iterativePenalty, + arrayView2d< real64 const > const & dispJump, + arrayView2d< real64 const > const & deltaDispJump, + arrayView1d< integer const > const & ghostRank, + arrayView1d< integer const > const & fractureState, + arrayView2d< real64 > const & traction ) +{ + array2d< real64 > traction_new; + std::ptrdiff_t const sizes[2] = {rsize, 3}; + traction_new.resize( 2, sizes ); + array1d< int > condConv; + condConv.resize( rsize ); + + // Update the traction field based on the displacement results from the nonlinear solve + constitutiveUpdatePassThru( frictionLaw, [&] ( auto & castedFrictionLaw ) + { + using FrictionType = TYPEOFREF( castedFrictionLaw ); + typename FrictionType::KernelWrapper frictionWrapper = castedFrictionLaw.createKernelUpdates(); + + if( isSimultaneous ) + { + solidMechanicsALMKernels::ComputeTractionSimultaneousKernel:: + launch< parallelDevicePolicy<> >( rsize, + iterativePenalty, + traction, + dispJump, + deltaDispJump, + traction_new.toView() ); + } + else + { + solidMechanicsALMKernels::ComputeTractionKernel:: + launch< parallelDevicePolicy<> >( rsize, + frictionWrapper, + iterativePenalty, + traction, + dispJump, + deltaDispJump, + traction_new.toView() ); + } + } ); + + // real64 const slidingCheckTolerance = m_slidingCheckTolerance; + + constitutiveUpdatePassThru( frictionLaw, [&] ( auto & castedFrictionLaw ) + { + using FrictionType = TYPEOFREF( castedFrictionLaw ); + typename FrictionType::KernelWrapper frictionWrapper = castedFrictionLaw.createKernelUpdates(); + + solidMechanicsALMKernels::ConstraintCheckKernel:: + launch< parallelDevicePolicy<> >( rsize, + frictionWrapper, + ghostRank, + traction, + dispJump, + deltaDispJump, + normalTractionTolerance, + normalDisplacementTolerance, + slidingTolerance, + slidingCheckTolerance, + fractureState, + condConv.toView() ); + } ); + + + return std::make_tuple( traction_new, condConv ); +} + + bool SolidMechanicsAugmentedLagrangianContact::updateConfiguration( DomainPartition & domain, integer const GEOS_UNUSED_PARAM( configurationLoopIter ) ) { @@ -1107,61 +1187,27 @@ bool SolidMechanicsAugmentedLagrangianContact::updateConfiguration( DomainPartit std::ptrdiff_t const sizes[ 2 ] = {subRegion.size(), 3}; traction_new.resize( 2, sizes ); - arrayView2d< real64 > const traction_new_v = traction_new.toView(); - - condConv.resize( subRegion.size()); + // arrayView2d< real64 > const traction_new_v = traction_new.toView(); + + // condConv.resize( subRegion.size()); + + std::tie( traction_new, condConv ) = updateTractionAndConstraintCheck( + subRegion.size(), + frictionLaw, + m_simultaneous, + m_slidingCheckTolerance, + normalDisplacementTolerance, + normalTractionTolerance, + slidingTolerance, + iterativePenalty, + dispJump, + deltaDispJump, + ghostRank, + fractureState, + traction + ); arrayView1d< int > const condConv_v = condConv.toView(); - // Update the traction field based on the displacement results from the nonlinear solve - constitutiveUpdatePassThru( frictionLaw, [&] ( auto & castedFrictionLaw ) - { - using FrictionType = TYPEOFREF( castedFrictionLaw ); - typename FrictionType::KernelWrapper frictionWrapper = castedFrictionLaw.createKernelUpdates(); - - if( m_simultaneous ) - { - solidMechanicsALMKernels::ComputeTractionSimultaneousKernel:: - launch< parallelDevicePolicy<> >( subRegion.size(), - iterativePenalty, - traction, - dispJump, - deltaDispJump, - traction_new_v ); - } - else - { - solidMechanicsALMKernels::ComputeTractionKernel:: - launch< parallelDevicePolicy<> >( subRegion.size(), - frictionWrapper, - iterativePenalty, - traction, - dispJump, - deltaDispJump, - traction_new_v ); - } - } ); - - real64 const slidingCheckTolerance = m_slidingCheckTolerance; - - constitutiveUpdatePassThru( frictionLaw, [&] ( auto & castedFrictionLaw ) - { - using FrictionType = TYPEOFREF( castedFrictionLaw ); - typename FrictionType::KernelWrapper frictionWrapper = castedFrictionLaw.createKernelUpdates(); - - solidMechanicsALMKernels::ConstraintCheckKernel:: - launch< parallelDevicePolicy<> >( subRegion.size(), - frictionWrapper, - ghostRank, - traction, - dispJump, - deltaDispJump, - normalTractionTolerance, - normalDisplacementTolerance, - slidingTolerance, - slidingCheckTolerance, - fractureState, - condConv_v ); - } ); RAJA::ReduceSum< parallelDeviceReduce, localIndex > localSum[5] = { RAJA::ReduceSum< parallelDeviceReduce, localIndex >( 0 ), @@ -1923,6 +1969,77 @@ void SolidMechanicsAugmentedLagrangianContact::addCouplingSparsityPattern( Domai } +void SolidMechanicsAugmentedLagrangianContact::computeTolerancePerFace( real64 const area, + real64 const (&volume)[2], + real64 const (&bulkModulus)[2], + real64 const (&shearModulus)[2], + arraySlice2d< real64 const > const & faceRotationMatrix, + real64 const tolJumpDispNFac, + real64 const tolJumpDispTFac, + real64 const tolNormalTracFac, + real64 const iterPenaltyNFac, + real64 const iterPenaltyTFac, + real64 & normalDisplacementTolerance, + real64 & slidingTolerance, + real64 & normalTractionTolerance, + real64 (& iterativePenalty)[2] ) +{ + // approximation of the stiffness along coordinate directions + // ( first, second ) index -> ( element index, direction ) + // 1. T -> top (index 0), B -> bottom (index 1) + // 2. the coordinate direction (x, y, z) + real64 stiffDiagApprox[ 2 ][ 3 ]; + real64 averageYoungModulus = 0.0; + real64 averageConstrainedModulus = 0.0; + real64 averageCharLength = 0.0; + + for( localIndex i = 0; i < 2; ++i ) + { + // Get linear elastic isotropic constitutive parameters for the element + real64 const K = bulkModulus[i]; + real64 const G = shearModulus[i]; + real64 const E = 9.0 * K * G / ( 3.0 * K + G ); + real64 const nu = ( 3.0 * K - 2.0 * G ) / ( 2.0 * ( 3.0 * K + G ) ); + real64 const M = K + 4.0 / 3.0 * G; + + real64 const charLength = pow( volume[i], 1.0 / 3.0 ); + + // Combine E and nu to obtain a stiffness approximation (like it was an hexahedron) + for( localIndex j = 0; j < 3; ++j ) + { + stiffDiagApprox[ i ][ j ] = E / ( ( 1.0 + nu )*( 1.0 - 2.0*nu ) ) * 4.0 / 9.0 * ( 2.0 - 3.0 * nu ) * charLength; + } + + averageYoungModulus += 0.5*E; + averageConstrainedModulus += 0.5*M; + averageCharLength += 0.5*charLength; + } + + // Average the stiffness and compute the inverse + real64 invStiffApprox[ 3 ][ 3 ] = { { 0 } }; + for( localIndex j = 0; j < 3; ++j ) + { + invStiffApprox[ j ][ j ] = ( stiffDiagApprox[ 0 ][ j ] + stiffDiagApprox[ 1 ][ j ] ) / ( stiffDiagApprox[ 0 ][ j ] * stiffDiagApprox[ 1 ][ j ] ); + } + + // Rotate in the local reference system, computing R^T * (invK) * R + real64 temp[ 3 ][ 3 ]; + LvArray::tensorOps::Rij_eq_AkiBkj< 3, 3, 3 >( temp, faceRotationMatrix, invStiffApprox ); + real64 rotatedInvStiffApprox[ 3 ][ 3 ]; + LvArray::tensorOps::Rij_eq_AikBkj< 3, 3, 3 >( rotatedInvStiffApprox, temp, faceRotationMatrix ); + LvArray::tensorOps::scale< 3, 3 >( rotatedInvStiffApprox, area ); + + // Finally, compute tolerances and iterative penalties for the given fracture element + normalDisplacementTolerance = rotatedInvStiffApprox[ 0 ][ 0 ] * averageYoungModulus * tolJumpDispNFac; + slidingTolerance = sqrt( pow( rotatedInvStiffApprox[ 1 ][ 1 ], 2 ) + + pow( rotatedInvStiffApprox[ 2 ][ 2 ], 2 ) ) * averageYoungModulus * tolJumpDispTFac; + normalTractionTolerance = tolNormalTracFac * ( averageConstrainedModulus / averageCharLength ) * + normalDisplacementTolerance; + + iterativePenalty[0] = iterPenaltyNFac * averageConstrainedModulus / averageCharLength; + iterativePenalty[1] = iterPenaltyTFac * averageConstrainedModulus / averageCharLength; +} + void SolidMechanicsAugmentedLagrangianContact::computeTolerances( DomainPartition & domain ) const { GEOS_MARK_FUNCTION; @@ -1975,20 +2092,19 @@ void SolidMechanicsAugmentedLagrangianContact::computeTolerances( DomainPartitio arrayView1d< integer const > const ghostRank = subRegion.ghostRank(); + real64 const tolJumpDispNFac = m_tolJumpDispNFac; + real64 const tolJumpDispTFac = m_tolJumpDispTFac; + real64 const tolNormalTracFac = m_tolNormalTracFac; + real64 const iterPenaltyNFac = m_iterPenaltyNFac; + real64 const iterPenaltyTFac = m_iterPenaltyTFac; + forAll< parallelHostPolicy >( subRegion.size(), [=] ( localIndex const kfe ) { - if( ghostRank[kfe] < 0 ) { - real64 const area = faceArea[kfe]; - // approximation of the stiffness along coordinate directions - // ( first, second ) index -> ( element index, direction ) - // 1. T -> top (index 0), B -> bottom (index 1) - // 2. the coordinate direction (x, y, z) - real64 stiffDiagApprox[ 2 ][ 3 ]; - real64 averageYoungModulus = 0.0; - real64 averageConstrainedModulus = 0.0; - real64 averageCharLength = 0.0; + real64 volume[2]{}; + real64 bulk[2]{}; + real64 shear[2]{}; for( localIndex i = 0; i < 2; ++i ) { @@ -1997,51 +2113,30 @@ void SolidMechanicsAugmentedLagrangianContact::computeTolerances( DomainPartitio localIndex const esr = faceToElemSubRegion[faceIndex][0]; localIndex const ei = faceToElemIndex[faceIndex][0]; - real64 const volume = elemVolume[er][esr][ei]; - - // Get linear elastic isotropic constitutive parameters for the element - real64 const K = bulkModulus[er][esr][ei]; - real64 const G = shearModulus[er][esr][ei]; - real64 const E = 9.0 * K * G / ( 3.0 * K + G ); - real64 const nu = ( 3.0 * K - 2.0 * G ) / ( 2.0 * ( 3.0 * K + G ) ); - real64 const M = K + 4.0 / 3.0 * G; - - real64 const charLength = pow( volume, 1.0 / 3.0 ); - - // Combine E and nu to obtain a stiffness approximation (like it was an hexahedron) - for( localIndex j = 0; j < 3; ++j ) - { - stiffDiagApprox[ i ][ j ] = E / ( ( 1.0 + nu )*( 1.0 - 2.0*nu ) ) * 4.0 / 9.0 * ( 2.0 - 3.0 * nu ) * charLength; - } - - averageYoungModulus += 0.5*E; - averageConstrainedModulus += 0.5*M; - averageCharLength += 0.5*charLength; - } - - // Average the stiffness and compute the inverse - real64 invStiffApprox[ 3 ][ 3 ] = { { 0 } }; - for( localIndex j = 0; j < 3; ++j ) - { - invStiffApprox[ j ][ j ] = ( stiffDiagApprox[ 0 ][ j ] + stiffDiagApprox[ 1 ][ j ] ) / ( stiffDiagApprox[ 0 ][ j ] * stiffDiagApprox[ 1 ][ j ] ); + volume[i] = elemVolume[er][esr][ei]; + bulk[i] = bulkModulus[er][esr][ei]; + shear[i] = shearModulus[er][esr][ei]; } - // Rotate in the local reference system, computing R^T * (invK) * R - real64 temp[ 3 ][ 3 ]; - LvArray::tensorOps::Rij_eq_AkiBkj< 3, 3, 3 >( temp, faceRotationMatrix[ kfe ], invStiffApprox ); - real64 rotatedInvStiffApprox[ 3 ][ 3 ]; - LvArray::tensorOps::Rij_eq_AikBkj< 3, 3, 3 >( rotatedInvStiffApprox, temp, faceRotationMatrix[ kfe ] ); - LvArray::tensorOps::scale< 3, 3 >( rotatedInvStiffApprox, area ); - - // Finally, compute tolerances for the given fracture element - normalDisplacementTolerance[kfe] = rotatedInvStiffApprox[ 0 ][ 0 ] * averageYoungModulus * m_tolJumpDispNFac; - slidingTolerance[kfe] = sqrt( pow( rotatedInvStiffApprox[ 1 ][ 1 ], 2 ) + - pow( rotatedInvStiffApprox[ 2 ][ 2 ], 2 )) * averageYoungModulus * m_tolJumpDispTFac; - normalTractionTolerance[kfe] = m_tolNormalTracFac * (averageConstrainedModulus / averageCharLength) * - (normalDisplacementTolerance[kfe]); - - iterativePenalty[kfe][0] = m_iterPenaltyNFac*averageConstrainedModulus/(averageCharLength); - iterativePenalty[kfe][1] = m_iterPenaltyTFac*averageConstrainedModulus/(averageCharLength); + real64 iterativePenaltyKfe[2]{}; + + computeTolerancePerFace( faceArea[kfe], + volume, + bulk, + shear, + faceRotationMatrix[kfe], + tolJumpDispNFac, + tolJumpDispTFac, + tolNormalTracFac, + iterPenaltyNFac, + iterPenaltyTFac, + normalDisplacementTolerance[kfe], + slidingTolerance[kfe], + normalTractionTolerance[kfe], + iterativePenaltyKfe ); + + iterativePenalty[kfe][0] = iterativePenaltyKfe[0]; + iterativePenalty[kfe][1] = iterativePenaltyKfe[1]; } } ); } diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp index 02af0550280..281c16d8734 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp @@ -22,6 +22,7 @@ #define GEOS_PHYSICSSOLVERS_CONTACT_SOLIDMECHANICSAUGMENTEDLAGRANGIANCONTACT_HPP_ #include "physicsSolvers/solidMechanics/contact/ContactSolverBase.hpp" +#include "constitutive/contact/FrictionSelector.hpp" namespace geos { @@ -215,8 +216,62 @@ class SolidMechanicsAugmentedLagrangianContact : public ContactSolverBase */ void createBubbleCellList( DomainPartition & domain ) const; + + //TODO (change return idiom to passed through interface) + static + GEOS_HOST_DEVICE + std::tuple< array2d< real64 >, array1d< int > > updateTractionAndConstraintCheck( std::ptrdiff_t const rsize, + constitutive::FrictionBase const & frictionLaw, + bool isSimultaneous, + real64 const slidingCheckTolerance, + arrayView1d< real64 const > const & normalDisplacementTolerance, + arrayView1d< real64 const > const & normalTractionTolerance, + arrayView1d< real64 const > const & slidingTolerance, + arrayView2d< real64 const > const & iterativePenalty, + arrayView2d< real64 const > const & dispJump, + arrayView2d< real64 const > const & deltaDispJump, + arrayView1d< integer const > const & ghostRank, + arrayView1d< integer const > const & fractureState, + arrayView2d< real64 > const & traction ); + + + /** + * @brief Compute the augmented-Lagrangian tolerances and iterative penalties for a single + * fracture (contact interface) element. + * @param area area of the fracture element + * @param volume volume of the two adjacent cells (index 0/1 -> top/bottom) + * @param bulkModulus bulk modulus of the two adjacent cells + * @param shearModulus shear modulus of the two adjacent cells + * @param faceRotationMatrix rotation matrix bringing global -> local (fault) reference frame + * @param tolJumpDispNFac factor for the normal displacement jump tolerance + * @param tolJumpDispTFac factor for the tangential displacement jump (sliding) tolerance + * @param tolNormalTracFac factor for the normal traction tolerance + * @param iterPenaltyNFac factor for the normal iterative penalty + * @param iterPenaltyTFac factor for the tangential iterative penalty + * @param[out] normalDisplacementTolerance normal displacement jump tolerance for this element + * @param[out] slidingTolerance sliding (tangential displacement jump) tolerance for this element + * @param[out] normalTractionTolerance normal traction tolerance for this element + * @param[out] iterativePenalty iterative penalties (normal, tangential) for this element + */ + static GEOS_HOST_DEVICE void computeTolerancePerFace( real64 const area, + real64 const (&volume)[2], + real64 const (&bulkModulus)[2], + real64 const (&shearModulus)[2], + arraySlice2d< real64 const > const &faceRotationMatrix, + real64 const tolJumpDispNFac, + real64 const tolJumpDispTFac, + real64 const tolNormalTracFac, + real64 const iterPenaltyNFac, + real64 const iterPenaltyTFac, + real64 & normalDisplacementTolerance, + real64 & slidingTolerance, + real64 & normalTractionTolerance, + real64 ( &iterativePenalty )[2] ); + private: + + /** * @brief Validate that tetrahedral meshes use high-order quadrature rules * @param meshBodies the group containing the mesh bodies diff --git a/src/coreComponents/schema/schema.xsd b/src/coreComponents/schema/schema.xsd index d678b835b53..84cd6e459ca 100644 --- a/src/coreComponents/schema/schema.xsd +++ b/src/coreComponents/schema/schema.xsd @@ -565,6 +565,10 @@ + + + + @@ -6513,6 +6517,7 @@ See note on referenceReservoirRegion for reservoir condition options--> + @@ -6627,6 +6632,29 @@ Information output from lower logLevels is added with the desired log level + + + + + + + + + + + + + + + + + + + + + + + +