tps1.py

### TPS1 ###
from . import _tpsx
from . import _details


################################################################################


class TPS1(_tpsx.TPSX):#, _details.interfaces.tpsPreconditionerJacobiStationary): TODO
  '''
  implement integrator from
  Alouges, F. A new finite element scheme for Landau–Lifchitz equations. Discrete Contin. Dyn. Syst. Ser. S 1, 2 (2008), 187–196.
  '''


#------------------------------------------------------------------------------#

  
  def __init__(self, geometry, parameters = None, output = None, numthreads = 1):
    
    _tpsx.TPSX.__init__(self, geometry, parameters, output, numthreads)

    self._theta = 0.5


#------------------------------------------------------------------------------#


  def _SetUpAuxiliaryFunctions(self):
    from commics._gridFunctions import Magnetization

    super()._SetUpAuxiliaryFunctions()

    self._MagUnprojected = Magnetization(self._geometry, self._V, self._X \
      , self._Mag, maxwellCoupling="none", solvetol=self._solvetol)


#------------------------------------------------------------------------------#


  def GetReady(self):

    self._solvetol = self.parameters.TOL
    self._theta = self.parameters.theta

    _tpsx.TPSX.GetReady(self)
    self._Mag.Project()


#------------------------------------------------------------------------------#


  def _AssembleRHS(self):
    self._Assemble_f()


#------------------------------------------------------------------------------#


  def _SetRHS_3D(self):
    self._RHS_3D = self._f.vec.FV().NumPy()


#------------------------------------------------------------------------------#


  def _StepUpdate(self):
    self._Mag.m.vec.data += self._k * self._Mag.v.gf.vec

    if self.output.traceUnprojectedExchangeEnergy:
      self._MagUnprojected.SetFromOther(self._Mag, withStrayField=False)

    self._Mag.Project()


#------------------------------------------------------------------------------#


  def _SetUpSystems(self):
    from commics._tools import ScipySparseFromNgMatrix
    from ngsolve import TaskManager, TET, IntegrationRule, SymbolicLFI \
      , LinearForm, SymbolicBFI, BilinearForm

    #blf
    lhs_stat_cf = self._MassTerm() + self._StiffTerm()
    lhs_instat_cf = self._SkewTerm()
    
    self._a_stat = BilinearForm(self._X)
    self._a_stat += SymbolicBFI( lhs_stat_cf.Compile(realcompile=True, wait=True) )
    
    with TaskManager():
      self._a_stat.Assemble()
    self._A_stat = ScipySparseFromNgMatrix(self._a_stat)

    ir3 = IntegrationRule(TET, order=3)
    
    self._a_instat = BilinearForm(self._X)
    self._a_instat += SymbolicBFI( lhs_instat_cf.Compile(realcompile=True, wait=True), intrule=ir3 )
    
    with TaskManager():
      self._a_instat.Assemble()
    self._A_instat = ScipySparseFromNgMatrix(self._a_instat) 

    #lf
    rhs_cf = self._ExchangeTerm()
    rhs_cf += self._LinearTerms_pi()
    rhs_cf += self._NonLinearTerms_Pi()
    rhs_cf += self._ExplicitTerms()

    self._f = LinearForm(self._X)
    self._f += SymbolicLFI(rhs_cf.Compile(realcompile=True, wait=True))


#------------------------------------------------------------------------------#


  def _Compute_v(self):
    import scipy.sparse.linalg
    
    self._solution, succ = scipy.sparse.linalg.gmres(self._LHS, self._RHS \
      , x0=self._guess, tol=self._solvetol, maxiter=4000 \
      , M=self._preconditioner, callback=self._callback)

    self._Mag.v.gf.vec.FV().NumPy()[:] = self._Q.transpose().dot(self._solution)

    self._UpdateGuess()


#------------------------------------------------------------------------------#


  def _Assemble_f(self):
    self._f.Assemble()


#------------------------------------------------------------------------------#

  
  def _MassTerm(self):
    return self.parameters.alpha * self._Mass()


#------------------------------------------------------------------------------#


  def _StiffTerm(self):
    if(self._theta == 0.0):
      return self._ZeroCF()
    else:
      return self._C_ex * self._theta * self._k * self._Stiff()


#------------------------------------------------------------------------------#


  @staticmethod
  def _IsMagProjected():
    return True


#------------------------------------------------------------------------------#


  def SetTheta(self, value):
    self.parameters.theta = value


#------------------------------------------------------------------------------#


  def TraceUnprojectedExchangeEnergy(self, value):
    self.output.traceUnprojectedExchangeEnergy = value


#------------------------------------------------------------------------------#


################################################################################