This piece of code does the parameter search that gives the maximum error estimator Click to go back
Download here : ParameterSearch.idp
// // Take a representative set for domain parameter mu // Want to find choice of mu that gives largest error estimator real errorestimator0=0; real bestk0; real ds=1;// particion fina // cout<<"Start Search... "<< n <<endl; for(int j=1;j<=(k0max-k0min);j++) { // // TO MAKE SURE DO NOT REPITE PARAMETER AND INDEX // I AM SURE IT CAN BE MADE BETTER! // //cout<<" Parameter index j = "<< j <<" Best parameter index jbest "<< jbest[cont]<<endl; //cout<< " Check (j!=jbest[cont]) "<<(j!=jbest[cont])<<endl; // //IDEA // WHEN sum=n THE INDEX j HAS NOT BEEN USED BEFORE // real sum=0; //for(int s=0;s<=cont;s++) for(int s=0;s<n;s++) {//cout<<" j= "<<j <<" jbest[s]= "<<jbest[s]<< " bool "<<(j!=jbest[s])<<endl; sum=sum+(j!=jbest[s]); //sum=cont ---> j!=jbest[s] cierto (igual a 1) para todos // sum<cont---> j!=jbest[s] falso (igual a 0) para alguno //cout<<" Valor de sum "<<sum <<endl; } if (sum==n) // NOW WE USE INDEX j { k0=k0min+ds*j; // PARAMETER ACCORDING TO j // // FIX ALL THE PARAMETERS ACCORDING TO THIS ONE // include "parameters.idp"; // //SOLVE THE PROBLEM ACCORDING TO THIS PARAMETER WITH REDUCED BASIS // include "ReduceBasisSolver.idp" // // SEE SOLUTION // // isovalues // for (int i=0;i<viso.n;i++) viso[i]=RBu[].min+i*(RBu[].max-RBu[].min)/mm; plot(cmm= " Search "+n+ ": RB Solution for k0 = "+k0,RBu,viso=viso(0:viso.n-1),hsv=colorhsv,value=1); // // COMPUTE ERROR ESTIMATE: // RIESZ REPRESENTATION OF L(*)-a(RBu,*) en Vh. // // WE HAVE ALREADY Riesz Representation of L(*) // // AND Riesz Representation of a(RBu,*) is given by // // (Riesz_a(RBu,*), v)_V=a(RBu,v), for all v en Vh // // COMPUTATION OF a(RBu,*) // varf aa0(unused,w)=int3d(Th3,0,2,4,6,8,10)( dx(RBu)*dx(w)+dy(RBu)*dy(w)+dz(RBu)*dz(w) ); real[int] aa0vec=aa0(0,Vh); varf aa1(unused,w)=int3d(Th3,1,5,9)( dx(RBu)*dx(w)+dy(RBu)*dy(w)+dz(RBu)*dz(w) ); real[int] aa1vec=aa1(0,Vh); varf aa2(unused,w)=int3d(Th3,3,7,11)( dx(RBu)*dx(w)+dy(RBu)*dy(w)+dz(RBu)*dz(w) ); real[int] aa2vec=aa2(0,Vh); // varf aa3(unused,w)=int2d(Th3,10,12,14,16,18,110,11,15,19,13,17,111,30,32 ,34,36,38,310,31,35,39,33,37,311,40,44,48,41,45,49 ,20,24,28,23,27,211,199,201,299,301,399,401,499,501 ,599,601,699,700,701)(RBu*w); real[int] aa3vec=aa3(0,Vh); // varf aa4(unused,w)=int3d(Th3)(RBu*w); real[int] aa4vec=aa4(0,Vh); // // VECTOR ASSOCIATED TO a(RBu,*) ACTING ON FEM BASIS // Vh rhsforRieszaRBu; rhsforRieszaRBu[]=k0*aa0vec; rhsforRieszaRBu[]=rhsforRieszaRBu[]+kleft*aa1vec; rhsforRieszaRBu[]=rhsforRieszaRBu[]+kright*aa2vec; rhsforRieszaRBu[]=rhsforRieszaRBu[]+aa*aa4vec; rhsforRieszaRBu[]=rhsforRieszaRBu[]+Bi*aa3vec; // // RIESZ REPRESENTATION COMES FROM INVERTING THE MATRIX // OF THE SCALAR PRODUCT // Vh RieszaRBu; RieszaRBu[]=Ascalarproduct^-1*rhsforRieszaRBu[]; // // VECTOR THAT HAS RIESZ REPRESENTATION OF RESIDUAL // Vh diff=RieszaRBu-Rieszf; // // ERROR ESTIMATOR IS THE NORM OF DIFF // real part1= int3d(Th3)( dx(diff)^2+dy(diff)^2+dz(diff)^2); real part2=int3d(Th3)(diff^2); real normdiff=sqrt(part1+part2*invvol); real errorestimator= normdiff/k0min; cout <<"Search "<<n<<" and iter j = "<<j<<". For k0 = "<<k0 <<": error estimator "<<errorestimator<<endl; // // SEARCH FOR THE LARGEST ERROR ESTIMATOR // if(errorestimator>errorestimator0) { errorestimator0=errorestimator; bestk0=k0; jbest[n]=j; } } } cout <<"Search "<<n<<". Best k0 = "<<bestk0 <<": error estimator "<<errorestimator0<<endl;