[CIG-SEISMO] Roundoff inconsistency in assembling routines

[Dimitri Komatitsch]

Dear Surendra,

Thank you very much. This is very interesting to know (in our 
implementation we always assume that a + b + c = c + b + a, even though 
on a computer as you say that is not true because of roundoff).

I agree that we should merge your contribution into the main SVN code of 
SPECFEM3D (assuming your modification has no negative impact on load 
balancing; could you confirm that?). I cc Joseph Charles, could you 
contact him directly to tell him how to merge your modifications (no 
need to cc the whole mailing list then I guess; just cc Daniel and I).

PS (also for Jean-Paul): what about merging the whole "3D/FAULT_SOURCE" 
branch into "3D/SPECFEM3D" to avoid having to maintain two different 
codes? A few months ago Tarje Nissen told me that you were about to do 
that, which would be great (if technically feasible; please let us know).

Thank you,
Best wishes,

Dimitri.

On 10/19/2012 05:07 AM, surendra at caltech.edu wrote:
> Dear Dimitri,
>
> We have noticed that different processors produce results with different
> round-off errors during the force assembly in SPECFEM3D. This apparently
> minor inconvenience actually lead to significant inconsistencies in our
> dynamic rupture simulations if parallel interfaces cut through faults.
> The non-linearity of the dynamic rupture problem likely amplified the
> impact of round-off inconsistencies, beyond what would be noticeable in
> a linear wave propagation simulation. However, as the issue might
> cumulate through the iterations in tomography and source inversions with
> SPECFEM3D, we thought it would be useful to share with you and other
> users what we found and how we fixed it.
>
> Attached are three figures to illustrate how the issue appears in our
> rupture simulation results. All are views of the fault plane:
> 1) fault_processor_partition.png: each color indicates a different processor
> 2) sliprate_snapshot_glitches.png: snapshot of slip velocity using the
> original assemble routines of SPECFEM3D
> 3) sliprate_snapshot_afterFixing_assemblingRoutine.png: snapshot of slip
> velocity using our modified assemble routine
> Comparing carefully figures (2) and (3), there are some glitches in
> figure (2) that disappeared after using the new assembling routine. One
> of the most perceivable glitches is close to (x,z)=(5,-5). Remarkably,
> all glitches are located at intersections between three processors as
> can be seen figure (1). With time, these glitches spread over
> neighboring nodes and pollute the solution.
>
> Here is the origin of the issue. In the assembling routines of SPECFEM3D
> each processor adds the contribution of its neighbors in the order of
> their rank. However, because that sum is added to the local processor
> value, the result is that the overall order in which the contribution
> from each processor (including self) is added is not the same on all
> processors. For instance, at a triple junction between processors 1, 2
> and 3, the sum is
> v1+v2+v3 in proc 1, but
> v2+v1+v3 in proc 2, and
> v3+v1+v2 in proc 3.
> Adding floating point numbers in different orders leads to different
> round-off errors.
>
> We fixed this issue with a few lines of code to insure that the order of
> sums is the same on all processors. We created a new routine
> "assemble_MPI_vector_ext_mesh_w_ordered" in assemble_MPI_vector.f90
> file. The file can be found on the SVN server from the following link:
> http://geodynamics.org/wsvn/cig/seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/assemble_MPI_vector.f90?op=file&rev=0&sc=0
> <http://geodynamics.org/wsvn/cig/seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/assemble_MPI_vector.f90?op=file&rev=0&sc=0>
>
> We suspect this version improves the accuracy and consistency of the
> code without affecting its performance. It might be worth using it in
> the main branch of SPECFEM3D too.
>
> Best regards,
>
> Surendra Somala and Jean-Paul Ampuero - Caltech
>
>
>
> _______________________________________________
> CIG-SEISMO mailing list
> CIG-SEISMO at geodynamics.org
> http://geodynamics.org/cgi-bin/mailman/listinfo/cig-seismo

-- 
Dimitri Komatitsch
CNRS Research Director (DR CNRS), Laboratory of Mechanics and Acoustics,
UPR 7051, Marseille, France    http://komatitsch.free.fr