Commit 89c8abf6 authored by Florian Wittkamp's avatar Florian Wittkamp

Fix conjugate gradient and low pass filtering

This fix will reset the conjugate gradient history each time a frequency or workflow stage is switched.
If such a stage is switched we will have a "new" objective function, thus, it makes no sense to use the old gradients for speeding up the convergence.
This bug was the reason why the objective function could not be reduced for 2-3 iterations after a new frequency stage was appplied.
parent dd778c01
......@@ -99,6 +99,9 @@ int main(int argc, char **argv){
int SOURCE_SHAPE_OLD=0;
int SOURCE_SHAPE_OLD_SH=0;
/* Variables for conjungate gradient */
int PCG_iter_start=1;
/* Variables for L-BFGS */
int LBFGS_NPAR=3;
int LBFGS_iter_start=1;
......@@ -222,9 +225,6 @@ int main(int argc, char **argv){
at the receiver positions to calculate the adjoint sources and to do
the backpropagation; look at function saveseis_glob.c to see that every
NDT sample for the forward modeled wavefield is written to su files*/
lsnap=iround(TSNAP1/DT); /* first snapshot at this timestep */
lsamp=NDT;
/* output of parameters to log-file or stdout */
if (MYID==0) write_par(FP);
......@@ -504,13 +504,7 @@ int main(int argc, char **argv){
cjm = vector(1,L);
}
/*nf=4;
nfstart=4;*/
NTST=20;
NTSTI=NTST/DTINV;
nxny=NX*NY;
nxnyi=(NX/IDXI)*(NY/IDYI);
/* Parameters for step length calculations */
......@@ -2900,17 +2894,25 @@ int main(int argc, char **argv){
/* ----------------------------------------------------------------------*/
/* ----------- Preconditioned Conjugate Gradient Method (PCG) ----------*/
/* ----------------------------------------------------------------------*/
if((GRAD_METHOD==1)){
if(GRAD_METHOD==1){
if( (iter-PCG_iter_start) < 2 ) {
fprintf(FP,"\n\n ----- Conjugate Gradient Method -----");
fprintf(FP,"\n Will not use second last gradient for conjugate gradient");
if( (iter-PCG_iter_start) < 1 ) {
fprintf(FP,"\n Will not use last gradient for conjugate gradient");
}
}
if (WAVETYPE==1 || WAVETYPE==3) PCG(waveconv, taper_coeff, nsrc, srcpos, recpos, ntr_glob, iter, C_vp, gradp, nfstart_jac, waveconv_u, C_vs, gradp_u, waveconv_rho, C_rho, gradp_rho,Vs_avg,F_LOW_PASS);
if (WAVETYPE==2 || WAVETYPE==3) PCG_SH(taper_coeff, nsrc, srcpos, recpos, ntr_glob, iter, nfstart_jac, waveconv_u_z, C_vs, gradp_u_z, waveconv_rho_z, C_rho, gradp_rho_z,Vs_avg,F_LOW_PASS);
if (WAVETYPE==1 || WAVETYPE==3) PCG(waveconv, taper_coeff, nsrc, srcpos, recpos, ntr_glob, iter, C_vp, gradp, nfstart_jac, waveconv_u, C_vs, gradp_u, waveconv_rho, C_rho, gradp_rho,Vs_avg,F_LOW_PASS,PCG_iter_start);
if (WAVETYPE==2 || WAVETYPE==3) PCG_SH(taper_coeff, nsrc, srcpos, recpos, ntr_glob, iter, nfstart_jac, waveconv_u_z, C_vs, gradp_u_z, waveconv_rho_z, C_rho, gradp_rho_z,Vs_avg,F_LOW_PASS,PCG_iter_start);
}
/* ---------------------------------------------------------*/
/* ----------- Beginn Joint Inversion PSV and SH ----------*/
/* ---------------------------------------------------------*/
if((FORWARD_ONLY==0)){
if(FORWARD_ONLY==0){
switch (WAVETYPE) {
case 2:
/* If only SH is inverted, set these gradients to the actual ones */
......@@ -2950,7 +2952,7 @@ int main(int argc, char **argv){
/* -------------------------------------------------------------------------*/
/* ----------- Limited Memory - Broyden-Fletcher-Goldfarb-Shanno ----------*/
/* -------------------------------------------------------------------------*/
if((GRAD_METHOD==2)){
if(GRAD_METHOD==2){
/*---------------------*/
/* TAPER */
......@@ -3903,6 +3905,9 @@ int main(int argc, char **argv){
/* Restart L-BFGS at next iteration */
LBFGS_iter_start=iter+1;
/* Restart conjugate gradient at next iteration */
PCG_iter_start=iter+1;
wolfe_SLS_failed=0;
}
......@@ -3933,6 +3938,10 @@ int main(int argc, char **argv){
/* Restart L-BFGS at next iteration */
LBFGS_iter_start=iter+1;
/* Restart conjugate gradient at next iteration */
PCG_iter_start=iter+1;
wolfe_SLS_failed=0;
alpha_SL_old=1;
}
......@@ -3964,6 +3973,10 @@ int main(int argc, char **argv){
/* Restart L-BFGS at next iteration */
LBFGS_iter_start=iter+1;
/* Restart conjungate gradient at next iteration */
PCG_iter_start=iter+1;
wolfe_SLS_failed=0;
alpha_SL_old=1;
}
......
......@@ -24,7 +24,7 @@
#include "fd.h"
void PCG(float ** waveconv, float ** taper_coeff, int nsrc, float ** srcpos, int ** recpos, int ntr_glob, int iter, float C_vp, float ** gradp, int nfstart_jac, float ** waveconv_u, float C_vs, float ** gradp_u, float ** waveconv_rho, float C_rho, float ** gradp_rho, float Vs_avg, float F_LOW_PASS){
void PCG(float ** waveconv, float ** taper_coeff, int nsrc, float ** srcpos, int ** recpos, int ntr_glob, int iter, float C_vp, float ** gradp, int nfstart_jac, float ** waveconv_u, float C_vs, float ** gradp_u, float ** waveconv_rho, float C_rho, float ** gradp_rho, float Vs_avg, float F_LOW_PASS, int PCG_iter_start){
extern int NX, NY, IDX, IDY, SPATFILTER, GRAD_FILTER;
extern int FORWARD_ONLY, SWS_TAPER_GRAD_VERT, SWS_TAPER_GRAD_HOR, SWS_TAPER_GRAD_SOURCES, SWS_TAPER_FILE;
......@@ -41,7 +41,13 @@ void PCG(float ** waveconv, float ** taper_coeff, int nsrc, float ** srcpos, int
int use_conjugate_2=1;
FILE *FP3, *FP4, *FP6, *FP5 = NULL;
/* Check if conjugate gradient can be used */
if( (iter-PCG_iter_start) < 2 ) {
use_conjugate_2=0;
if( (iter-PCG_iter_start) < 1 ) {
use_conjugate_1=0;
}
}
/* ===================================================================================================================================================== */
/* ===================================================== GRADIENT ZP ================================================================================== */
......
......@@ -24,7 +24,7 @@
#include "fd.h"
void PCG_SH(float ** taper_coeff, int nsrc, float ** srcpos, int ** recpos, int ntr_glob, int iter, int nfstart_jac, float ** waveconv_u, float C_vs, float ** gradp_u, float ** waveconv_rho, float C_rho, float ** gradp_rho, float Vs_avg, float F_LOW_PASS){
void PCG_SH(float ** taper_coeff, int nsrc, float ** srcpos, int ** recpos, int ntr_glob, int iter, int nfstart_jac, float ** waveconv_u, float C_vs, float ** gradp_u, float ** waveconv_rho, float C_rho, float ** gradp_rho, float Vs_avg, float F_LOW_PASS, int PCG_iter_start){
extern int NX, NY, IDX, IDY, SPATFILTER, GRAD_FILTER;
extern int FORWARD_ONLY, SWS_TAPER_GRAD_VERT, SWS_TAPER_GRAD_HOR, SWS_TAPER_GRAD_SOURCES, SWS_TAPER_FILE;
......@@ -40,6 +40,14 @@ void PCG_SH(float ** taper_coeff, int nsrc, float ** srcpos, int ** recpos, int
int use_conjugate_2=1;
/* Check if conjugate gradient can be used */
if( (iter-PCG_iter_start) < 2 ) {
use_conjugate_2=0;
if( (iter-PCG_iter_start) < 1 ) {
use_conjugate_1=0;
}
}
/* ===================================================================================================================================================== */
/* ===================================================== GRADIENT Zs ================================================================================== */
/* ===================================================================================================================================================== */
......
......@@ -186,9 +186,9 @@ void taper(float *section, int ns, float fc);
void output_source_signal(FILE *fp, float **signals, int ns, int seis_form);
void PCG(float ** waveconv, float ** taper_coeff, int nsrc, float ** srcpos, int ** recpos, int ntr_glob, int iter, float C_vp, float ** gradp, int nfstart_jac, float ** waveconv_u, float C_vs, float ** gradp_u, float ** waveconv_rho, float C_rho, float ** gradp_rho, float Vs_avg, float F_LOW_PASS);
void PCG(float ** waveconv, float ** taper_coeff, int nsrc, float ** srcpos, int ** recpos, int ntr_glob, int iter, float C_vp, float ** gradp, int nfstart_jac, float ** waveconv_u, float C_vs, float ** gradp_u, float ** waveconv_rho, float C_rho, float ** gradp_rho, float Vs_avg, float F_LOW_PASS, int PCG_iter_start);
void PCG_SH(float ** taper_coeff, int nsrc, float ** srcpos, int ** recpos, int ntr_glob, int iter, int nfstart_jac, float ** waveconv_u, float C_vs, float ** gradp_u, float ** waveconv_rho, float C_rho, float ** gradp_rho, float Vs_avg, float F_LOW_PASS);
void PCG_SH(float ** taper_coeff, int nsrc, float ** srcpos, int ** recpos, int ntr_glob, int iter, int nfstart_jac, float ** waveconv_u, float C_vs, float ** gradp_u, float ** waveconv_rho, float C_rho, float ** gradp_rho, float Vs_avg, float F_LOW_PASS, int PCG_iter_start);
void PML_pro(float * d_x, float * K_x, float * alpha_prime_x, float * a_x, float * b_x,
float * d_x_half, float * K_x_half, float * alpha_prime_x_half, float * a_x_half, float * b_x_half,
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment