/* $Header: /Users/ikriest/CVS/mops/tmm_misfit.c,v 1.1 2015/11/17 14:18:51 ikriest Exp $ */ /* $Name: mops-2_0 $*/ #define MIN(x, y) (((x) < (y)) ? (x) : (y)) #define MAX(x, y) (((x) > (y)) ? (x) : (y)) #include #include #include #include "petscmat.h" #include "petsc_matvec_utils.h" #include "tmm_main.h" #include "tmm_forcing_utils.h" #include "tmm_profile_utils.h" #include "tmm_profile_data.h" #include "tmm_timer.h" #include "uvok.h" #include "uvok_biogeochem_misfit_data.h" #include "UVOK_TMM_OPTIONS.h" #define TR v[0] #undef __FUNCT__ #define __FUNCT__ "iniMisfit" PetscErrorCode iniMisfit(PetscScalar tc, PetscInt Iter, PetscInt numTracers, Vec *v) { PetscErrorCode ierr; PetscBool flg; PetscScalar zero = 0.0; PetscViewer fd; /* Add your code here */ averageCost = PETSC_FALSE; /* IK: added for misfit function - start */ /* Type of cost function; so far, only annual averages to PO4, O2 and NO3 available */ ierr = PetscOptionsGetString(PETSC_NULL,"-misfit_file",misfitFile,PETSC_MAX_PATH_LEN-1,&flg);CHKERRQ(ierr); if (!flg) { strcpy(misfitFile,""); sprintf(misfitFile,"%s","misfit.txt"); } ierr = PetscPrintf(PETSC_COMM_WORLD,"misfit will be written to %s\n",misfitFile);CHKERRQ(ierr); ierr = PetscFOpen(PETSC_COMM_WORLD,misfitFile,"w",&misfitf);CHKERRQ(ierr); ierr = PetscOptionsHasName(PETSC_NULL,"-average_cost",&averageCost);CHKERRQ(ierr); if (averageCost) { ierr = iniStepTimer("cost_", Iter0, &costTimer);CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD,"Cost for misfit function will be computed starting at (and including) time step: %d\n", costTimer.startTimeStep);CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD,"Cost for misfit function will be computed over %d time steps\n", costTimer.numTimeSteps);CHKERRQ(ierr); /* Read fractional volume of each box for scaling according to box size */ ierr = VecDuplicate(TR,&globVolFrac);CHKERRQ(ierr); ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"volume_fraction.petsc",FILE_MODE_READ,&fd);CHKERRQ(ierr); ierr = VecLoad(globVolFrac,fd);CHKERRQ(ierr); /* IntoVector */ ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD,"Read volume file for weighting\n");CHKERRQ(ierr); /* wyao 20171204 read in volume fruction for iron data*/ #ifdef O_TMM_opt_iron ierr = VecDuplicate(TR,&globVolFrac_iron);CHKERRQ(ierr); ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"volume_fraction_iron.petsc",FILE_MODE_READ,&fd);CHKERRQ(ierr); ierr = VecLoad(globVolFrac_iron,fd);CHKERRQ(ierr); /* IntoVector */ ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD,"Read volume file for weighting of iron\n");CHKERRQ(ierr); #endif /* Read vector of weights for individual tracer types and gridboxes */ ierr = VecDuplicate(TR,&wbgc1);CHKERRQ(ierr); ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"weights.petsc",FILE_MODE_READ,&fd);CHKERRQ(ierr); ierr = VecLoad(wbgc1,fd);CHKERRQ(ierr); /* IntoVector */ ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr); ierr = VecDuplicate(TR,&wbgc2);CHKERRQ(ierr); ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"weights.petsc",FILE_MODE_READ,&fd);CHKERRQ(ierr); ierr = VecLoad(wbgc2,fd);CHKERRQ(ierr); /* IntoVector */ ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr); ierr = VecDuplicate(TR,&wbgc3);CHKERRQ(ierr); ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"weights.petsc",FILE_MODE_READ,&fd);CHKERRQ(ierr); ierr = VecLoad(wbgc3,fd);CHKERRQ(ierr); /* IntoVector */ ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr); /*wyao 20171204*/ #ifdef O_TMM_opt_iron ierr = VecDuplicate(TR,&wbgc4);CHKERRQ(ierr); ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"weights_iron.petsc",FILE_MODE_READ,&fd);CHKERRQ(ierr); ierr = VecLoad(wbgc4,fd);CHKERRQ(ierr); /* IntoVector */ ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr); #endif ierr = PetscPrintf(PETSC_COMM_WORLD,"Read misfit file(s) for weighting\n");CHKERRQ(ierr); /* construct total weight: for now, multiply weight with fractional volume */ ierr = VecDuplicate(TR,&w1);CHKERRQ(ierr); ierr = VecDuplicate(TR,&w2);CHKERRQ(ierr); ierr = VecDuplicate(TR,&w3);CHKERRQ(ierr); ierr = VecPointwiseMult(w1,wbgc1,globVolFrac); ierr = VecPointwiseMult(w2,wbgc2,globVolFrac); ierr = VecPointwiseMult(w3,wbgc3,globVolFrac); /*wyao 20171204*/ #ifdef O_TMM_opt_iron ierr = VecDuplicate(TR,&w4);CHKERRQ(ierr); ierr = VecPointwiseMult(w4,wbgc4,globVolFrac_iron); #endif ierr = PetscPrintf(PETSC_COMM_WORLD,"Set up weight arrays\n");CHKERRQ(ierr); /* Read the observations */ ierr = VecDuplicate(TR,&obgc1);CHKERRQ(ierr); ierr = VecDuplicate(TR,&obgc2);CHKERRQ(ierr); ierr = VecDuplicate(TR,&obgc3);CHKERRQ(ierr); /*wyao 20171204*/ #ifdef O_TMM_opt_iron ierr = VecDuplicate(TR,&obgc4);CHKERRQ(ierr); #endif ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"po4ini.petsc",FILE_MODE_READ,&fd);CHKERRQ(ierr); ierr = VecLoad(obgc1,fd);CHKERRQ(ierr); /* IntoVector */ ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD,"Read PO4 observations\n");CHKERRQ(ierr); ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"o2ini.petsc",FILE_MODE_READ,&fd);CHKERRQ(ierr); ierr = VecLoad(obgc2,fd);CHKERRQ(ierr); /* IntoVector */ ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD,"Read O2 observations\n");CHKERRQ(ierr); ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"no3ini.petsc",FILE_MODE_READ,&fd);CHKERRQ(ierr); ierr = VecLoad(obgc3,fd);CHKERRQ(ierr); /* IntoVector */ ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD,"Read NO3 observations\n");CHKERRQ(ierr); /*wyao 20171204*/ #ifdef O_TMM_opt_iron ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"Fe_obs.petsc",FILE_MODE_READ,&fd);CHKERRQ(ierr); ierr = VecLoad(obgc4,fd);CHKERRQ(ierr); /* IntoVector */ ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD,"Read dfe observations\n");CHKERRQ(ierr); #endif /* Initialize the cost function vectors */ ierr = VecDuplicate(TR,&mbgc1);CHKERRQ(ierr); ierr = VecSet(mbgc1,zero);CHKERRQ(ierr); ierr = VecGetArray(mbgc1,&localmbgc1);CHKERRQ(ierr); ierr = VecDuplicate(TR,&mbgc1avg);CHKERRQ(ierr); ierr = VecSet(mbgc1avg,zero);CHKERRQ(ierr); ierr = VecDuplicate(TR,&mbgc2);CHKERRQ(ierr); ierr = VecSet(mbgc2,zero);CHKERRQ(ierr); ierr = VecGetArray(mbgc2,&localmbgc2);CHKERRQ(ierr); ierr = VecDuplicate(TR,&mbgc2avg);CHKERRQ(ierr); ierr = VecSet(mbgc2avg,zero);CHKERRQ(ierr); ierr = VecDuplicate(TR,&mbgc3);CHKERRQ(ierr); ierr = VecSet(mbgc3,zero);CHKERRQ(ierr); ierr = VecGetArray(mbgc3,&localmbgc3);CHKERRQ(ierr); ierr = VecDuplicate(TR,&mbgc3avg);CHKERRQ(ierr); ierr = VecSet(mbgc3avg,zero);CHKERRQ(ierr); /*wyao 20171204*/ #ifdef O_TMM_opt_iron ierr = VecDuplicate(TR,&mbgc4);CHKERRQ(ierr); ierr = VecSet(mbgc4,zero);CHKERRQ(ierr); ierr = VecGetArray(mbgc4,&localmbgc4);CHKERRQ(ierr); ierr = VecDuplicate(TR,&mbgc4avg);CHKERRQ(ierr); ierr = VecSet(mbgc4avg,zero);CHKERRQ(ierr); #endif ierr = VecDuplicate(TR,&cbgc1);CHKERRQ(ierr); ierr = VecSet(cbgc1,zero);CHKERRQ(ierr); ierr = VecDuplicate(TR,&cbgc2);CHKERRQ(ierr); ierr = VecSet(cbgc2,zero);CHKERRQ(ierr); ierr = VecDuplicate(TR,&cbgc3);CHKERRQ(ierr); ierr = VecSet(cbgc3,zero);CHKERRQ(ierr); /*wyao 20171204*/ #ifdef O_TMM_opt_iron ierr = VecDuplicate(TR,&cbgc4);CHKERRQ(ierr); ierr = VecSet(cbgc4,zero);CHKERRQ(ierr); #endif /* ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"cbgc1.petsc",FILE_MODE_WRITE,&fdcbgc1avg);CHKERRQ(ierr);*/ /* ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"cbgc2.petsc",FILE_MODE_WRITE,&fdcbgc2avg);CHKERRQ(ierr);*/ /* ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"cbgc3.petsc",FILE_MODE_WRITE,&fdcbgc3avg);CHKERRQ(ierr);*/ /*wyao 20171204*/ #ifdef O_TMM_opt_iron /* ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"cbgc4.petsc",FILE_MODE_WRITE,&fdcbgc4avg);CHKERRQ(ierr); */ #endif // costCount=0; } else { ierr = PetscPrintf(PETSC_COMM_WORLD,"Only average cost function available, currently.\n");CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD,"Specify -average_cost, and start time step and number of time steps to average.\n");CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD,"No misfit will be computed. \n");CHKERRQ(ierr); } /* IK: added for misfit function - end */ return 0; } /* -----------------------------------------------------------------------------------------------------------*/ /* -----------------------------------------------------------------------------------------------------------*/ /* -----------------------------------------------------------------------------------------------------------*/ /* -----------------------------------------------------------------------------------------------------------*/ #undef __FUNCT__ #define __FUNCT__ "calcMisfit" PetscErrorCode calcMisfit(PetscScalar tc, PetscInt iLoop, PetscInt numTracers, Vec *v) { /* Note: tc and iLoop are the time and step at the end of the current time step. */ PetscErrorCode ierr; PetscScalar zero = 0.0, one = 1.0, minusone = -1.0 ; if (Iter0+iLoop>=costTimer.startTimeStep) { /* note: startTimeStep is ABSOLUTE time step */ /* Add your code here */ if (averageCost) { /* only option, currently */ /* IK : Add all tracer snapshots, and increase counter by one */ if (costTimer.count<=costTimer.numTimeSteps) { ierr = VecAXPY(mbgc1avg,one,mbgc1);CHKERRQ(ierr); ierr = VecAXPY(mbgc2avg,one,mbgc2);CHKERRQ(ierr); ierr = VecAXPY(mbgc3avg,one,mbgc3);CHKERRQ(ierr); /*wyao 20171204*/ #ifdef O_TMM_opt_iron ierr = VecAXPY(mbgc4avg,one,mbgc4);CHKERRQ(ierr); #endif costTimer.count++; // costCount+1; } } } return 0; } /* -----------------------------------------------------------------------------------------------------------*/ /* -----------------------------------------------------------------------------------------------------------*/ /* -----------------------------------------------------------------------------------------------------------*/ /* -----------------------------------------------------------------------------------------------------------*/ #undef __FUNCT__ #define __FUNCT__ "writeMisfit" PetscErrorCode writeMisfit(PetscScalar tc, PetscInt iLoop, PetscInt numTracers, Vec *v) { /* Note: tc and iLoop are the time and step at the end of the current time step. */ PetscErrorCode ierr; PetscScalar zero = 0.0, one = 1.0, minusone = -1.0 ; if (Iter0+iLoop>=costTimer.startTimeStep) { /* note: startTimeStep is ABSOLUTE time step */ /* Add your code here */ if (averageCost){ /* only option, currently */ /* IK : added for misfit function : start */ if (costTimer.count==costTimer.numTimeSteps) { /* time to compute cost function and write to file */ ierr = PetscPrintf(PETSC_COMM_WORLD,"Computing cost function of time average over %d steps at time %10.5f, step %d\n", costTimer.count, tc, Iter0+iLoop);CHKERRQ(ierr); /* IK : average model over a year; this will be stored again in mbgc1avg, ... */ ierr = VecScale(mbgc1avg,1.0/costTimer.count);CHKERRQ(ierr); ierr = VecScale(mbgc2avg,1.0/costTimer.count);CHKERRQ(ierr); ierr = VecScale(mbgc3avg,1.0/costTimer.count);CHKERRQ(ierr); /*wyao 20171204*/ #ifdef O_TMM_opt_iron ierr = VecScale(mbgc4avg,1.0/costTimer.count);CHKERRQ(ierr); #endif ierr = PetscPrintf(PETSC_COMM_WORLD,"Computed annual average concentrations \n");CHKERRQ(ierr); /* IK : (later) write vectors of cost functions for a posteriori analysis */ /* ierr = VecView(mbgc1avg,fdcbgc1avg);CHKERRQ(ierr);*/ /* ierr = VecView(mbgc2avg,fdcbgc2avg);CHKERRQ(ierr);*/ /* ierr = VecView(mbgc3avg,fdcbgc3avg);CHKERRQ(ierr);*/ /*wyao 20171204*/ #ifdef O_TMM_opt_iron /* ierr = VecView(mbgc4avg,fdcbgc4avg);CHKERRQ(ierr);*/ #endif /* ierr = PetscPrintf(PETSC_COMM_WORLD,"Wrote annual average concentrations \n");CHKERRQ(ierr); */ /* IK : get global average model tracers, e.g., if I ever want to account for the bias Gavembgc1, ... */ ierr = VecDot(mbgc1avg,globVolFrac,&Gavembgc1);CHKERRQ(ierr); ierr = VecDot(mbgc2avg,globVolFrac,&Gavembgc2);CHKERRQ(ierr); ierr = VecDot(mbgc3avg,globVolFrac,&Gavembgc3);CHKERRQ(ierr); /*wyao 20171204*/ #ifdef O_TMM_opt_iron ierr = VecDot(mbgc4avg,globVolFrac_iron,&Gavembgc4);CHKERRQ(ierr); #endif ierr = PetscPrintf(PETSC_COMM_WORLD,"Computed global mean model tracers \n");CHKERRQ(ierr); /* copy global mean tracer to array cbgc1, ... and perform subsequent operations on this */ ierr = VecCopy(mbgc1avg,cbgc1);CHKERRQ(ierr); ierr = VecCopy(mbgc2avg,cbgc2);CHKERRQ(ierr); ierr = VecCopy(mbgc3avg,cbgc3);CHKERRQ(ierr); /*wyao 20171204*/ #ifdef O_TMM_opt_iron ierr = VecCopy(mbgc4avg,cbgc4);CHKERRQ(ierr); #endif /* IK : subtract observations from average model concentrations; this will be stored in cbgc1avg, ...., again */ ierr = VecAXPY(cbgc1,minusone,obgc1);CHKERRQ(ierr); ierr = VecAXPY(cbgc2,minusone,obgc2);CHKERRQ(ierr); ierr = VecAXPY(cbgc3,minusone,obgc3);CHKERRQ(ierr); /*wyao 20171204*/ #ifdef O_TMM_opt_iron ierr = VecAXPY(cbgc4,minusone,obgc4);CHKERRQ(ierr); #endif ierr = PetscPrintf(PETSC_COMM_WORLD,"Subtracted model-observations \n");CHKERRQ(ierr); /* IK : square the difference between observations and average model concentrations */ /* Note: I use the more VecPow instead of VecNorm, because I want to be as generic as possible, e.g., use functions proposed by Evans, 2003; result is again in cbgcavg1, ..., again */ ierr = VecPow(cbgc1,2.0);CHKERRQ(ierr); ierr = VecPow(cbgc2,2.0);CHKERRQ(ierr); ierr = VecPow(cbgc3,2.0);CHKERRQ(ierr); /*wyao 20171204*/ #ifdef O_TMM_opt_iron ierr = VecPow(cbgc4,2.0);CHKERRQ(ierr); #endif ierr = PetscPrintf(PETSC_COMM_WORLD,"Squared the difference \n");CHKERRQ(ierr); /* IK : sum all local misfits, and scale with fractional volume at the same time; result is scalar value Gavecost1, ... */ ierr = VecDot(cbgc1,w1,&Gavecost1);CHKERRQ(ierr); ierr = VecDot(cbgc2,w2,&Gavecost2);CHKERRQ(ierr); ierr = VecDot(cbgc3,w3,&Gavecost3);CHKERRQ(ierr); /*wyao 20171204*/ #ifdef O_TMM_opt_iron ierr = VecDot(cbgc4,w4,&Gavecost4);CHKERRQ(ierr); #endif ierr = PetscPrintf(PETSC_COMM_WORLD,"Computed overall sum, weighted by volume and weight \n");CHKERRQ(ierr); /* IK : get global average observed concentrations Gavebgc1, ... for normalization (needed to sum up different quantities; could use any other scaling, e.g., global variance */ ierr = VecDot(obgc1,globVolFrac,&Gaveobgc1);CHKERRQ(ierr); ierr = VecDot(obgc2,globVolFrac,&Gaveobgc2);CHKERRQ(ierr); ierr = VecDot(obgc3,globVolFrac,&Gaveobgc3);CHKERRQ(ierr); /*wyao 20171204*/ #ifdef O_TMM_opt_iron ierr = VecDot(obgc4,globVolFrac_iron,&Gaveobgc4);CHKERRQ(ierr); #endif ierr = PetscPrintf(PETSC_COMM_WORLD,"Computed global mean observations \n");CHKERRQ(ierr); /* IK : Get the total cost */ RGavecost1=sqrt(Gavecost1); RGavecost2=sqrt(Gavecost2); RGavecost3=sqrt(Gavecost3); /*wyao 20171204*/ #ifdef O_TMM_opt_iron RGavecost4=sqrt(Gavecost4); Gcost = RGavecost1/Gaveobgc1+RGavecost2/Gaveobgc2+RGavecost3/Gaveobgc3+RGavecost4/Gaveobgc4; #else Gcost = RGavecost1/Gaveobgc1+RGavecost2/Gaveobgc2+RGavecost3/Gaveobgc3; #endif ierr = PetscPrintf(PETSC_COMM_WORLD,"Computed overall cost function \n");CHKERRQ(ierr); /* write misfit to file - perhaps better binary? ierr = writeBinaryScalarData("misfit.bin",&Gcost,1,PETSC_TRUE); */ /* some diagnostic output for checking */ ierr = PetscPrintf(PETSC_COMM_WORLD,"Check translation\n");CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD,"Writing cost function at time %10.5f, step %d:: %10.5f\n", tc, Iter0+iLoop, Gcost);CHKERRQ(ierr); /*wyao 20171204*/ #ifdef O_TMM_opt_iron ierr = PetscPrintf(PETSC_COMM_WORLD,"Components are PO4 %10.8f O2 %10.8f NO3 %10.8f DFE %10.8f\n", RGavecost1,RGavecost2,RGavecost3,RGavecost4);CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD,"Average observations are PO4 %10.8f O2 %10.8f NO3 %10.8f DFE %10.8f\n", Gaveobgc1,Gaveobgc2,Gaveobgc3,Gaveobgc4);CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD,"Average model tracers are PO4 %10.8f O2 %10.8f NO3 %10.8f DFE %10.8f\n", Gavembgc1,Gavembgc2,Gavembgc3,Gavembgc4);CHKERRQ(ierr); #else ierr = PetscPrintf(PETSC_COMM_WORLD,"Components are PO4 %10.8f O2 %10.8f NO3 %10.8f\n", RGavecost1,RGavecost2,RGavecost3);CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD,"Average observations are PO4 %10.8f O2 %10.8f NO3 %10.8f\n", Gaveobgc1,Gaveobgc2,Gaveobgc3);CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD,"Average model tracers are PO4 %10.8f O2 %10.8f NO3 %10.8f\n", Gavembgc1,Gavembgc2,Gavembgc3);CHKERRQ(ierr); #endif ierr = PetscFPrintf(PETSC_COMM_WORLD,misfitf,"%10.8f\n",Gcost);CHKERRQ(ierr); ierr = updateStepTimer("cost_", Iter0+iLoop, &costTimer);CHKERRQ(ierr); // costCount = 0; } /* costCount==costNumTimeSteps */ } else { /* no other option that average cost, currently */ ierr = PetscPrintf(PETSC_COMM_WORLD,"No cost function available \n");CHKERRQ(ierr); } } /* Iter0+iLoop>=costTimer.startTimeStep */ /* IK : added for misfit function : end */ return 0; } /* -----------------------------------------------------------------------------------------------------------*/ /* -----------------------------------------------------------------------------------------------------------*/ /* -----------------------------------------------------------------------------------------------------------*/ /* -----------------------------------------------------------------------------------------------------------*/ #undef __FUNCT__ #define __FUNCT__ "finalizeMisfit" PetscErrorCode finalizeMisfit(PetscScalar tc, PetscInt Iter, PetscInt numTracers) { PetscErrorCode ierr; /* Add your code here */ /* IK : added for misfit function */ ierr = VecDestroy(&mbgc1);CHKERRQ(ierr); ierr = VecDestroy(&mbgc1avg);CHKERRQ(ierr); ierr = VecDestroy(&obgc1);CHKERRQ(ierr); ierr = VecDestroy(&cbgc1);CHKERRQ(ierr); ierr = PetscViewerDestroy(&fdcbgc1avg);CHKERRQ(ierr); ierr = VecDestroy(&mbgc2);CHKERRQ(ierr); ierr = VecDestroy(&mbgc2avg);CHKERRQ(ierr); ierr = VecDestroy(&obgc2);CHKERRQ(ierr); ierr = VecDestroy(&cbgc2);CHKERRQ(ierr); ierr = PetscViewerDestroy(&fdcbgc2avg);CHKERRQ(ierr); ierr = VecDestroy(&mbgc3);CHKERRQ(ierr); ierr = VecDestroy(&mbgc3avg);CHKERRQ(ierr); ierr = VecDestroy(&obgc3);CHKERRQ(ierr); ierr = VecDestroy(&cbgc3);CHKERRQ(ierr); ierr = PetscViewerDestroy(&fdcbgc3avg);CHKERRQ(ierr); /*wyao 20171204*/ #ifdef O_TMM_opt_iron ierr = VecDestroy(&mbgc4);CHKERRQ(ierr); ierr = VecDestroy(&mbgc4avg);CHKERRQ(ierr); ierr = VecDestroy(&obgc4);CHKERRQ(ierr); ierr = VecDestroy(&cbgc4);CHKERRQ(ierr); ierr = PetscViewerDestroy(&fdcbgc4avg);CHKERRQ(ierr); #endif ierr = PetscFClose(PETSC_COMM_WORLD,misfitf);CHKERRQ(ierr); /* IK : added for misfit function */ return 0; }