|
173 | 173 | "from examples.seismic import AcquisitionGeometry\n", |
174 | 174 | "\n", |
175 | 175 | "t0 = 0.\n", |
176 | | - "tn = 1000. \n", |
| 176 | + "tn = 1000.\n", |
177 | 177 | "f0 = 0.010\n", |
178 | 178 | "# First, position source centrally in all dimensions, then set depth\n", |
179 | 179 | "src_coordinates = np.empty((1, 2))\n", |
|
419 | 419 | "outputs": [], |
420 | 420 | "source": [ |
421 | 421 | "# Create FWI gradient kernel \n", |
422 | | - "from devito import Function, TimeFunction, norm\n", |
| 422 | + "from devito import Function, norm\n", |
423 | 423 | "from examples.seismic import Receiver\n", |
424 | 424 | "\n", |
425 | | - "import scipy\n", |
426 | | - "def fwi_gradient(vp_in): \n", |
| 425 | + "def fwi_gradient(vp_in):\n", |
427 | 426 | " # Create symbols to hold the gradient\n", |
428 | 427 | " grad = Function(name=\"grad\", grid=model.grid)\n", |
429 | 428 | " # Create placeholders for the data residual and data\n", |
|
440 | 439 | " for i in range(nshots):\n", |
441 | 440 | " # Update source location\n", |
442 | 441 | " geometry.src_positions[0, :] = source_locations[i, :]\n", |
443 | | - " \n", |
| 442 | + "\n", |
444 | 443 | " # Generate synthetic data from true model\n", |
445 | 444 | " _, _, _ = solver.forward(vp=model.vp, rec=d_obs)\n", |
446 | | - " \n", |
| 445 | + "\n", |
447 | 446 | " # Compute smooth data and full forward wavefield u0\n", |
448 | 447 | " _, u0, _ = solver.forward(vp=vp_in, save=True, rec=d_syn)\n", |
449 | | - " \n", |
| 448 | + "\n", |
450 | 449 | " # Compute gradient from data residual and update objective function \n", |
451 | 450 | " compute_residual(residual, d_obs, d_syn)\n", |
452 | | - " \n", |
| 451 | + "\n", |
453 | 452 | " objective += .5*norm(residual)**2\n", |
454 | 453 | " solver.gradient(rec=residual, u=u0, vp=vp_in, grad=grad)\n", |
455 | | - " \n", |
| 454 | + "\n", |
456 | 455 | " return objective, grad" |
457 | 456 | ] |
458 | 457 | }, |
|
571 | 570 | "name": "stdout", |
572 | 571 | "output_type": "stream", |
573 | 572 | "text": [ |
574 | | - "Objective value is 39292.605833 at iteration 1\n", |
575 | | - "Objective value is 24506.628229 at iteration 2\n", |
576 | | - "Objective value is 14386.573665 at iteration 3\n", |
577 | | - "Objective value is 7907.616850 at iteration 4\n", |
578 | | - "Objective value is 3960.106497 at iteration 5\n" |
| 573 | + "Objective value is 39293.304688 at iteration 1\n", |
| 574 | + "Objective value is 24506.697266 at iteration 2\n", |
| 575 | + "Objective value is 14386.468750 at iteration 3\n", |
| 576 | + "Objective value is 7907.354492 at iteration 4\n", |
| 577 | + "Objective value is 3959.922607 at iteration 5\n" |
579 | 578 | ] |
580 | 579 | } |
581 | 580 | ], |
|
590 | 589 | " # Compute the functional value and gradient for the current\n", |
591 | 590 | " # model estimate\n", |
592 | 591 | " phi, direction = fwi_gradient(model0.vp)\n", |
593 | | - " \n", |
| 592 | + "\n", |
594 | 593 | " # Store the history of the functional values\n", |
595 | 594 | " history[i] = phi\n", |
596 | | - " \n", |
| 595 | + "\n", |
597 | 596 | " # Artificial Step length for gradient descent\n", |
598 | 597 | " # In practice this would be replaced by a Linesearch (Wolfe, ...)\n", |
599 | 598 | " # that would guarantee functional decrease Phi(m-alpha g) <= epsilon Phi(m)\n", |
600 | 599 | " # where epsilon is a minimum decrease constant\n", |
601 | 600 | " alpha = .05 / mmax(direction)\n", |
602 | | - " \n", |
| 601 | + "\n", |
603 | 602 | " # Update the model estimate and enforce minimum/maximum values\n", |
604 | 603 | " update_with_box(model0.vp , alpha , direction)\n", |
605 | | - " \n", |
| 604 | + "\n", |
606 | 605 | " # Log the progress made\n", |
607 | 606 | " print('Objective value is %f at iteration %d' % (phi, i+1))" |
608 | 607 | ] |
|
0 commit comments