Theory questions

[kerim371]

Hi,

I have some little experience of writing finite-difference schemes and I used to solve 2D wave eq by composing finite difference matrix and solving the equation using Matlab's mldivide() function. I also read the documentation for provided by JUDI autors at wiki SEG, JUDI docs. This helps a lot, thank you, but few things I can't resolve.

Some of the examples mention Linear and Nonlinear modeling. Could you please explain what is the difference? For now I think that linear modeling uses linear wave equation and nonlinear modeling uses some variation of it and takes into account nonlinear effects. Most likely my understanding is wrong.

Also just to make sure: does JUDI (or Devito i think) solves wave equation in time domain only or it is also possible to solve equations in frequency domain (Helmholtz equation)?

[mloubout]

Hi, thanks for the interest. Concerning your two questions:

• What we call linear and non-linear modelling refer to the inversion operators. The non-linear modeling is the standard wave-equation solve i.e for the acoustic case m * u.dt2 - u.laplace = q solved for u. The linear modeling is the linearization of this equation around m for a perturbation dm , it is also usually called Born modelling. In this case is is actually linear because you can also see it as the action of the Jacobian J*dm. So to summarize, if you are interested in modelling the solution of the wave equation you would solve the non-linear case.

• Yes we solve the equation in the time-domain via explicit time-stepping. We do not have an Helmoltz solver at the time which requires an implicit solver. We however provide some methods for imaging that mimic the Helmoltz equation via computing "on-the-fly DFT" to compute the frequency modes from the wave propagation. You can check the related paper here

[kerim371]

Thank you for explaining this,

JUDI Options contains space_order parameter wich is documented as:
finite difference space order for wave equation (default is 8, needs to be multiple of 4)

Is space_order equal to the number of points that are used to approximate finite difference schema? I.e. as I remember we may get the finite-difference formulas via Taylor series. Thus does the space_order is defined as the number of points in XZ (or XYZ) space that are used to approximate u.laplace at x_i, y_i, z_i?

One more thing.
Does the solution of wave equation (I guess I mean non-linear case) is given by composing the matrix that approximates u.laplace for a shot (for each x,y,z in this shot) and then we solve the pretty big system of linear equations to get the solution at some t_i0 at each x, y, z. To get the soluetion at t_i1 we simply solve the same system of equations with intial condition equal to the previous solution u(t_i0) (plus source Ricker signal I think). And so on.
Or we solve the wave equation without composing matrix u.laplace for the whole xyz space domain but rather iteratively point by point approximate u.laplace and get the solution for some specified x_i, y_i, z_i and t_i.
I'm asking this mostly because I expect that my source and receiver positions are given in double (float64) type while the velocity model is given in float (float32). So if u.laplace is approximated by matrix at each xyz (for a given shot) then probably this matrix will have type double (in my case) and I will get excessive RAM consumption. In this case I should better convert double coordinates to float32 before making any computations. But if RAM consumed by JUDI (Devito) is weakly depend of source/rec coordinates then I would prefer to leave them as double.

[mloubout]

Is space_order equal to the number of points that are used to approximate finite difference schema? I

Not exactly. This defines the order of the FD scheme. In practice this leads to a space_order+1 wide stencil so for example with space_order=4 it will use 5 points (two on each side of the center)

Concerning your second point, JUDI uses Devito to solve the equations. Devito is a stencil DSL, meaning it will not manipulate any matrix but generate stencil code that basically loop over all the gridpoint and apply the stencil update (see here for more details).
As for the data type, JUDI does all the computation in Float32 and will convert your input into it if provided in other data type.

JUDI is a linear-algebra abstraction but we do not work with any explicit matrices. The whole framework is based on JOLI and Matalb's SPOT before that to create matrix-free linear operator.

[kerim371]

Thank you,

I've seen these seismic examples in Devito and I need to make choice whether to use JUDI or directly Devito. These are the questions that are important:

1. Concerning wave modeling in JUDI, lets suppose I want to get snapshots (for displaying right after each is calculated or writing them in custom fileformat) for a given shot and the best case would be to get snapshots directly by making wavefiled calculations in loop. Is there a way to achieve that? Though I tried Options::save_data_to_disk but it I'm getting error (I will report later in issues, some examples doesn't work).
2. I remember I've read (probably from a paper in Geophysics) that JUDI is tested on big data 1.6 TB but I have not found the PC (or cluster) parameters (at least used cores/RAM). For me it would be good to have some assumtion how much RAM/cores I need to calculate wafield in 3D space for a one shot model say 100x100x100 points. Would it be enough my 10-15 GB of RAM. For now I keep in mind that in one shot model weigh 100 MB, then each core should have probably space 120-150 MB to use JUDI.
   Also am I right that JUDI uses outermost parallelization that distributes sources over cores (each core calculates its source)? If there is more detailed info about JUDI parallelization (and GPU computing) please share it.

[mloubout]

Hi

• Concerning 1. you can run the modeling without receivers, for example see wavefield example. The save_data_to_disk option is to save the receiver data (pointwise measurements) into a SEGY file.

• The large scale 3D example was run on a cluster we do not have access to anymore with slurm. The data is 1.6TB but this is the full data survey (measurement), not the wavefield. We used about 500 nodes with 128Gb of RAM per node for this case study. In your case for a 100^3 model you would need 100^3 * nt * 4 Gb of RAM (model size X number of computational time steps X Float32. So you can estimate it for your own case before hand. JUDI will parallelize over the resources you use for julia distributed` whether it is core or nodes. This will be managed by your own cluster manager

• Concerning the actual compute parallelism, Devito generates openmp code that is run from JUDI (openacc on GPU). The GPU computing side is based on offloading meaning the wave equation will run on the GPU (with the Devito flags from the README) ut the rest of JUDI will run on the host CPU.

[kerim371]

@mloubout hi,

I'm looking for a way to store calculated shots to the custom format (based on hdf5 but somehow structured).
For now I do that in loop (by iterating over shots) like in the snippet below but I guess in this case I won't be able to use JUDI distributed computing:

# iterate over shots
for i in 1:nsrc
    # srcGeometry is setup for a single shot 'i'
    # recGeometry is configured in accordance with the given shot 'i'

    # Setup operators
    Pr = judiProjection(info, recGeometry)
    F = judiModeling(info, model; options=opt)
    Ps = judiProjection(info, srcGeometry)

    # Nonlinear modeling
    dobs = Pr*F*adjoint(Ps)*q
  
    # now I need to write the calculated data to the custom format
end

Another way is to write the data to the SEGY and then translate it to custom format.
I thought maybe you have something to propose..

[mloubout]

Yes so in this case you have to do it for each shot. Currently we support in the distributed case the SEGY format, see for example this 3D example:

JUDI.jl/examples/software_paper/examples/generate_data_overthrust_3D.jl

Line 70 in f65e339

file_path = "path/to/shot/records", # replace w/ path to directory in which data will be stored (~total of 1.7 TB)

This write to disk each shot independently when computed. If you wish to work with HDF5 rather than SEGY, then you will need to basically extend the current out-of-core behavior of judiVector to write to HDF5 rather than SEGY. In more detail this function:

JUDI.jl/src/TimeModeling/Types/judiVector.jl

Line 472 in f65e339

function write_shot_record(srcGeometry, srcData, recGeometry, recData,options)

would need to be extended to accept an hdf5 file name an a write to hdf5 function.

[kerim371]

@mloubout hi,

Is it possible to request a status of a computation? For example if the user has big data to compute say for forward modeling and I would like to show a GUI window with progress bar where I could show the number of processed shots.