diff --git a/src/SiconosSolver.py b/src/SiconosSolver.py index cbd863b..9cfb6e8 100644 --- a/src/SiconosSolver.py +++ b/src/SiconosSolver.py @@ -58,13 +58,13 @@ def __call__(self, problem, reactions, velocities, global_velocities= None): if self._mpi_comm is not None: N.NM_MPI_set_comm(problem.M, self._mpi_comm) - print ('SET MPI:', problem, problem.M, self._mpi_comm) + #print('SET MPI:', problem, problem.M, self._mpi_comm) if self._mumps_id is not None: print('MUMPS_id:', self._mumps_id) N.NM_MUMPS_set_id(problem.M, self._mumps_id) - print ('SET MUMPS id:', problem, problem.M, self._mumps_id) - print ('MUMPS id icntl 14:', N.NM_MUMPS_icntl(problem.M, 14)) + #print ('SET MUMPS id:', problem, problem.M, self._mumps_id) + #print ('MUMPS id icntl 14:', N.NM_MUMPS_icntl(problem.M, 14)) # N.frictionContact_display(problem) real_time = c_float() @@ -139,3 +139,112 @@ def guess(self, filename): class SiconosWrappedSolver(SiconosSolver): def __call__(self, problem, reactions, velocities): return self._API(problem, reactions, velocities, self._SO) + + + +class DummySolverOptions(): + def __init__(self): + self.iparam = np.empty(10).astype(int) + self.dparam = np.empty(10) + +def ACfun(reactions, velocities, mu, rho): + n = len(mu) + r = reactions.reshape((n, 3)) + v = velocities.reshape((n, 3)) + (f, A, B) = numerics.fc3d_AlartCurnierFunction(reactions, velocities, mu, rho) + return (f, A, B) + +def NMfun(reactions, velocities, mu, rho): + n = len(mu) + r = reactions.reshape((n, 3)) + v = velocities.reshape((n, 3)) + (f, A, B) = numerics.fc3d_NaturalMapFunction(reactions, velocities, mu, rho) + return (f, A, B) + + +def FBfun(reactions, velocities, mu, rho): + n = len(mu) + r = reactions.reshape((n, 3)) + v = velocities.reshape((n, 3)) + (f, A, B) = numerics.fc3d_FischerBurmeisterFunction(reactions, velocities, mu, rho) + return (f, A, B) + +def AC2(Xfun,W,q,mu): + nc=len(mu) + rho = np.full((nc,3),1.) + #print(W.shape) + def ACi(r): + return Xfun(r, W.dot(r)+q, mu, rho) + + def ACf(r): + return ACi(r)[0].flatten() + + def ACj(r): + Ablks = ACi(r)[1].reshape(nc,9) + Bblks = ACi(r)[2].reshape(nc,9) + A = scipy.sparse.block_diag([Ablks[k,:].reshape(3,3) for k in range(0,nc)]) + B = scipy.sparse.block_diag([Bblks[k,:].reshape(3,3) for k in range(0,nc)]) + D=scipy.sparse.csc_matrix(scipy.sparse.diags(np.random.choice([1-1e-3,1+1e-3], 3*len(mu)))) + J = W*A+B + return scipy.sparse.csc_matrix(J) + return (ACf, ACj) + +#from module import symbol + +from scipy.optimize.nonlin import BroydenFirst, DiagBroyden, InverseJacobian, Anderson, LinearMixing, ExcitingMixing + +import scipy + +import siconos.numerics as numerics + + +def krylov_solv(fcp, reactions, velocities, solver_options): + + tolerance = solver_options.dparam[0] + rdiff = solver_options.dparam[2] + W = numerics.NM_triplet(fcp.M) + mu = fcp.mu + q = fcp.q + it = 0 + convergence = False + func = AC2(ACfun, W, q, mu)[0] + + if solver_options.iparam[9] == 0: + kjac = None + elif solver_options.iparam[9] == 1: + jac = BroydenFirst() + kjac = InverseJacobian(jac) + elif solver_options.iparam[9] == 2: + jac = DiagBroyden() + kjac = InverseJacobian(jac) + elif solver_options.iparam[9] == 3: + jac = Anderson() + kjac = InverseJacobian(jac) + + try: + while ((not convergence) and (it < solver_options.iparam[0])): + + res1 = scipy.optimize.root(fun=func, x0=reactions, method='krylov', + options={'jac_options': {'rdiff': rdiff, + 'method': 'lgmres', + 'inner_M': kjac}, + 'fatol': 1e-10, 'disp': 0, + 'maxiter': 1, + 'line_search': 'wolfe'}) + reactions = res1.x + ce1 = numerics.fc3d_compute_error(fcp, reactions, velocities, + tolerance, None, + np.linalg.norm(fcp.q)) + error = ce1[1] + convergence = (error <= tolerance) + it += 1 + except Exception: + it = 0 + error = np.inf + convergence = False + solver_options.iparam[1] = it + solver_options.dparam[1] = error + if convergence: + return 0 + else: + return 1 diff --git a/src/comp.py b/src/comp.py index 8b02bc7..c751495 100755 --- a/src/comp.py +++ b/src/comp.py @@ -228,8 +228,11 @@ def __call__(self, tpl): print('Exception in internal call') - traceback.print_exc(e) - + try: + traceback.print_exc(e) + except: + print(e) + try: os.remove(output_filename) except: @@ -286,7 +289,7 @@ def __call__(self, tpl): print (list_print, file=report_file) -# @timeout(utimeout) + @timeout(utimeout) def _internal_call(self, solver, problem, filename, pfilename, output_filename): #print("_internal_call") @@ -465,7 +468,10 @@ def pffff(r, v, e): attrs.create('nc', numberOfDegreeofFreedomContacts(filename)) attrs.create('nds', numberOfDegreeofFreedom(filename)) attrs.create('cond_nc', cond_problem(filename)) - attrs.create('digest', digest) + try: + attrs.create('digest', digest) + except: + pass attrs.create('info', info) attrs.create('iter', iter) attrs.create('err', err) diff --git a/src/faf_solvers.py b/src/faf_solvers.py index fbe6fa7..3b63487 100644 --- a/src/faf_solvers.py +++ b/src/faf_solvers.py @@ -1450,6 +1450,60 @@ def fc3d_nsn_ac_eg(problem, reactions, velocities, _SO): admm_asym_sbr.SolverOptions().iparam[N.SICONOS_FRICTION_3D_ADMM_IPARAM_RHO_STRATEGY]= N.SICONOS_FRICTION_3D_ADMM_RHO_STRATEGY_SCALED_RESIDUAL_BALANCING admm_asym_sbr.SolverOptions().iparam[N.SICONOS_FRICTION_3D_ADMM_IPARAM_SYMMETRY]=N.SICONOS_FRICTION_3D_ADMM_FORCED_ASYMMETRY + ############################################## + # krylov + ############################################## + + krylov_s0 = SiconosSolver(name="KRYLOV-BASIC-0", + gnuplot_name="KRYLOV-BASIC-0", + TAG=N.SICONOS_FRICTION_3D_NSN_AC, + iparam_iter=1, + dparam_err=1, + maxiter=self._maxiter, + precision=self._precision, + with_guess=self._with_guess, + API=krylov_solv) + krylov_s0.SolverOptions().iparam[9] = 0 # preconditioning : None + krylov_s0.SolverOptions().dparam[2] = 1e-10 + + krylov_s1 = SiconosSolver(name="KRYLOV-BASIC-1", + gnuplot_name="KRYLOV-BASIC-1", + TAG=N.SICONOS_FRICTION_3D_NSN_AC, + iparam_iter=1, + dparam_err=1, + maxiter=self._maxiter, + precision=self._precision, + with_guess=self._with_guess, + API=krylov_solv) + krylov_s1.SolverOptions().iparam[9] = 1 # preconditioning : BroydenFirst + krylov_s1.SolverOptions().dparam[2] = 1e-10 + + krylov_s2 = SiconosSolver(name="KRYLOV-BASIC-2", + gnuplot_name="KRYLOV-BASIC-2", + TAG=N.SICONOS_FRICTION_3D_NSN_AC, + iparam_iter=1, + dparam_err=1, + maxiter=self._maxiter, + precision=self._precision, + with_guess=self._with_guess, + API=krylov_solv) + krylov_s2.SolverOptions().iparam[9] = 2 # preconditioning : DiagBroyden + krylov_s2.SolverOptions().dparam[2] = 1e-10 + + + krylov_s3 = SiconosSolver(name="KRYLOV-BASIC-3", + gnuplot_name="KRYLOV-BASIC-3", + TAG=N.SICONOS_FRICTION_3D_NSN_AC, + iparam_iter=1, + dparam_err=1, + maxiter=self._maxiter, + precision=self._precision, + with_guess=self._with_guess, + API=krylov_solv) + krylov_s3.SolverOptions().iparam[9] = 3 # preconditioning : Anderson + krylov_s3.SolverOptions().dparam[2] = 1e-10 + + ############################################## # solvers selection @@ -1483,6 +1537,7 @@ def fc3d_nsn_ac_eg(problem, reactions, velocities, _SO): #ProxNSGS, #Proxfixed, regul_series[0], + krylov_s0, krylov_s1, krylov_s2, krylov_s3, #ACLMFixedPoint, #ACLMFixedPoint_vi_fpp, #ACLMFixedPoint_vi_eg]