diff --git a/bin/direct_comparison_small_strain.py b/bin/direct_comparison_small_strain.py index 8ac568a..12acd0e 100644 --- a/bin/direct_comparison_small_strain.py +++ b/bin/direct_comparison_small_strain.py @@ -1,161 +1,228 @@ ## Most of the following is copied from https://github.com/tdegeus/GooseFFT (MIT license) import numpy as np import scipy.sparse.linalg as sp import itertools import os import sys sys.path.append(os.path.join(os.getcwd(), "language_bindings/python")) import pyMuSpectre as µ # ----------------------------------- GRID ------------------------------------ ndim = 2 # number of dimensions -N = 99 #31 # number of voxels (assumed equal for all directions) +N = 31 #31 # number of voxels (assumed equal for all directions) offset = 3 #9 ndof = ndim**2*N**2 # number of degrees-of-freedom cell = µ.SystemFactory(µ.get_2d_cube(N), µ.get_2d_cube(1.), µ.Formulation.small_strain) # ---------------------- PROJECTION, TENSORS, OPERATIONS ---------------------- # tensor operations/products: np.einsum enables index notation, avoiding loops # e.g. ddot42 performs $C_ij = A_ijkl B_lk$ for the entire grid trans2 = lambda A2 : np.einsum('ij... ->ji... ',A2 ) ddot42 = lambda A4,B2: np.einsum('ijkl...,lk... ->ij... ',A4,B2) ddot44 = lambda A4,B4: np.einsum('ijkl...,lkmn...->ijmn...',A4,B4) dot22 = lambda A2,B2: np.einsum('ij... ,jk... ->ik... ',A2,B2) dot24 = lambda A2,B4: np.einsum('ij... ,jkmn...->ikmn...',A2,B4) dot42 = lambda A4,B2: np.einsum('ijkl...,lm... ->ijkm...',A4,B2) dyad22 = lambda A2,B2: np.einsum('ij... ,kl... ->ijkl...',A2,B2) shape = tuple((N for _ in range(ndim))) # identity tensor [single tensor] i = np.eye(ndim) def expand(arr): new_shape = (np.prod(arr.shape), np.prod(shape)) ret_arr = np.zeros(new_shape) ret_arr[:] = arr.reshape(-1)[:, np.newaxis] return ret_arr.reshape((*arr.shape, *shape)) # identity tensors [grid of tensors] I = expand(i) I4 = expand(np.einsum('il,jk',i,i)) I4rt = expand(np.einsum('ik,jl',i,i)) I4s = (I4+I4rt)/2. II = dyad22(I,I) # projection operator [grid of tensors] # NB can be vectorized (faster, less readable), see: "elasto-plasticity.py" # - support function / look-up list / zero initialize delta = lambda i,j: np.float(i==j) # Dirac delta function freq = np.arange(-(N-1)/2.,+(N+1)/2.) # coordinate axis -> freq. axis Ghat4 = np.zeros([ndim,ndim,ndim,ndim,N,N]) # zero initialize # - compute for i,j,l,m in itertools.product(range(ndim),repeat=4): for x,y in itertools.product(range(N), repeat=ndim): q = np.array([freq[x], freq[y]]) # frequency vector if not q.dot(q) == 0: # zero freq. -> mean Ghat4[i,j,l,m,x,y] = -(q[i]*q[j]*q[l]*q[m])/(q.dot(q))**2+\ (delta(j,l)*q[i]*q[m]+delta(j,m)*q[i]*q[l]+\ delta(i,l)*q[j]*q[m]+delta(i,m)*q[j]*q[l])/(2.*q.dot(q)) # (inverse) Fourier transform (for each tensor component in each direction) fft = lambda x: np.fft.fftshift(np.fft.fftn (np.fft.ifftshift(x),shape)) ifft = lambda x: np.fft.fftshift(np.fft.ifftn(np.fft.ifftshift(x),shape)) # functions for the projection 'G', and the product 'G : K : eps' G = lambda A2 : np.real( ifft( ddot42(Ghat4,fft(A2)) ) ).reshape(-1) K_deps = lambda depsm: ddot42(K4,depsm.reshape(ndim,ndim,N,N)) G_K_deps = lambda depsm: G(K_deps(depsm)) # ------------------- PROBLEM DEFINITION / CONSTITIVE MODEL ------------------- # phase indicator: cubical inclusion of volume fraction (9**3)/(31**3) +E2, E1 = 75e10, 70e9 +poisson = .33 + +hard = µ.material.MaterialHooke2d.make(cell, "hard", + E2, poisson) +soft = µ.material.MaterialHooke2d.make(cell, "soft", + E1, poisson) + +#for pixel in cell: +# if ((pixel[0] >= N-offset) and +# (pixel[1] < offset)): +# hard.add_pixel(pixel) +# else: +# soft.add_pixel(pixel) +# phase = np.zeros(shape); phase[-offset:,:offset,] = 1. phase = np.zeros(shape); phase[0,:] = 1. + +phase = np.zeros(shape); +for i, j in itertools.product(range(N), repeat=ndim): + c = N//2 + if (i-c)**2 + (j-c)**2 < (N//5)**2: + #if j<10: + phase[i,j] = 1. + hard.add_pixel([i,j]) + else: + soft.add_pixel([i,j]) # material parameters + function to convert to grid of scalars param = lambda M0,M1: M0*np.ones(shape)*(1.-phase)+M1*np.ones(shape)*phase # K = param(0.833,8.33) # bulk modulus [grid of scalars] # mu = param(0.386,3.86) # shear modulus [grid of scalars] -E2, E1 = 210e9, 70e9 -poisson = .33 K2, K1 = (E/(3*(1-2*poisson)) for E in (E2, E1)) m2, m1 = (E/(2*(1+poisson)) for E in (E2, E1)) K = param(K1, K2) mu = param(m1, m2) # stiffness tensor [grid of tensors] K4 = K*II+2.*mu*(I4s-1./3.*II) -hard = µ.material.MaterialHooke2d.make(cell, "hard", - E2, poisson) -soft = µ.material.MaterialHooke2d.make(cell, "soft", - E1, poisson) - -for pixel in cell: - if ((pixel[0] >= N-offset) and - (pixel[1] < offset)): - hard.add_pixel(pixel) - else: - soft.add_pixel(pixel) # ----------------------------- NEWTON ITERATIONS ----------------------------- # initialize stress and strain tensor [grid of tensors] sig = np.zeros([ndim,ndim,N,N]) eps = np.zeros([ndim,ndim,N,N]) # set macroscopic loading DE = np.zeros([ndim,ndim,N,N]) DE[0,1] += 0.01 DE[1,0] += 0.01 +delEps0 = DE[:ndim, :ndim, 0, 0] µDE = np.zeros([ndim**2, cell.size()]) cell.evaluate_stress_tangent(µDE); µDE[:,:] = DE[:,:,0,0].reshape([-1, 1]) # initial residual: distribute "DE" over grid using "K4" b = -G_K_deps(DE) -G_K_deps2 = lambda de: cell.directional_stiffness(de) -b2 = -G_K_deps2(µDE) +G_K_deps2 = lambda de: cell.directional_stiffness(de.reshape(µDE.shape)).ravel() +b2 = -G_K_deps2(µDE).ravel() +print("b2.shape = {}".format(b2.shape)) eps += DE En = np.linalg.norm(eps) iiter = 0 # iterate as long as the iterative update does not vanish class accumul(object): def __init__(self): self.counter = 0 - def __call__(self, dummy): + def __call__(self, x): self.counter += 1 + print("at step {}: ||Ax-b||/||b|| = {}".format( + self.counter, + np.linalg.norm(G_K_deps(x)-b)/np.linalg.norm(b))) acc = accumul() +cg_tol = 1e-8 +maxiter = 60 +solver = µ.solvers.SolverCGEigen(cell, cg_tol, maxiter, verbose=True) +solver = µ.solvers.SolverCG(cell, cg_tol, maxiter, verbose=True) +try: + r = µ.solvers.newton_cg(cell, delEps0, solver, 1e-5, verbose=True) +except Exception as err: + print(err) while True: - depsm,_ = sp.cg(tol=1.e-8, + depsm,_ = sp.cg(tol=cg_tol, A = sp.LinearOperator(shape=(ndof,ndof),matvec=G_K_deps,dtype='float'), b = b, callback=acc ) # solve linear system using CG + #depsm2,_ = sp.cg(tol=1.e-8, + # A = sp.LinearOperator(shape=(ndof,ndof),matvec=G_K_deps2,dtype='float'), + # b = b2, + # callback=acc + #) # solve linear system using CG eps += depsm.reshape(ndim,ndim,N,N) # update DOFs (array -> tens.grid) sig = ddot42(K4,eps) # new residual stress b = -G(sig) # convert residual stress to residual print('%10.2e'%(np.max(depsm)/En)) # print residual to the screen + print(np.linalg.norm(depsm)/np.linalg.norm(En)) if np.linalg.norm(depsm)/En<1.e-5 and iiter>0: break # check convergence iiter += 1 print("nb_cg: {0}".format(acc.counter)) import matplotlib.pyplot as plt -plt.pcolormesh(eps[0,0,:,:]) + +s11 = sig[0,0, :,:] +s22 = sig[1,1, :,:] +s21_2 = sig[0,1, :, :]*sig[1,0,:, :] +vonM1 = np.sqrt(.5*((s11-s22)**2) + s11**2 + s22**2 + 6*s21_2) + +#vonM1 = np.sqrt(sig[0, 0, :, :]**2 + sig[1, 1, :, :]**2 - sig[0, 0, :, :] * sig[1, 1, :, :] + 3 * sig[0,1]**2) + +plt.figure() +img = plt.pcolormesh(vonM1)#eps[0,1,:,:]) +plt.title("goose") +plt.colorbar(img) + +try: + print(r.stress.shape) + arr = r.stress.T.reshape(N, N, ndim, ndim) + s11 = arr[:,:,0,0] + s22 = arr[:,:,1,1] + s21_2 = arr[:,:,0,1]*arr[:,:,1,0] + vonM2 = np.sqrt(.5*((s11-s22)**2) + s11**2 + s22**2 + 6*s21_2) + + plt.figure() + plt.title("µSpectre") + img = plt.pcolormesh(vonM2)#eps[0,1,:,:]) + plt.colorbar(img) + print("err = {}".format (np.linalg.norm(vonM1-vonM2))) + print("err2 = {}".format (np.linalg.norm(vonM1-vonM1.T))) + print("err3 = {}".format (np.linalg.norm(vonM2-vonM2.T))) + + plt.figure() + plt.title("diff") + img = plt.pcolormesh(vonM1-vonM2.T)#eps[0,1,:,:]) + plt.colorbar(img) +except Exception as err: + print(err) + + plt.show() diff --git a/src/solver/solvers.cc b/src/solver/solvers.cc index fd90dcf..2e261ee 100644 --- a/src/solver/solvers.cc +++ b/src/solver/solvers.cc @@ -1,376 +1,378 @@ /** * file solvers.cc * * @author Till Junge * * @date 20 Dec 2017 * * @brief implementation of solver functions * * @section LICENSE * * Copyright © 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include "solver/solvers.hh" #include "solver/solver_cg.hh" #include "common/iterators.hh" #include #include #include namespace muSpectre { template std::vector de_geus (SystemBase & sys, const GradIncrements & delFs, SolverBase & solver, Real newton_tol, Dim_t verbose) { using Field_t = typename MaterialBase::StrainField_t; auto solver_fields{std::make_unique>()}; solver_fields->initialise(sys.get_resolutions()); // Corresponds to symbol δF or δε auto & incrF{make_field("δF", *solver_fields)}; // Corresponds to symbol ΔF or Δε auto & DeltaF{make_field("ΔF", *solver_fields)}; // field to store the rhs for cg calculations auto & rhs{make_field("rhs", *solver_fields)}; solver.initialise(); if (solver.get_maxiter() == 0) { solver.set_maxiter(sys.size()*DimM*DimM*10); } size_t count_width{}; const auto form{sys.get_formulation()}; std::string strain_symb{}; if (verbose > 0) { //setup of algorithm 5.2 in Nocedal, Numerical Optimization (p. 111) std::cout << "de Geus-" << solver.name() << " for "; switch (form) { case Formulation::small_strain: { strain_symb = "ε"; std::cout << "small"; break; } case Formulation::finite_strain: { strain_symb = "F"; std::cout << "finite"; break; } default: throw SolverError("unknown formulation"); break; } std::cout << " strain with" << std::endl << "newton_tol = " << newton_tol << ", cg_tol = " << solver.get_tol() << " maxiter = " << solver.get_maxiter() << " and Δ" << strain_symb << " =" <(tup)}; auto && grad{std::get<1>(tup)}; std::cout << "Step " << counter + 1 << ":" << std::endl << grad << std::endl; } count_width = size_t(std::log10(solver.get_maxiter()))+1; } // initialise F = I or ε = 0 auto & F{sys.get_strain()}; switch (form) { case Formulation::finite_strain: { F.get_map() = Matrices::I2(); break; } case Formulation::small_strain: { F.get_map() = Matrices::I2().Zero(); for (const auto & delF: delFs) { if (!check_symmetry(delF)) { throw SolverError("all Δε must be symmetric!"); } } break; } default: throw SolverError("Unknown formulation"); break; } // initialise return value std::vector ret_val{}; // initialise materials constexpr bool need_tangent{true}; sys.initialise_materials(need_tangent); Grad_t previous_grad{Grad_t::Zero()}; for (const auto & delF: delFs) { //incremental loop Real incrNorm{2*newton_tol}, gradNorm{1}; auto convergence_test = [&incrNorm, &gradNorm, &newton_tol] () { return incrNorm/gradNorm <= newton_tol; }; Uint newt_iter{0}; for (; (newt_iter < solver.get_maxiter()) && (!convergence_test() || (newt_iter==1)); ++newt_iter) { // obtain material response auto res_tup{sys.evaluate_stress_tangent(F)}; auto & P{std::get<0>(res_tup)}; auto & K{std::get<1>(res_tup)}; auto tangent_effect = [&sys, &K] (const Field_t & dF, Field_t & dP) { sys.directional_stiffness(K, dF, dP); }; if (newt_iter == 0) { DeltaF.get_map() = -(delF-previous_grad); // neg sign because rhs tangent_effect(DeltaF, rhs); incrF.eigenvec() = solver.solve(rhs.eigenvec(), incrF.eigenvec()); F.eigen() -= DeltaF.eigen(); } else { rhs.eigen() = -P.eigen(); sys.project(rhs); + incrF.eigen() = 0; incrF.eigenvec() = solver.solve(rhs.eigenvec(), incrF.eigenvec()); } F.eigen() += incrF.eigen(); incrNorm = incrF.eigen().matrix().norm(); gradNorm = F.eigen().matrix().norm(); if (verbose>0) { std::cout << "at Newton step " << std::setw(count_width) << newt_iter << ", |δ" << strain_symb << "|/|Δ" << strain_symb << "| = " << std::setw(17) << incrNorm/gradNorm << ", tol = " << newton_tol << std::endl; if (verbose-1>1) { std::cout << "<" << strain_symb << "> =" << std::endl << F.get_map().mean() << std::endl; } } } // update previous gradient previous_grad = delF; ret_val.push_back(OptimizeResult{F.eigen(), sys.get_stress().eigen(), convergence_test(), Int(convergence_test()), "message not yet implemented", newt_iter, solver.get_counter()}); //!store history variables here } return ret_val; } template std::vector de_geus (SystemBase & sys, const GradIncrements& delF0, SolverBase & solver, Real newton_tol, Dim_t verbose); // template typename SystemBase::StrainField_t & // de_geus (SystemBase & sys, const GradIncrements& delF0, // const Real cg_tol, const Real newton_tol, Uint maxiter, // Dim_t verbose); template std::vector de_geus (SystemBase & sys, const GradIncrements& delF0, SolverBase & solver, Real newton_tol, Dim_t verbose); /* ---------------------------------------------------------------------- */ template std::vector newton_cg (SystemBase & sys, const GradIncrements & delFs, SolverBase & solver, Real newton_tol, Dim_t verbose) { using Field_t = typename MaterialBase::StrainField_t; auto solver_fields{std::make_unique>()}; solver_fields->initialise(sys.get_resolutions()); // Corresponds to symbol δF or δε auto & incrF{make_field("δF", *solver_fields)}; // field to store the rhs for cg calculations auto & rhs{make_field("rhs", *solver_fields)}; solver.initialise(); if (solver.get_maxiter() == 0) { solver.set_maxiter(sys.size()*DimM*DimM*10); } size_t count_width{}; const auto form{sys.get_formulation()}; std::string strain_symb{}; if (verbose > 0) { //setup of algorithm 5.2 in Nocedal, Numerical Optimization (p. 111) std::cout << "Newton-" << solver.name() << " for "; switch (form) { case Formulation::small_strain: { strain_symb = "ε"; std::cout << "small"; break; } case Formulation::finite_strain: { strain_symb = "Fy"; std::cout << "finite"; break; } default: throw SolverError("unknown formulation"); break; } std::cout << " strain with" << std::endl << "newton_tol = " << newton_tol << ", cg_tol = " << solver.get_tol() << " maxiter = " << solver.get_maxiter() << " and Δ" << strain_symb << " =" <(tup)}; auto && grad{std::get<1>(tup)}; std::cout << "Step " << counter + 1 << ":" << std::endl << grad << std::endl; } count_width = size_t(std::log10(solver.get_maxiter()))+1; } // initialise F = I or ε = 0 auto & F{sys.get_strain()}; switch (sys.get_formulation()) { case Formulation::finite_strain: { F.get_map() = Matrices::I2(); break; } case Formulation::small_strain: { F.get_map() = Matrices::I2().Zero(); for (const auto & delF: delFs) { if (!check_symmetry(delF)) { throw SolverError("all Δε must be symmetric!"); } } break; } default: throw SolverError("Unknown formulation"); break; } // initialise return value std::vector ret_val{}; // initialise materials constexpr bool need_tangent{true}; sys.initialise_materials(need_tangent); Grad_t previous_grad{Grad_t::Zero()}; for (const auto & delF: delFs) { //incremental loop // apply macroscopic strain increment for (auto && grad: F.get_map()) { grad += delF - previous_grad; } Real incrNorm{2*newton_tol}, gradNorm{1}; auto convergence_test = [&incrNorm, &gradNorm, &newton_tol] () { return incrNorm/gradNorm <= newton_tol; }; Uint newt_iter{0}; for (; newt_iter < solver.get_maxiter() && !convergence_test(); ++newt_iter) { // obtain material response auto res_tup{sys.evaluate_stress_tangent(F)}; auto & P{std::get<0>(res_tup)}; rhs.eigen() = -P.eigen(); sys.project(rhs); + incrF.eigen() = 0; incrF.eigenvec() = solver.solve(rhs.eigenvec(), incrF.eigenvec()); F.eigen() += incrF.eigen(); incrNorm = incrF.eigen().matrix().norm(); gradNorm = F.eigen().matrix().norm(); if (verbose > 0) { std::cout << "at Newton step " << std::setw(count_width) << newt_iter << ", |δ" << strain_symb << "|/|Δ" << strain_symb << "| = " << std::setw(17) << incrNorm/gradNorm << ", tol = " << newton_tol << std::endl; if (verbose-1>1) { std::cout << "<" << strain_symb << "> =" << std::endl << F.get_map().mean() << std::endl; } } } // update previous gradient previous_grad = delF; ret_val.push_back(OptimizeResult{F.eigen(), sys.get_stress().eigen(), convergence_test(), Int(convergence_test()), "message not yet implemented", newt_iter, solver.get_counter()}); //store history variables here } return ret_val; } template std::vector newton_cg (SystemBase & sys, const GradIncrements& delF0, SolverBase & solver, Real newton_tol, Dim_t verbose); // template typename SystemBase::StrainField_t & // newton_cg (SystemBase & sys, const GradIncrements& delF0, // const Real cg_tol, const Real newton_tol, Uint maxiter, // Dim_t verbose); template std::vector newton_cg (SystemBase & sys, const GradIncrements& delF0, SolverBase & solver, Real newton_tol, Dim_t verbose); /* ---------------------------------------------------------------------- */ bool check_symmetry(const Eigen::Ref& eps, Real rel_tol){ return rel_tol >= (eps-eps.transpose()).matrix().norm()/eps.matrix().norm(); } } // muSpectre