MUAworkers now within inversion

GHilmarG · Jan 10, 2024 · c263468 · c263468
1 parent 9dbbeab
commit c263468
Show file tree

Hide file tree

Showing 8 changed files with 60 additions and 22 deletions.
diff --git a/HessianAC.m b/HessianAC.m
@@ -1,7 +1,12 @@
-function Hessian=HessianAC(p,lambda,plb,pub,UserVar,CtrlVar,MUA,BCs,F,l,InvStartValues,Priors,Meas,BCsAdjoint,RunInfo)
 
 
-[J,dJdp,Hessian,JGHouts,F,RunInfo]=JGH(p,plb,pub,UserVar,CtrlVar,MUA,BCs,F,l,InvStartValues,Priors,Meas,BCsAdjoint,RunInfo) ; 
+
+
+
+function Hessian=HessianAC(p,lambda,plb,pub,UserVar,CtrlVar,MUA,BCs,F,l,InvStartValues,Priors,Meas,BCsAdjoint,RunInfo,MUAworkers)
+
+
+[J,dJdp,Hessian,JGHouts,F,RunInfo,MUAworkers]=JGH(p,plb,pub,UserVar,CtrlVar,MUA,BCs,F,l,InvStartValues,Priors,Meas,BCsAdjoint,RunInfo,MUAworkers) ; 
 
 
 

diff --git a/InvertForModelParameters.m b/InvertForModelParameters.m
@@ -16,6 +16,8 @@
 %% Define inverse parameters and anonymous function returning objective function, directional derivative, and Hessian
 %
 
+MUAworkers=[]; % replace later with MUA field
+
 
 F=InvStartValues2F(CtrlVar,MUA,F,InvStartValues,Priors,Meas) ;
 [F.b,F.s,F.h,F.GF]=Calc_bs_From_hBS(CtrlVar,MUA,F.h,F.S,F.B,F.rho,F.rhow);
@@ -35,15 +37,18 @@
 
 
 CtrlVar.Inverse.ResetPersistentVariables=1;
-[J0,dJdp,Hessian,JGHouts,F,RunInfo]=JGH(p0,plb,pub,UserVar,CtrlVar,MUA,BCs,F,l,InvStartValues,Priors,Meas,BCsAdjoint,RunInfo);
+
+
+
+[J0,dJdp,Hessian,JGHouts,F,RunInfo,MUAworkers]=JGH(p0,plb,pub,UserVar,CtrlVar,MUA,BCs,F,l,InvStartValues,Priors,Meas,BCsAdjoint,RunInfo,MUAworkers);
 CtrlVar.Inverse.ResetPersistentVariables=0;
 % The parameters passed in the anonymous function are those that exist at the time the anonymous function is created.
 
 
 CtrlVar.WriteRunInfoFile=0;
-func=@(p) JGH(p,plb,pub,UserVar,CtrlVar,MUA,BCs,F,l,InvStartValues,Priors,Meas,BCsAdjoint,RunInfo);
+func=@(p) JGH(p,plb,pub,UserVar,CtrlVar,MUA,BCs,F,l,InvStartValues,Priors,Meas,BCsAdjoint,RunInfo,MUAworkers);
 
-Hfunc=@(p,lambda) HessianAC(p,lambda,plb,pub,UserVar,CtrlVar,MUA,BCs,F,l,InvStartValues,Priors,Meas,BCsAdjoint,RunInfo);
+Hfunc=@(p,lambda) HessianAC(p,lambda,plb,pub,UserVar,CtrlVar,MUA,BCs,F,l,InvStartValues,Priors,Meas,BCsAdjoint,RunInfo,MUAworkers);
 
 fprintf('\n +++++++++++ At start of inversion:  \t J=%-g \t I=%-g \t R=%-g  |grad|=%g \n \n',J0,JGHouts.MisfitOuts.I,JGHouts.RegOuts.R,norm(dJdp))
 
@@ -111,8 +116,10 @@
         error('what case? ')
     end
 
+
     F=p2F(CtrlVar,MUA,p,F,Meas,Priors);
-    [J,dJdp,Hessian,JGHouts,F]=JGH(p,plb,pub,UserVar,CtrlVar,MUA,BCs,F,l,InvStartValues,Priors,Meas,BCsAdjoint,RunInfo);
+
+    [J,dJdp,Hessian,JGHouts,F,RunInfo,MUAworkers]=JGH(p,plb,pub,UserVar,CtrlVar,MUA,BCs,F,l,InvStartValues,Priors,Meas,BCsAdjoint,RunInfo,MUAworkers);
     fprintf('\n +++++++++++ At end of inversion:  \t J=%-g \t I=%-g \t R=%-g  |grad|=%g \n \n',J,JGHouts.MisfitOuts.I,JGHouts.RegOuts.R,norm(dJdp))
 
 

diff --git a/JGH.m b/JGH.m
@@ -1,4 +1,9 @@
-function [J,dJdp,Hessian,JGHouts,F,RunInfo]=JGH(p,plb,pub,UserVar,CtrlVar,MUA,BCs,F,l,InvStartValues,Priors,Meas,BCsAdjoint,RunInfo)
+
+
+
+
+
+function [J,dJdp,Hessian,JGHouts,F,RunInfo,MUAworkers]=JGH(p,plb,pub,UserVar,CtrlVar,MUA,BCs,F,l,InvStartValues,Priors,Meas,BCsAdjoint,RunInfo,MUAworkers)
 
 
 
@@ -7,7 +12,7 @@
 
 persistent ubP vbP
 
-narginchk(14,14)
+narginchk(15,15)
 CtrlVar.nargoutJGH=nargout;
 
 if nargout==1
@@ -47,7 +52,7 @@
     error( ' C nan ') ; 
 end
 
-[UserVar,RunInfo,F,l,dFduv]= uv(UserVar,RunInfo,CtrlVar,MUA,BCs,F,l);
+[UserVar,RunInfo,F,l,dFduv,~,~,MUAworkers]= uv(UserVar,RunInfo,CtrlVar,MUA,BCs,F,l,MUAworkers);
 
 
 if nargout==1

diff --git a/SSTREAM2dNR2.m b/SSTREAM2dNR2.m
@@ -1,11 +1,15 @@
-function  [UserVar,F,l,Kuv,Ruv,RunInfo,L]=SSTREAM2dNR2(UserVar,CtrlVar,MUA,BCs,F,l,RunInfo)
+% function  [UserVar,F,l,Kuv,Ruv,RunInfo,L]=SSTREAM2dNR2(UserVar,CtrlVar,MUA,BCs,F,l,RunInfo)
+
+
+
+function  [UserVar,RunInfo,F,l,Kuv,Ruv,L,MUAworkers]=SSTREAM2dNR2(UserVar,RunInfo,CtrlVar,MUA,BCs,F,l,MUAworkers)
 
 
     % Solves SSA/SSTREAM for u and v
 
 
-    nargoutchk(7,7)
-    narginchk(7,7)
+    nargoutchk(8,8)
+    narginchk(8,8)
 
 
     tStart=tic;
@@ -20,7 +24,7 @@
    % RunInfo.Forward.uvIterations(CtrlVar.CurrentRunStepNumber)=NaN;  
    % RunInfo.Forward.Residual=NaN; BackTrackInfo.iarm=NaN;
 
-    Kuv=[] ; Ruv=[]; MUAworkers=[]; 
+    Kuv=[] ; Ruv=[]; 
 
 
     % MLC=BCs2MLC(CtrlVar,MUA,BCs);
@@ -114,7 +118,7 @@
 
 
 
-    CtrlVar.uvAssembly.ZeroFields=true;   CtrlVar.uvMatrixAssembly.Ronly=true ; MUAworkers=[]; 
+    CtrlVar.uvAssembly.ZeroFields=true;   CtrlVar.uvMatrixAssembly.Ronly=true ; 
     fext0=KRTFgeneralBCs(CtrlVar,MUA,F,MUAworkers); % RHS with velocities set to zero, i.e. only external forces
 
     %% New normalization idea, 10 April 2023
@@ -130,7 +134,7 @@
     %%
 
 
-    CtrlVar.uvAssembly.ZeroFields=false;   CtrlVar.uvMatrixAssembly.Ronly=true ; MUAworkers=[]; 
+    CtrlVar.uvAssembly.ZeroFields=false;   CtrlVar.uvMatrixAssembly.Ronly=true ; 
     Ruv=KRTFgeneralBCs(CtrlVar,MUA,F,MUAworkers);     % RHS with calculated velocities, i.e. difference between external and internal forces
 
     RunInfo.CPU.Solution.uv=0;

diff --git a/Ua2D.m b/Ua2D.m
@@ -277,7 +277,8 @@
     F.m=InvFinalValues.m ; 
     F.n=InvFinalValues.n ;
 
-    [UserVar,RunInfo,F,l,drdu,Ruv,Lubvb]= uv(UserVar,RunInfo,CtrlVar,MUA,BCs,F,l);
+    MUAworkers=[] ; % replace once this becomes` a part of MUA
+    [UserVar,RunInfo,F,l,drdu,Ruv,Lubvb,MUAworkers]= uv(UserVar,RunInfo,CtrlVar,MUA,BCs,F,l,MUAworkers);
 
 
     if CtrlVar.Inverse.WriteRestartFile

diff --git a/uv.m b/uv.m
@@ -1,9 +1,13 @@
-function [UserVar,RunInfo,F,l,Kuv,Ruv,Lubvb]= uv(UserVar,RunInfo,CtrlVar,MUA,BCs,F,l)
 
 
 
-nargoutchk(4,7);
-narginchk(7,7)
+
+
+function [UserVar,RunInfo,F,l,Kuv,Ruv,Lubvb,MUAworkers]= uv(UserVar,RunInfo,CtrlVar,MUA,BCs,F,l,MUAworkers)
+
+
+nargoutchk(4,8);
+narginchk(8,8)
 
 tdiagnostic=tic;
 
@@ -80,7 +84,8 @@
 
         if CtrlVar.InfoLevel >= 10 ; fprintf(CtrlVar.fidlog,' Starting SSTREAM diagnostic step. \n') ;  end
 
-        [UserVar,F,l,Kuv,Ruv,RunInfo,Lubvb]=SSTREAM2dNR2(UserVar,CtrlVar,MUA,BCs,F,l,RunInfo);
+        [UserVar,RunInfo,F,l,Kuv,Ruv,Lubvb,MUAworkers]=SSTREAM2dNR2(UserVar,RunInfo,CtrlVar,MUA,BCs,F,l,MUAworkers);
+        %[UserVar,F,l,Kuv,Ruv,RunInfo,Lubvb]=SSTREAM2dNR2(UserVar,CtrlVar,MUA,BCs,F,l,RunInfo);
 
         if ~RunInfo.Forward.uvConverged
             fprintf('uv forward calculation did not converge. Resetting ub and vb and solving again.\n')

diff --git a/uvMatrixAssembly.m b/uvMatrixAssembly.m
@@ -45,7 +45,7 @@
 
 
                 if CtrlVar.Parallel.uvAssembly.spmd.isOn
-                    MUAworkers=[]; 
+                    % MUAworkers=[]; 
                     tSPMD=tic ;  [RuvSPMD,KuvSPMD,Tint,Fext,MUAworkers]=uvMatrixAssemblySSTREAM_SPMD(CtrlVar,MUA,F,MUAworkers); tSPMD=toc(tSPMD) ;
                     fprintf(' tSeq=%f sec     tSPMD=%f sec \t MUA.Nnodes=%i \t     norm(Ruv-RuvSPMD)/norm(Ruv)=%g \t     norm(diag(Kuv)-diag(KuvSPMD))/norm(diag(Kuv))=%g \n',...
                         tSeq,tSPMD,MUA.Nnodes,norm(full(Ruv-RuvSPMD))/norm(full(Ruv)),norm(diag(Kuv)-diag(KuvSPMD))/norm(diag(Kuv)))

diff --git a/uvMatrixAssemblySSTREAM_SPMD.m b/uvMatrixAssemblySSTREAM_SPMD.m
@@ -6,6 +6,12 @@
 
 narginchk(4,4)
 
+persistent iCount
+
+if isempty(iCount)
+    iCount=0;
+end
+
 
 if isempty(CtrlVar.Parallel.uvAssembly.spmd.nWorkers)
     poolobj = gcp;
@@ -90,13 +96,18 @@
 tSum=toc(tSum) ;
 
 
+spmd ; rr=[] ; kk=[] ; rrsum=[] ; kksum=[] ; end
 
+% for ii=1:nW ; rr{ii}=[] ; kk{ii}=[] ; rrsum{ii}=[] ; kksum{ii}= 0 ; end
 
 Tint=[] ; Fext=[];
 
 
 if CtrlVar.Parallel.isTest
-    fprintf("uvMatrixAssemblySSTREAM_SPMD: Creating partition arrays %f sec. \t Building arrays on workers %f sec. \t SPMD Assembly %f sec. \t Summing up results from workers %f sec.  \n",tPartition,tBuild,tAssembly,tSum)
+
+    iCount=iCount+1;
+    fprintf("uvMatrixAssemblySSTREAM_SPMD (%i): Creating partition arrays %f sec. \t Building arrays on workers %f sec. \t SPMD Assembly %f sec. \t Summing up results from workers %f sec.  \n",...
+        iCount,tPartition,tBuild,tAssembly,tSum)
 end