MA Alexander Grunewald/Code-Fragmente

Dämpfungzonen

in init.c:UserDefBoundary()

{
    if(side == 0){ // is called *before* the other sides
        double Ymin = g_domBeg[JDIR],Ymax = g_domEnd[JDIR],
            dzs = g_inputParam[Dampingzones_size_cu],
            tau = g_inputParam[Damping_time_cu];
        if (dzs > 0 && g_intStage == 1) {
            int kl,jl,il;
            RBox dom_box;
            DOM_LOOP(kl,jl,il) {
                double damp = 0,y = grid->x[JDIR][jl];
                if (y < Ymin+dzs)
                    damp = Ymin+dzs-y;
                else if (y > Ymax-dzs)
                    damp = y-(Ymax-dzs);
                damp *= 1/(dzs*tau);
                d->Vc[VX2][kl][jl][il] *= MAX(1-damp*g_dt,0.0);
            }
            RBoxDefine (IBEG,IEND,JBEG,JEND,KBEG,KEND,CENTER,&dom_box);
            PrimToCons3D(d->Vc, d->Uc, &dom_box);
        }

        ApplyRestrictions(d, grid);
    }
    else if(side == X2_END) {

⋮

in user_defined_parameters.h

⋮

#define Dampingzones_size_cu            5
#define Damping_time_cu                 6

//
// Specify number of user defined parameters:
//
#define USER_DEF_PARAMETERS             7

Für die beiden neuen Parameter sind natürlich in pluto.ini noch Werte einzusetzen:

• Sinnvoll sind \(0.02\) für die Dampingzones_size_cu, das ist die Breite einer jeden Zone, es bleibt also eine \(0.06\) breite, ungedämpfte Zone in der Mitte. Zum Innenrand einer Zone wird die Dämpfung linear schwächer.
• Für die charakteristische Zeit Damping_time_cu gilt je kleiner, desto stärker die Dämpfung; \(5\cdot 10^{-7}\) ist schon recht stark.

Störung in der Beschleunigung

in init.c

Die Zeile

data->Vc[VX2][kl][jl][il] += perturbation_amplitude * sin(perturbation_mode * 2.0*M_PI * x1 / (Xmax-Xmin)) * ampy;

muss natürlich raus. (Und alles was dazugehört, auch raus.)

in body_force.c

    //
    // Costant gravitational field in the x2-direction:
    //
    g[JDIR] = - G_GravityConstant * SolarMass / (SolarRadius*SolarRadius);

    // perturbation
    if (g_time < g_inputParam[IC_Perturbation_Duration_cu]) {
      // code time unit = 1 year (see user_defined_parameters.h)

      double g_grav = G_GravityConstant * SolarMass / (SolarRadius*SolarRadius),
	    perturbation_mode  = g_inputParam[IC_Perturbation_mode], // dimensionless and integer
	    perturbation_amplitude  = g_inputParam[IC_Perturbation_Amplitude_g]*g_grav,
	    Xmin = grid->xl_glob[IDIR][grid->gbeg[IDIR]],
	    Xmax = grid->xr_glob[IDIR][grid->gend[IDIR]],
	    Ymin = grid->xl_glob[JDIR][grid->gbeg[JDIR]],
	    Ymax = grid->xr_glob[JDIR][grid->gend[JDIR]],
	    ampy = perturbation_amplitude * exp(- 10.0 * fabs(0.5 - (x2-Ymin) / (Ymax-Ymin)));

      g[JDIR] += ampy * sin(perturbation_mode * 2.0*M_PI * x1 / (Xmax-Xmin));
    }

in user_defined_parameters.h

//
// Specify names of user defined parameters:
//
#define IC_Temperature_X2BEG_K          0
#define IC_Temperature_X2END_K          1
#define IC_Perturbation_mode            2
#define IC_Perturbation_Amplitude_g     3
#define IC_Perturbation_Duration_cu     4


//
// Specify number of user defined parameters:
//
#define USER_DEF_PARAMETERS             5

analytische Anfangsprofile

init.c

relevante Zeilen für die initialen Druck- und Dichte-Profile:

        if(jl <= 2) data->Vc[PRS][kl][jl-1][il] = P_0;
        data->Vc[PRS][kl][jl][il] = (1.0 - xi2) / (1.0 + xi2) * data->Vc[PRS][kl][jl-1][il];
        data->Vc[RHO][kl][jl][il] = data->Vc[PRS][kl][jl][il] / (R_UniversalGasConstant / MolarMass(data, kl, jl, il) * Tgas);

Lothar: Die erste davon kann entfallen, die dritte bleibt so, und in der zweiten soll der Druckwert direkt zugewiesen werden: data->Vc[PRS][kl][jl][il] = P_0*pow(1-w*(x2-Ymin),q/w); aber es müssen \(q\) und \(w\) zuvor noch deklariert werden (double q,w; weiter oben) und vor der DOM_LOOP aus den Parametern berechnet werden.