[CIG-MC] MC Teleconference notes from May 15, 2006

[Ariel Shoresh]

MC Telecon 5-15-06
 
Thorsten Becker 
Shijie Zhong
Carolina Lithgow-Bertelloni
Clint Conrad
Bernhard Steinberger
Michael Gurnis
 
 
Becker opened the meeting, referring to email of topics he¹d like to
discuss. He wanted to talk briefly on what the community has done, and then
discuss the plan for going forward.
 
In terms of benchmarks, they have put together a document with integrations
following the AGU meeting. Bernhard used his Fortran code to run the
benchmarks, and then they used the results to make few revisions to the
benchmarks. So far they are Density and viscosity models, and it made
predictions for the test cases. He set this up to run automatically for
himself, and completed his revision of Bernhard¹s code in C, so now its more
modular and can eventually be integrated with a Citcom type code. It¹s
called HC. It¹s basic but working, and it was able to replicate the velocity
benchmarks of the test case, and it only had a problem with one of the
velocity benchmarks. But they should match, so that¹s why more work needs to
be done. Zhong did a test with an analytic result, and both codes have
confirmed the result.
 
This HC code is designed for a modular approach, as it would be easy to
replace routines in code, and expands easily to take on other functions. It
has instantaneous flow for compression without boundary. The HC code is
currently available in the CIG subversion for download, and he¹d appreciate
if people would take a look at the code
(http://www.geodynamics.org/websvn/). It would be good to have more
benchmarks designed and to make the decision on whether or not HC should be
developed to form the basis for a community code that can be integrated with
Citcom. 
 
Is this something the community needs, or is it too narrow? It won¹t need to
be a formal decision; the group may just want to work on benchmarks instead.
The group also should discuss what is wanted in this community code, since
it would be slightly duplicating things like Bernhard¹s code.
 
The call moved on to discussing benchmark results.
 
Lithgow-Bertelloni hadn¹t done the benchmarks, as her code doesn¹t blow up
to that degree. Becker said that there was one high-density case that he
would like to run through her code, and she agreed to that. Becker suggested
that Bernhard might run them, as he has fixes already designed. His team
tried to remove the properties, but he would be interested in how her code
treats the benchmarks.
Bernhard was also curious as to why it won¹t blow up with LB¹s code, but
with theirs. He wondered if it was because her code doesn¹t match the
boundary the way their code does. LB responded that she had gone up to
degree 50 with no problems.
Bernhard decided that they should wait to see how her code runs the
benchmark before they think of editing anything.
 
LB said that Rick was going to put together something mathematical to show
what the problem was, and Becker said that he was also planning to write a
small code for it. Craig O¹Neil has looked at it, and reaffirmed his
commitment to do some coding, but Becker agreed that the group should see
how LB¹s code handles the benchmarks first.
 
LB asked if Becker had the plate algorithm written into his code, and
responded that he had plate velocity, and ability to set up slip. The input
files have spherical harmonics, and you can download the conventions from
the web. The output should be given spatially for spherical harmonic
overkill, as if it works once, it will work all the time.
 
LB asked if the conventions were specified in the PDF, and Becker responded
that it¹s just using a standard theoretical physics convention, and there¹s
a conversion factor that goes to -1. Becker thought that if the benchmark
will work with LB¹s code, and match Tromp¹s solution and the two numerical
solutions, then the group can reasonably assume things are going in the
right direction
 
Zhong thought this code was important not only for research, but for
teaching. Others agreed. While he has a code, it¹s a variant.  LB said that
Steve Grand and other seismologists wouldn¹t mind having something that
works with their tomographic inversions.
 
Becker said that even though the community has codes that do all that, it
would be good to have a C version of this, as it would be easier for users
to understand. He suggested everyone look at code to see if they can
understand what it can do. He and O¹Neil will hold off on further
development until we get more feedback.
 
Question was raised about replacing the propagator with a numerical
solution. Conrad asked if that would make the code run more slowly, but
Becker said no, and that using a different propagator or solver would help.
O¹Neil was interested in the schemes used. He was also thinking of adding
his contributions back into the code, which LB agreed was a good idea.
 
CIG has an initiative to work on a framework for this. Zhong said that it¹s
the numerical analysis for compressibility, and that hasn¹t gotten too far
yet, but there have been some 2-D tests, and those were okay. LB said that
she was thinking about other propagator codes, such as one from Hagar¹s
student, or Yanik¹s code.
 
Zhong said the community really wants a 1-D code, but we also want to
implement compressibility in a 1-D code without gravity to compare it to the
numerical code. Some of the compressible code calculations include things
that make it more complicated. Bernhard¹s code has the geoid implemented.
Bernhard said that in his implementation, he followed Svetlana¹s code. He
didn¹t have a dynamic pressure term; he programmed her equations directly
into his version of the code. He wasn¹t sure if that was the right method,
and perhaps similar codes could be looked at to see what equations they use.
It is more complicated in incompressible versions, both in equation writing
and in coding. It wasn¹t something he had thought about.

Becker asked if the group agreed that in terms of priority, the geoid comes
first, and then compressibility is a research issue that can be done later.
All agreed. LB asked if the benchmarks would provide information like
density, velocity fields, etc., and Becker said that he would upload his
data and input files to his benchmark website. In terms of seismology, there
is an old benchmark website that has models in a format similar to his code,
so they can be used to test as well. There is an effort in SPICE to update
their repository, and trying to get people to use these models.
 
Zhong asked if the code would include plate motion, like LB¹s code, but
Becker said that it wouldn¹t be in the first release, but he could add it to
the list for future versions. He currently does these equations externally,
puts the data into a table, and feeds it into the code. LB said that her
code has this included, and the inversion is done afterwards, but all the
other variables/calculations are done in the code, like unit rotation. The
final inversion is at the end.
 
Becker said that the HC code is set up to use the same subroutines, but
writing a different main program around it. There¹s one routine called
solve, and you might want to call this into a code that does torque, or into
Citcom for velocity at edges. So following up on a design idea the HC code
would have single function, then loop into another code to do things like
torque rotations. It would require a bit more coding in that plan.
 
LB said keeping track of torque direction and which plate is affected is
necessary. Becker agreed, but that right now he uses an external Fortran
program so that would have to incorporated in C. LB said her code has this
already, and she has maps to put in grids to use the code. But if the user
wanted to do specific research, they would have to develop their own input
file. Becker said he does the same for torque for his code.
 
Priorities for code then are geoid, then plate motions. LB said these are
what most people are interested in.
 
Becker said that in terms of time dedicated to the project, he¹d be happy to
add geoid, aspect documentations, then run this code against Bernhard¹s. The
plate motions would take a bit longer ­ a couple of months perhaps ­ and
he¹d be happy to work on this. There is also O¹Neil¹s idea, but he thought
the priorities are set as they stand.
 
LB said that since this code can interface with Citcom, it may be more
effective for doing the convection versions O¹Neil mentioned. Bernhard said
that Alex Forte made a convection version, and that he¹s used it. The user
would only need to add a few variables to make it compressible, but I don¹t
think it¹s a priority for the community either, as there are a lot of codes
that do this. Others agreed.
 
Becker then encouraged everyone to check the code out, and that he¹d add a
subdirectory to his subversion to add benchmarks. He will send an email on
how to access these files. Zhong asked how soon could HC be put out as a
release. Becker said that was too ³intimate² a question. Currently the code
is function for velocity and benchmarks, but the spherical harmonic density
has a problem, like Bernhard¹s code. Once he gets the geoid coded in, he
would be happy to release it as an alpha. He doesn¹t have any documentation
to speak of, but since its coded in C, pretty much anyone can figure it out.
LB agreed that releasing it without the geoid wasn¹t useful. Zhong was
hoping it would also include plate motions, which Conrad agreed, but that
working on plate motions in C would take a good week of coding, and would
take time for testing etc. The routines I have are all external, and the
geoid currently looks at Bernhard¹s code and put into C. Bernhard said he
added his own thing from the Hager code, and used the Matlab code equations
from Sveltana. It would be useful to compare it with other geoid codes.
 
LB said she had two version of her code ­ one that does plate motions, and
one that does the geoid. The geoid was from Hager, and does slab motions at
the same time. Her team doesn¹t have an integrated code either.  LB said
that the boundary conditions for the geoid and plates are different, so
that¹s why the distinction needs to put in the documentation. If we have one
code doing both, the user has to know what variables are placed where. If
they do a plate, the velocity would look terrible. Bernhard agreed, saying
that people would think the code is bad.
 
LB thought that before the release we should talk about fundamental issues
in the code, more than just needs of community. Becker didn¹t want to build
in lots of functions, as users should do some thing themselves. It also
should be modular. Bernhard agreed, and said the limits should be
documented. The code is merely a tool.
 
Becker was wondering if help were available, and Zhong said that CIG would
help either financially or man-hours once we get a plan together. O¹Neil
declined funding previously. Becker noted that he was also very busy, which
also concerned Zhong.
Becker would like to get plate motion help for the alpha release, and maybe
LB could share her approach using C.
 
LB didn¹t think it was all that different than Becker¹s. You sum up
tractions on edge of plates and balance them, which is different than
Conrad¹s. Her code is specific to give plate functions. Once you have
tractions from external field and the internal part, you get a plate map to
do the integration on outside of code. She would give the code to Becker, as
³it¹s not a secret².  Becker said that in his code the user has to switch
the traction calculation back on and add some integration routines for the
plates. LB said her code doesn¹t do that. She calculates the curve and
Cartesian stuff beforehand and the map for the torque is given as
information to read in as well. Becker said that was similar to his code.
 
Bernhard offered to send some routines that he has from Rick.  LB thought
that might take too much time, and may be better just to supply maps.  She
did note that her code is in Fortran, as she doesn¹t know C.
 
Becker hoped CIG could help with spatial to spectral conversions, and with
linking it to Citcom. He needed to go from spectral to velocity in space. He
asked if CIG¹s developers had any subroutines that have representational
fields on a sphere.
 
Gurnis said no, but that Luis Armendariz has been tasked with a related
problem, which is to go from one mesh formulation to another, in terms of
building a generalized tool to compare various models¹ output for benchmark
and code comparisons. The idea is for one code does this, so that mantle
convection codes could be compared against each other, and then for the
comparison of long term and short term tectonic codes. CIG plans to use
shared libraries to compute these. We can put this functionality into the
code, as I think there is some value for this. To compare the outputs of HC
to Citcom or Terra, you would need to extrapolate this down into a mesh.
Gurnis suggested Becker meet with Armendariz to hammer out some kind of
concrete set of tasks for the SSC to chew on.

Becker thought it sounded like a good plan. Not in terms of spherical
harmonic conversion, but going from one point on the sphere to another.
 
Conrad offered to help LB on benchmark testing of Becker¹s code. LB said she
would do the geo, and Conrad said he¹d do velocity
 
Becker was happy that there was a plan in place, and that Plate formulation
is something the group has to think about a little more to get the code done
as quickly and robustly as possible. LB offered to iterate with him on this.
 
Becker asked Gurnis if he should send an email to Armendariz, and Gurnis
said that was fine. He said he¹d also mention it to him when he returns to
the office.