Self.formulation.fields

PyLith
1

Hi,

In the main time loop, after a certain number of time steps, I am trying to switch the formulation from Implicit to Explicit. In lieu of that, I am adding two fields i.e. disp(t-dt) and acceleration(t) to self.formulation.fields. Now I am aware that I am not changing self.formulation from Implicit to Explicit, but I am transfering self.formulation.fields directly as an argument wherever I need the Explicit computations. For some reason, the code does not recognize the appended variables disp(t-dt) and acceleration(t). This is what I am doing

175 fields.add(“disp(t-dt)”, “displacement”)
176 fields.add(“acceleration(t)”, “acceleration”)
177
178 dispTmdt = fields.get(“disp(t-dt)”)
179 dispTmdt.zero()
180
181 lengthScale = normalizer.lengthScale()
182 timeScale = normalizer.timeScale()
183 velocityScale = lengthScale / timeScale
184 accelerationScale = velocityScale / timeScale
185 accelerationT = fields.get(“acceleration(t)”)
186 accelerationT.scale(accelerationScale.value)
187 accelerationT.zero()
188

Why does the code not recognize these appended variables and do the Explicit computations as if we started with the Explicit formulation?

Saumik

2

It would be helpful to have more information. Can you describe in broader terms what you are trying to do? What version of the code are you working with? Can you point me to a Git branch that has these changes?

3

I start the simulation with Implicit formulation, and want to switch to Explicit formulation somewhere down the line in the simulation at the point at which the fault starts to slip. I am using pylith version 1.8.0. After I sent the previous message, I added those lines of code to def initialize in /problems/Implicit.py to ensure the solution fields have disp(t-dt) and acceleration(t) right off the bat even though we start with Implicit.

The problem I have is now, when I do self.updateSettings inside _reformjacobian, I get the following error:

Fatal error. Calling MPI_Abort() to abort PyLith application.
Traceback (most recent call last):
File “/home/ec2-user/pylithTrial2/pylith/lib/python2.7/site-packages/pylith/apps/PetscApplication.py”, line 65, in onComputeNodes
self.main(*args, **kwds)
File “/home/ec2-user/pylithTrial2/pylith/lib/python2.7/site-packages/pylith/apps/PyLithApp.py”, line 126, in main
self.problem.run(self)
File “/home/ec2-user/pylithTrial2/pylith/lib/python2.7/site-packages/pylith/problems/TimeDependent.py”, line 215, in run
exp_obj.hybrid_step(t, dt, self.normalizer, self.dimension, self.formulation.fields, self.formulation.gprs, self.mesh, self.materials, self.bc, self.interfaces)
File “/home/ec2-user/pylithTrial2/pylith/lib/python2.7/site-packages/pylith/problems/Explicit.py”, line 320, in hybrid_step
self.hybrid_prestep(t, dt, gprs, fields, mesh, materials, boundaryConditions, interfaceConditions)
File “/home/ec2-user/pylithTrial2/pylith/lib/python2.7/site-packages/pylith/problems/Explicit.py”, line 240, in hybrid_prestep
self.hybrid_reformJacobian(self.jacobian, t, dt, gprs, fields)
File “/home/ec2-user/pylithTrial2/pylith/lib/python2.7/site-packages/pylith/problems/Explicit.py”, line 665, in hybrid_reformJacobian
self.updateSettings(jacobian, fields, gprs, t, dt) # BJ: add self.gprs
File “/home/ec2-user/pylithTrial2/pylith/lib/python2.7/site-packages/pylith/problems/problems.py”, line 88, in updateSettings
def updateSettings(self, *args): return _problems.Formulation_updateSettings(self, *args)
NotImplementedError: Wrong number or type of arguments for overloaded function ‘Formulation_updateSettings’.
Possible C/C++ prototypes are:
updateSettings(pylith::problems::Formulation *,pylith::topology::Jacobian *,pylith::topology::SolutionFields *,PylithScalar const,PylithScalar const)
updateSettings(pylith::problems::Formulation *,pylith::topology::Jacobian *,pylith::topology::SolutionFields *,pylith::problems::GPRS *,double const,double const)
updateSettings(pylith::problems::Formulation *,pylith::topology::Field< pylith::topology::Mesh > *,pylith::topology::SolutionFields *,PylithScalar const,PylithScalar const)
updateSettings(pylith::problems::Formulation *,pylith::topology::Field< pylith::topology::Mesh > *,pylith::topology::SolutionFields *,pylith::problems::GPRS *,double const,double const)

which is strange, given that I have given the right data structures as arguments. The code framework has coupling with flow, and unfortunately will not be able to share it here because we have not yet made it open source :expressionless:

Saumik

4

Make sure the SWIG interface file is up to date with the C++ header file. If the SWIG interface has any differences with the public C++ interface definition, then the cast to the C++ type won’t work and SWIG will say the arguments don’t match.

5

SWIG interface is up to date with the C++ header file. The updatesettings does work when I running Implicit initially with the following routines in /problems/Formulation.py

595 def _reformJacobian(self, t, dt):
596 “”"
597 Reform Jacobian matrix for operator.
598 “”"
599 from pylith.mpi.Communicator import mpi_comm_world
600 comm = mpi_comm_world()
601
602 self._debug.log(resourceUsageString())
603 if 0 == comm.rank:
604 self._info.log(“Integrating Jacobian operator.”)
605 self._eventLogger.stagePush(“Reform Jacobian”)
606 if self.gprs == None:
607 self.updateSettings(self.jacobian, self.fields, t, dt)
608 else:
609 self.updateSettings(self.jacobian, self.fields, self.gprs, t, dt) # BJ: add self.gprs
610 self.reformJacobian()
611
612 self._eventLogger.stagePop()
613
614 if self.viewJacobian:
615 from pylith.mpi.Communicator import Communicator
616 comm = Communicator(self.mesh.comm())
617 self.jacobianViewer.view(self.jacobian, t, comm)
618
619 self._debug.log(resourceUsageString())
620 return
621
622 def _reformResidual(self, t, dt):
623 “”"
624 Reform residual vector for operator.
625 “”"
626 from pylith.mpi.Communicator import mpi_comm_world
627 comm = mpi_comm_world()
628
629 if 0 == comm.rank:
630 self._info.log(“Integrating residual term in operator.”)
631 self._eventLogger.stagePush(“Reform Residual”)
632 if self.gprs == None:
633 self.updateSettings(self.jacobian, self.fields, t, dt)
634 else:
635 self.updateSettings(self.jacobian, self.fields, self.gprs, t, dt) # BJ: add self.gprs
636 self.reformResidual()
637
638 self._eventLogger.stagePop()
639 self._debug.log(resourceUsageString())
640 return
641
642

But it is not working after switching to Explicit with the following routines in /problems/Explicit.py

595 def _reformJacobian(self, t, dt):
596 “”"
597 Reform Jacobian matrix for operator.
598 “”"
599 from pylith.mpi.Communicator import mpi_comm_world
600 comm = mpi_comm_world()
601
602 self._debug.log(resourceUsageString())
603 if 0 == comm.rank:
604 self._info.log(“Integrating Jacobian operator.”)
605 self._eventLogger.stagePush(“Reform Jacobian”)
606 if self.gprs == None:
607 self.updateSettings(self.jacobian, self.fields, t, dt)
608 else:
609 self.updateSettings(self.jacobian, self.fields, self.gprs, t, dt) # BJ: add self.gprs
610 self.reformJacobian()
611
612 self._eventLogger.stagePop()
613
614 if self.viewJacobian:
615 from pylith.mpi.Communicator import Communicator
616 comm = Communicator(self.mesh.comm())
617 self.jacobianViewer.view(self.jacobian, t, comm)
618
619 self._debug.log(resourceUsageString())
620 return
621
622 def _reformResidual(self, t, dt):
623 “”"
624 Reform residual vector for operator.
625 “”"
626 from pylith.mpi.Communicator import mpi_comm_world
627 comm = mpi_comm_world()
628
629 if 0 == comm.rank:
630 self._info.log(“Integrating residual term in operator.”)
631 self._eventLogger.stagePush(“Reform Residual”)
632 if self.gprs == None:
633 self.updateSettings(self.jacobian, self.fields, t, dt)
634 else:
635 self.updateSettings(self.jacobian, self.fields, self.gprs, t, dt) # BJ: add self.gprs
636 self.reformResidual()
637
638 self._eventLogger.stagePop()
639 self._debug.log(resourceUsageString())
640 return
641
642

That is where I am surprised that the updatesettings is not working after switching to Explicit

Saumik

6

Sorry, I meant the following lines in /problems/Explicit.py

663 def hybrid_reformJacobian(self, jacobian, t, dt, gprs, fields):
664
665 self.updateSettings(jacobian, fields, gprs, t, dt) # BJ: add self.gprs
666 ModuleExplicit.reformJacobianLumped(self)
667
668 return
669
670 def hybrid_reformResidual(self, jacobian, t, dt, gprs, fields):
671
672 self.updateSettings(jacobian, fields, gprs, t, dt) # BJ: add self.gprs
673 self.reformResidual()
674
675 return

Saumik

7

PyLith v1.8.0 was designed to use only a single time-stepping formulation within a given simulation. There are several pieces in the initialization where we assume that the code is called once and only once. it looks to me like updateSettings() is one of them. The _fields data member in the C++ Formulation object is not deallocated prior to assigning a new fields object. I am guessing that other places are holding onto a copy of the original _fields object. I think changing the time stepping formulation in the middle of a simulation would require significant restructuring and reorganization of the top-level code to ensure everything gets initialized and updated properly. Unfortunately, if you want to go down this path, I don’t think we can provide support given that PyLith v1.8.0 is ten years old and the code is outdated. There are many, many changes since then to both PETSc and PyLith. We have significantly restructured the code so that the current development version is much better suited for earthquake cycle modeling (a work in progress using coupled problems, not switching the time stepping within a problem) and coupled physics.