time_step.c Source File

00001 /*------------------------------------------------------------------*/
00002 /* SCCS Information: %W%    %G% */
00003 /*------------------------------------------------------------------*/
00004 #include "utilf.h"
00005 #include <readSIunits.h>
00006 #include "markers.h"
00007 #include "surfaces.h"
00008 #include "extra.h"
00009 #include "stresses.h"
00010 #include "global.h"
00011 #include "printxplot.h"
00012 #include "printfree.h"
00013 #include "momentum.h"
00014 #include "tension.h"
00015 #include "advect.h"
00016 #include "make_bc.h"
00017 
00018 extern real minau(real2D au, real2D ap, int nx, int ny, int *imin, int *jmin);
00019 extern real minav(real2D av, real2D ap, int nx, int ny, int *imin, int *jmin);
00020 
00021 /* #define RECONNECT_AXIS */
00022 
00023 void timestep(int t)
00024 {
00025   int imax = 0, jmax = 0, iumax = 0, jumax = 0, i1, i, j;
00026   real umax, vmax, h, hi2;
00027   static real pe, pg;
00028 #ifdef SONO
00029   static real R, Rp;
00030   real a;
00031 #endif
00032   FILE *msgptr;
00033 /*cld*/
00034   real uaxis = 0.0;
00035   int iaxis;
00036   int test = 0;
00037 /*cld*/
00038 
00039   h = 1.0 / (num.nx - 2);
00040   hi2 = 1.0 / (h * h);
00041 
00042   /* ---------------------------------------------------------------
00043      Control CFL condition
00044      ---------------------------------------------------------------- */
00045   if (t % 1 == 0) {
00046     umax = fmodmax(u, num.nx, num.ny, &iumax, &jumax);
00047     vmax = fmodmax(v, num.nx, num.ny, &i, &j);
00048     if (vmax > umax) { umax = vmax; iumax = i; jumax = j; }
00049     if (umax > 0.5) {
00050       printf("Warning> t=%d velmax=%g x=%d y=%d\n", t, 
00051              umax, iumax, jumax);
00052       printf("CFL number too large. Try smaller tau\n");
00053       printxpl(u, v, ap, p, num.nx, num.ny, 
00054                num.tau, num.du, num.dp, "cfl.datb");
00055       msgptr = fopen("cfl.xmgr", "wt");
00056       for (i = 0; i < Nint; i++)
00057         interface_write_markers(interfaces[i], msgptr);
00058       fclose(msgptr);
00059       exit(1);
00060     }
00061   }
00062 
00063   /* ---------------------------------------------------------------
00064      Iteration
00065      ---------------------------------------------------------------- */
00066   /*------------------- Advection ------------------- */
00067   for (i = 0; i < Nint; i++) {
00068     advect(interfaces[i], u, v, ap, num.nx, num.ny, &uaxis); 
00069     /*    repulse(interfaces[i], 0.5, 1e-5, num.lredis); */
00070 #ifdef RECONNECT_AXIS
00071     if (reconnect_axis(&interfaces[i], num.nx, num.ny, 0.1)) {
00072       printf("t = %d: reconnect interface %d\n", t, i);
00073 
00074       msgptr = fopen("reconnect.xmgr", "wt");
00075       interface_write_markers(interfaces[i], msgptr);
00076       fclose(msgptr);
00077 #if 0
00078       interface_splines(interfaces[i], num.nx, num.ny);
00079       printf("Filling fractions ... "); 
00080       interface_fill(interfaces[i], ap, num.nx, num.ny);
00081       for (i1 = 2; i1 <= num.nx; i1++)
00082         for (j = 2; j <= num.ny; j++) {
00083           au[i1][j] = (ap[i1][j] + ap[i1-1][j])/2.;
00084           av[i1][j] = (ap[i1][j] + ap[i1][j-1])/2.;
00085         }
00086       printf("Ok\n"); fflush(stdout);
00087 #endif
00088     }
00089 #endif
00090     interface_splines(interfaces[i], num.nx, num.ny);
00091 #ifdef AXISYMMETRIC_SMOOTHING
00092     interface_axisymmetry(&interfaces[i], num.lredis, num.nx, num.ny);
00093 #else
00094     interface_redis(&interfaces[i], num.lredis);
00095 #endif
00096     interface_splines(interfaces[i], num.nx, num.ny);
00097   }
00098 
00099   pg = phys.po*exp(phys.gamma*log(phys.vo/
00100                    fabs(interfaces_axivolume(interfaces, Nint))));
00101   pe = phys.pext - phys.pe*sin(phys.omega*phys.t);
00102 
00103   a = exp(1/3.*log(fabs(interfaces_axivolume(interfaces, Nint))*3./(4.*PI)));
00104   if (t == 0) Rp = 0.0;
00105   else Rp = a - R;
00106   R = a;
00107 
00108   /* compute interfacial pressure */
00109   interfacial_pressure(interfaces[0], u, v, ap, phys.sigma, phys.visliq, pg,
00110                        num.tau, num.nx, num.ny);
00111 
00112   cut_interface(interfaces[0], &cutp, 0.0, 0.0);
00113   cut_interface(interfaces[0], &cutu, -0.5, 0.0);
00114   cut_interface(interfaces[0], &cutv, 0.0, -0.5);
00115 
00116   interface_surfaces(interfaces[0], cutp, sxp, syp, ap, cc,
00117                      S12, work, num.nx, num.ny);
00118   interface_surfaces(interfaces[0], cutu, sxu, syu, au, a1,
00119                      S12, work, num.nx, num.ny);
00120   interface_surfaces(interfaces[0], cutv, sxv, syv, av, a2,
00121                      S12, work, num.nx, num.ny);
00122 
00123   sx_fill(sxp, ap, 0.0, 0.0, interfaces[0], num.nx, num.ny);
00124   sy_fill(syp, ap, 0.0, 0.0, interfaces[0], num.nx, num.ny);
00125   sx_fill(sxu, au, -0.5, 0.0, interfaces[0], num.nx, num.ny);  
00126   sy_fill(syu, au, -0.5, 0.0, interfaces[0], num.nx, num.ny);
00127   for (j = 2; j <= num.ny; j++)
00128     syu[num.nx+1][j] = (real)j - 1.5;
00129   sx_fill(sxv, av, 0.0, -0.5, interfaces[0], num.nx, num.ny);
00130   for (i = 2; i <= num.nx; i++)
00131     sxv[i][num.ny+1] = (real)num.ny - 1.5;
00132   sy_fill(syv, av, 0.0, -0.5, interfaces[0], num.nx, num.ny);
00133 
00134   bc_vector_div(u, v, ap, interfaces[0], num.nx, num.ny);
00135 
00136   /*--------------------- Momentum equation -------------------------*/
00137   viscous_tensor(u, v, S11, S22, S12,
00138                  phys.visliq,
00139                  num.tau, num.nx, num.ny);
00140 
00141   momentum(u, v,
00142            a1, a2, ap, av,
00143            sxu, syu, sxv, syv,
00144            interfaces[0], cutu, cutv,
00145            rz1dr, rz2dr, rr1dz, rr2dz,
00146            work,
00147            S11, S22, S12,
00148            divs,
00149            phys.visliq,
00150            phys.gx, phys.gy,
00151            num.tau,
00152            num.nx, num.ny);
00153 
00154   interface_tensionu(interfaces[0], cutu, u, v, a1, ap, syu, p, S11, S22, 
00155                      num.nx, num.ny);
00156   interface_tensionv(interfaces[0], cutv, u, v, a2, ap, sxv, av, p, S11, S22,
00157                      num.nx, num.ny);
00158 
00159   bc_vector_bound(u, v, num.nx, num.ny);
00160 
00161   bc_pressure(p, ap, pe*sq(num.tau/h), 
00162               R,
00163               num.nx, num.ny);
00164 
00165   interface_rzdr(interfaces[0], cutp, -0.5, rz1dr, num.nx, num.ny);
00166   interface_rzdr(interfaces[0], cutp, 0.5, rz2dr, num.nx, num.ny);
00167   interface_rrdz(interfaces[0], cutp, -0.5, rr1dz, num.nx, num.ny);
00168   interface_rrdz(interfaces[0], cutp, 0.5, rr2dz, num.nx, num.ny);
00169 
00170   pressure(u, v, p, a1, a2, ap, av, sxp, syp, syu, sxv, rz1dr, rz2dr, 
00171            rr1dz, rr2dz,
00172            divs, work1,
00173            S11, S22, S12, a3, work, 
00174            interfaces[0], num.lredis,
00175            &num.cycles, num.npre, num.npost, num.maxdiv,
00176            num.nx, num.ny, num.ng);
00177   
00178   bc_vector_div(u, v, ap, interfaces[0], num.nx, num.ny);
00179 
00180   /* ---------------------------------------------------------------
00181      Outputs
00182      ---------------------------------------------------------------- */
00183   if (t % out.txplot == 0) {
00184     char s[256];
00185     /*
00186       sprintf(s, "free.%d.fig", t / out.txplot);
00187       msgptr = fopen(s, "wt");
00188       printfree(msgptr, u, v, ap, interfaces[0], num.nx, num.ny);
00189       fclose(msgptr);
00190     */
00191     sprintf(s, "data.%d.datb", t/out.txplot);
00192     printxpl(u, v, ap, p, num.nx, num.ny, num.tau, num.du, num.dp, s);
00193   }
00194   if (t % out.tprint == 0)
00195     {
00196       if (t == 0)
00197         msgptr = fopen("msg", "wt");
00198       else
00199         msgptr = fopen("msg", "at");
00200 
00201       printf("t=%d umax=%g at (%d,%d) ncycle=%d", t, umax, 
00202              iumax, jumax, num.cycles);
00203       if (phys.sigma > 0.0)
00204         printf(" cwn=%g\n", num.tau/sqrt(h*h*h/phys.sigma));
00205       else
00206         printf("\n");
00207       divergence(u, v, ap, sxp, syp, rz1dr, rz2dr, rr1dz, rr2dz, work, 
00208                  num.nx, num.ny);
00209       printf("t=%d div max before=%g div max after=%g\n",
00210              t, 
00211              rdivmax(divs, num.nx, num.ny, &imax, &jmax),
00212              rdivmax(work, num.nx, num.ny, &imax, &jmax));
00213       fprintf(msgptr, "%d %g %g %g %g %g %g %g %g %g %g",
00214               t,
00215               phys.t,
00216               umax * h / num.tau * num.du,
00217               interfaces_axivolume(interfaces, Nint),
00218               findmax(ap, num.nx, num.ny, &imax, &jmax),
00219               rdivmax(work, num.nx, num.ny, &imax, &jmax),
00220               num.tau,
00221               interfaces_x_amp(interfaces, Nint)/(num.nx - 2),
00222               interfaces_y_amp(interfaces, Nint)/(num.ny - 2),
00223               /*              interfaces_ymoment(interfaces, Nint), */
00224               pg*num.dp, 
00225               pe*num.dp);
00226       if (interfaces[0].n - 1 % 2 == 0) {
00227         i = (interfaces[0].n - 1)/2;
00228         fprintf(msgptr, " %g %g\n", 
00229                 sqrt(sq(interfaces[0].x[i] - 0.5*(num.nx - 2) - 2.0) + 
00230                      sq(interfaces[0].y[i] - 2.0)) -
00231                 exp(1./3.*log(3./(4.*PI)*fabs(interfaces_axivolume(interfaces, Nint)))),
00232                 exp(1./3.*log(3./(4.*PI)*fabs(interfaces_axivolume(interfaces, Nint)))));
00233       }
00234       else {
00235         i = (interfaces[0].n - 1)/2;
00236         fprintf(msgptr, " %g %g\n", 
00237                 0.5*(sqrt(sq(interfaces[0].x[i] - 0.5*(num.nx - 2) - 2.0) + 
00238                           sq(interfaces[0].y[i] - 2.0)) +
00239                      sqrt(sq(interfaces[0].x[i+1] - 0.5*(num.nx - 2) - 2.0) + 
00240                           sq(interfaces[0].y[i+1] - 2.0))) -
00241                 exp(1./3.*log(3./(4.*PI)*fabs(interfaces_axivolume(interfaces, Nint)))),
00242                 exp(1./3.*log(3./(4.*PI)*fabs(interfaces_axivolume(interfaces, Nint)))));
00243       }
00244       printf("mass=%g\n", sumfield(cc, num.nx, num.ny));
00245       printf("min c: %g  max c: %g\n", 
00246              findmin(ap, num.nx, num.ny, 
00247                      &imax, &jmax),
00248              findmax(ap, num.nx, num.ny,
00249                      &imax, &jmax)); 
00250       a = minau(au, ap, num.nx, num.ny, &imax, &jmax);
00251       printf("min au: %f at %d,%d\n", a, imax, jmax);
00252       a = minav(av, ap, num.nx, num.ny, &imax, &jmax);
00253       printf("min av: %f at %d,%d\n", a, imax, jmax);
00254       printf("dp (SI): %g\n",
00255              dp(p, cc, num.nx, num.ny) * h * h / sq(num.tau) * num.dp);
00256       fflush(stdout);
00257       fclose(msgptr);
00258     }
00259 
00260   if (t % out.tmark == 0)
00261     {
00262       if (t == 0)
00263         msgptr = fopen("markers.xmgr", "wt");
00264       else
00265         msgptr = fopen("markers.xmgr", "at");
00266       printf("volume: %g  Nint: %d\n", interfaces_axivolume(interfaces, Nint), 
00267              Nint);
00268       for (i = 0; i < Nint; i++)
00269         interface_write_markers(interfaces[i], msgptr);
00270 //      interface_write_markers_speed(phys.t, interfaces[i], u, v, 
00271 //                      num.tau, num.du, num.nx, num.ny, msgptr);
00272       fclose(msgptr);
00273 #if 0
00274       if (t == 0)
00275         msgptr = fopen("splines.xmgr", "wt");
00276       else
00277         msgptr = fopen("splines.xmgr", "at");
00278       for (i = 0; i < Nint; i++)
00279         interface_write_splines(interfaces[i], msgptr, 1000);
00280       fclose(msgptr);
00281 #endif
00282       if (t == 0)
00283         msgptr = fopen("curvature.xmgr", "wt");
00284       else
00285         msgptr = fopen("curvature.xmgr", "at");
00286       for (i = 0; i < Nint; i++)
00287         curvature_write(interfaces[i], msgptr);
00288       fclose(msgptr);
00289 
00290       if (t == 0)
00291         msgptr = fopen("ipressure.xmgr", "wt");
00292       else
00293         msgptr = fopen("ipressure.xmgr", "at");
00294       for (i = 0; i < interfaces[0].n; i++) 
00295         fprintf(msgptr, "%g %g\n", 
00296                 interfaces[0].t[i]/interfaces[0].t[interfaces[0].n-1], 
00297                 interfaces[0].p[i]);
00298       fprintf(msgptr, "\n");
00299       fclose(msgptr);
00300 
00301       if (t == 0)
00302         msgptr = fopen("radius.xmgr", "wt");
00303       else
00304         msgptr = fopen("radius.xmgr", "at");
00305       for (i = 0; i < Nint; i++) {
00306         fprintf(msgptr, "\n");
00307         for (j = 0; j < interfaces[i].n - 1; j++)
00308           fprintf(msgptr, "%g %g\n", interfaces[i].t[j]/
00309                   interfaces[i].t[interfaces[i].n - 1], 
00310                   sqrt(sq(interfaces[i].x[j] - 0.5*(num.nx - 2) - 2.0) + 
00311                        sq(interfaces[i].y[j] - 2.0)));
00312       }
00313       fclose(msgptr);
00314 /*cld*/
00315       if (t == 0)
00316           msgptr = fopen("axis.xmgr", "wt");
00317       else
00318           msgptr = fopen("axis.xmgr", "at");
00319 
00320       iaxis = interfaces[0].x[0];
00321 //      printf("%g %d \n", phys.t*num.tausi, iaxis);
00322 
00323       fprintf(msgptr, "%g %g %g %g %g\n", phys.t, 
00324               (interfaces[0].x[0]-2.0)*h, 
00325               u[iaxis-1][1] * h / num.tau * num.du, 
00326               u[iaxis-1][2] * h / num.tau * num.du, 
00327               u[iaxis-1][3] * h / num.tau * num.du);
00328 
00329       fclose(msgptr);
00330 
00331       if (t == 0)
00332               msgptr = fopen("velocity.xmgr", "wt");
00333       else
00334               msgptr = fopen("velocity.xmgr", "at");
00335 
00336 //      if ((interfaces[0].x[0]-2.0)*h > 0.616406 && test == 0){
00337 /* Double Size Bubble */
00338 //      if ((interfaces[0].x[0]-2.0)*h > 0.918281 && test == 0){
00339 /* Double Size Bubble */
00340               iaxis = interfaces[0].x[0];
00341               test = 1;
00342               fprintf(msgptr, "%g %g %g %g %g\n", phys.t, 
00343                               (interfaces[0].x[0]-2.0)*h,
00344                               u[iaxis-1][1] * h / num.tau * num.du,
00345                               u[iaxis-1][2] * h / num.tau * num.du,
00346                               u[iaxis-1][3] * h / num.tau * num.du);
00347 //            exit(1);
00348 //      }
00349       fclose(msgptr);
00350               
00351 
00352 /*cld*/
00353     }
00354 
00355   phys.t += num.tau;
00356 
00357 /*cld
00358         for (j = 1; j < interfaces[0].n ; j++){
00359                 if (interfaces[0].y[j] <= 3.0 && 
00360                     interfaces[0].y[j] < interfaces[0].y[j-1] &&
00361                     interfaces[0].x[j] < interfaces[0].x[j-1])
00362                 exit(1);
00363         }
00364 cld*/
00365 
00366 #ifndef NO_ADAPTATIVE
00367   adaptative(u, v, p, &num.tau, num.cycles, num.nx, num.ny);
00368 #endif
00369 }
00370 
00371 
00372 
00373 
00374 
00375