Code en dimension 3

General Discussion
1

Bonjour
j’ai le code suivant qui est en dimension deux d’espace. Quels changements faire pour qu’il soit en dimension 3 d’espace, c’est à dire (x,y,z) \in (0,1)^3 et t \in (0,1)? S’il vous plaît


load "UMFPACK64"
defaulttoUMFPACK64();

load "iovtk"

verbosity=1;

func real qdelta (real t, real s, real w) {
  real t1 = t - rint(t);
  real s1 = s - rint(s);
  real delta=1.0/sqrt(w*pi);
  if ( square(s1)+square(t1) <= square(delta)) 
    return w;
  else  
    return 0.0;
  };


func real q (real t, real s, real w) {
  real t1 = t - rint(t);
  real s1 = s - rint(s);
  real delta=1.0/sqrt(w*pi);
  if ( square(s1)+square(t1) <= square(delta)) 
    return 1;
  else  
    return 0.0;
  };



func real f (real x)
  {
  return 0.25*sin(2*pi*x);
  };

// initial values for the spatial domain
int N=200;
real eps=0.5, w=3.0;

real t=0.00;
real dt=0.1;

int ni=3;       // number of values of h=eps^2 w
real[int] h(ni);     // value of h at each step
real[int] error(ni); // error in norm H1 t each step
real[int] alpha(ni); // alpha =1/log(2). log(error(h))/log(error(h/2))
ofstream fout("error.txt");
fout << "# step h error alpha" << endl;



for (int i=0;i<=ni-1;++i) {

 ofstream fouti("max_"+i+".txt");
fouti << "# t max" << endl;

  mesh Th=square(N,N,[x,y]);
  fespace Vh(Th,P1,periodic=[[2,y],[4,y],[1,x],[3,x]]);
  fespace Wh(Th,P2,periodic=[[2,y],[4,y],[1,x],[3,x]]);
  Vh uh=0.0;
  Vh wh=0.0;
  Vh vh;
  Vh oldU=0.0;
  Vh oldW=0.0;
  real T=1.0;// arbitraire
int[int] Order= [1];

string DataName;

  Vh fh=f;

  Vh nith = 0.0; // for nonlinear iterations
  Vh g0,g1,g2,g3;
  Wh uplot;

  // Nonlinear part is calculated in previous iteration (nith) and becomes a part of the right hand side
  problem het(uh,vh,init=t) =
    int2d(Th)((dx(uh)*dx(vh)+dy(uh)*dy(vh))*dt)
    +int2d(Th)(uh*vh)
    -int2d(Th)(oldU*vh)
    +int2d(Th)(qdelta(x/eps,y/eps,w)*g0*vh*dt) +int2d(Th)(qdelta(x/eps,y/eps,w)*g1*uh*vh*dt)
    -int2d(Th)(fh*vh*dt)
    ;

  problem hom(wh,vh,init=t) =
    int2d(Th)((dx(wh)*dx(vh)+dy(wh)*dy(vh))*dt)
    +int2d(Th)(wh*vh)
    -int2d(Th)(oldW*vh)
    +int2d(Th)(g2*vh*dt) +int2d(Th)(g3*wh*vh*dt)
    -int2d(Th)(fh*vh*dt)
    ;
    
  uplot=q(x/eps,y/eps,w); 
//plot(uplot,fill=true,value=true,ps="q".eps");
savevtk("q_"+i+".vtu",Th,uplot,dataname="qdelta", order=Order);
 cout << "mean value of qdelta: " << int2d(Th)(uplot);
 




real umax = -1e100; //initialisation de max

  for (t=0;t<=T;t+=dt) {
    cout << "Time = " << t << endl;
    fh=f(t);
    oldU=uh;
    oldW=wh;
    int k=0;        
    int noit = 10; // maximum no of nonlinear iterations
    real errnlhet = 0;  // error of nonlinear iteration of het
    real normuh = 0;
    real errnlhom =0; // error of nonlinear iteration of hom
    real normwh=0;
    
    // Nonlinear iterations for heterogeneous solution
    nith = uh; // start nonlinear iterations  with last solution for het
    g0=uh*uh/(1.0+uh)^2;
    g1=1.0/(1.0+uh)^2;
    cout << "    Solving nonlinear heterogeneous problem.\n";
    for(k=1;k<noit; ++k) {
      het;
      normuh = sqrt(int2d(Th)(square(uh))); // L^2 nor of the solution of het
      errnlhet = sqrt(int2d(Th)(square(uh - nith))); // L^2 norm of the error in nonlinear iteration of het
      cout << "        " << k << ": error = " << errnlhet << endl;
      if(errnlhet <= 1E-5*normuh) // stopping criteria. If relative error
      break;        // is less than 1E-5 stop nonlinear iterations of het
      nith = uh; // prepare new iteration
      g0=uh*uh/(1.0+uh)^2;
      g1=1.0/(1.0+uh)^2;
      }
    // Check convergence for het  
    if(k >= noit) {
      cout << "Nonlinear solver for het did not converged.\n";
      cout << "Error = " << errnlhet << endl;
      cout << "|| uh || = " << normuh << endl;
      exit(1); // exit the program with error code 1
      }

    nith= wh; // start nonlinear iteration with last solution for hom
    g2=wh*wh/(1.0+wh)^2;
    g3=1.0/(1.0+wh)^2;
    cout << "    Solving nonlinear homogeneous problem.\n";
    for(k=1;k<noit; ++k) {
      hom;
      normwh = sqrt(int2d(Th)(square(wh))); // L^2 nor of the solution of hom
      errnlhom = sqrt(int2d(Th)(square(wh - nith))); // L^2 norm of the error in nonlinear iteration of hom
      cout << "        " << k << ": error = " << errnlhom << endl;
      if (errnlhom <= 1E-5*normwh) // stopping criteria. If relative error
        break;        // is less than 1E-5 stop nonlinear iterations of het
      nith = wh; // prepare new iteration
      g2=wh*wh/(1.0+wh)^2;
      g3=1.0/(1.0+wh)^2;
      }
    // Check convergence for het  
    if(k >= noit) {
      cout << "Nonlinear solver for hom did not converged.\n";
      cout << "Error = " << errnlhom << endl;
      cout << "|| wh || = " << normwh << endl;
      exit(1); // exit the program with error code 1
      }
  
 }


 



savevtk("u_"+i+".vtu",Th,uh,dataname="uh", order=Order);
savevtk("w_"+i+".vtu",Th,wh,dataname="vh", order=Order);

  real errH1 = sqrt(int2d(Th)(square(uh - wh)) + int2d(Th)(square(dx(uh)-dx(wh)) + square(dy(uh)-dy(wh))));
	
  cout << " ==== errH1 ( eps = " << eps <<", omega = " << w << ", eps^2*omega = " << eps^2*w << ") = "<< errH1 <<endl;
  
  h(i)=eps^2*w;
  error(i)=errH1;
  alpha(i)=i>0?1.0/log(2.0)*log(error(i))/log(error(i-1)):1.0;
  fout << i <<" "<< h(i) <<" "<< error(i) <<" "<< alpha(i) << endl;
  fout.flush;


  eps=eps/2; 
  w=2*w;

  }
'''freefem
1 Like