Simulation of a real capacitor to calculate its capacitance

Question

The goal is to simulate a real flat plate capacitor in 2D. The potential of the armature is fixed to +V and -V, on the edges of the simulation cell the potential is zero. Between the plates of the capacitor a dielectric material is placed. The capacity is also calculated by numerically determining the charge on the plates and integrating the flow of the electric field. The result is compared with a formula that does not take into account edge effects.

I am looking for ways to make it more efficient, and how to use some functions, because I don't think it looks really clear.

#include <iostream>
#include <string>
#include <cmath>
#include <fstream>

using namespace std;

int main(int argc, const char * argv[])
{
    double Lx,Ly;//lunghezze lati cella di simulazione
    int Nx,Ny;//numero punti lungo x e y
    double Px,Py;//posizione rispetto all'origine del vertice della prima piastra
    double Zx,Zy;//posizione rispetto all'angolo del bordo a destra del vertice della seconda piastra
    double Lc;//spessore della piastra
    double Hc;//altezza piastra
    int    dist_type;//tipo di distribuzione di carica: 1-uniforme
    double thrs;//valore soglia per convergenza
    double V;
    double Nstampa;

    int   iCount;
    int   phi_bordo=0.;//funzione phi sul bordo
    double hx,hy;
    double diff=1e10;

    ofstream fphi,fcharge, ffield;


    cout.precision(10);
    cout << "Lato cella lungo x\n";
    cin >>  Lx;
    cout << "Lato cella lungo y\n";
    cin >>  Ly;
    cout << "Numero punti lungo x\n";
    cin >>  Nx;
    cout << "Numero punti lungo y\n";
    cin >>  Ny;
    cout << "Vertice piastra 1 in x\n";
    cin >>  Px;
    cout << "Vertice piastra 2 in y\n";
    cin >>  Py;
    cout << "Vertice piastra 1 in x\n";
    cin >>  Zx;
    Zx = Lx - Zx;
    cout << "Vertice piastra 2 in y\n";
    cin >>  Zy;
    cout << "Spessore piastra \n";
    cin >>  Lc;
    cout << "Altezza piastra \n";
    cin >>  Hc;
    cout << "Potenziale piastre (in Volt) \n";
    cin >>  V;
    cout << "Soglia per convergenza \n";
    cin >>  thrs;

    double* phi0= new double[Nx*Ny];//funzione da trovare
    double* phi1= new double[Nx*Ny];//funzione da trovare
    double* phi2= new double[Nx*Ny];//funzione da trovare
    double* ftilde=new double[Nx*Ny];//distribuzione di carica e fattori bordo
    double* mphi= new double[Nx*Ny];

    hx=Lx/(Nx-1.);
    hy=Ly/(Ny-1.);


    iCount=0; //iCount corrisponde all'indice l delle dispense l = j*Nx + i con i,j indici rispettivamente di xi e yj
              //infatti stiamo calcolando ftilde che è l'array indicizzato con l
    for(int j=0;j<Ny;j++){
        double y;
        y=Ly*j/(Ny-1.);
        for(int i=0;i<Nx;i++){
            double x;
            x=Lx*i/(Nx-1.);
            if(x>=Px && x<=(Px+Lc) && y>=Py && y<=(Py+Hc)) {
                ftilde[iCount]=V;
            }
            else if(x>=Zx && x<=(Zx+Lc) && y>=Zy && y<=(Zy+Hc)){
              ftilde[iCount]=-V;
            }
            else{
                ftilde[iCount]=0.;
            }
            //potenziale sul bordo della cella
            if(i==0 || i==(Nx-1) || j==0 || j==(Ny-1)){
                ftilde[iCount]+=phi_bordo;
            }
            iCount++;

        }
    }

    //metti phi0 in condizioni iniziali quindi pari a 0
    for(int i=0;i<Nx*Ny;i++){
        phi1[i]=ftilde[i]; //vettore che costruiamo man mano e che è quello che ci serve: phi = ftilde + Mftilde + M^2ftilde ...
        phi0[i]=ftilde[i]; //vettore che moltiplicato per Mtilde iterativamente ci dà M^n*ftilde
        phi2[i]=0.; //vettore che poniamo ogni volta uguale a phi1 del ciclo precedente in modo da poter calcolare la differenza |phi1-phi2|
    }


    while(diff>thrs){
         for(int i=0;i<Nx*Ny;i++){
             mphi[i]=0; //mphi corrisponde al vettore dato dal prodotto M*phi
         }
        iCount=0.; //indice l di ftilde e di mphi
        //aggiungi Mtilde/phi0
        for(int j=0;j<Ny;j++){
            for(int i=0;i<Nx;i++){
                int jCount; //indice m che scorre lungo phi
                double div;
            //termine diagonale
          //   mphi[iCount]=phi0[iCount]*(-2.*(1./pow(hx,2)+1./pow(hy,2)));
           //termini +1, -1

//punti ai bordi della cella sono esclusi dalla prime 4 condizioni, mentre i punti interni alle due piastre
//che sono SEMPRE DI BORDO e su cui non voglio calcolare l'equazione di Poisson li escludo con la condizione ftilde[icount]=0,
//che comprende tutti i punti non contenuti nelle piastre

                 if(i!=0 && i!=(Nx-1) && j!=0 && j !=(Ny-1) && ftilde[iCount] == 0){
                     div=(-2.*(1./pow(hx,2)+1./pow(hy,2))); //corrisponde a Mll, ovvero denominatore di tutti i termini di matrice Mtilde
// Ci ricordiamo che l = j*Nx + i;
                     if(i>0){
                         jCount=j*Nx+(i-1); // m = l-1
                         mphi[iCount]+=-phi0[jCount]/pow(hx,2)/div;
                     }
                     if(i<Nx-1){
                         jCount=j*Nx+(i+1);// m = l+1
                         mphi[iCount]+=-phi0[jCount]/pow(hx,2)/div;
                     }
                     if(j>0){
                         jCount=(j-1)*Nx+i;// m = l - Nx
                         mphi[iCount]+=-phi0[jCount]/pow(hy,2)/div;
                     }
                     if(j<Ny-1){
                         jCount=(j+1)*Nx+i;// m = l + Nx
                          mphi[iCount]+=-phi0[jCount]/pow(hy,2)/div;
                     }
                }
                 iCount++;

            }
        }

//il motivo per cui non calcolo in ogni caso mphi per l indice che corrisponde a indici di punti di bordo è perchè
//la matrice Mtilde sarebbe pari a zero per tutti i termini diagonali e per tutte le righe in cui
//l è indice di un punto di bordo: infatti i valori che phi assume ai bordi sono già definiti in ftilde a cui aggiungo
//Mftilde + M^2ftilde.. che hanno invece componente di indice l di bordo nulla

     //metti in phi1
        for(int i=0;i<Nx*Ny;i++){
            phi1[i]+=mphi[i];
            phi0[i]=mphi[i];
        }

        //calcola differenza
        diff=0;
        for(int i=0;i<Nx*Ny;i++){
            diff+=pow(phi1[i]-phi2[i],2);
        }
        diff=sqrt(diff);
       cout << "Differenza: " << diff << "\n";
        //copia phi0 in phi1
        for(int i=0;i<Nx*Ny;i++){
            phi2[i]=phi1[i];
        }

    }


    double * Emat = new double [Nx*Ny*2];
    double ** E = new double* [2];

    for (int i=0;i<2*Nx*Ny;i++){
      E[i] = &Emat[Nx*Ny*i];
    }

//trovo componente x del campo elettrico
for(int j=0;j<Ny;j++){
    for(int i=2;i<Nx-2;i++){
          E[0][Nx*j+i] = (-phi1[Nx*j+(i+2)]+8*phi1[Nx*j+(i+1)]-8*phi1[Nx*j+(i-1)]+phi1[Nx*j+(i-2)])/(12*hx);
      }
  }
//trovo componente y del campo elettrico
for(int j=2;j<Ny-2;j++){
    for(int i=0;i<Nx;i++){
      E[1][Nx*j+i] = (-phi1[Nx*(j+2)+i]+8*phi1[Nx*(j+1)+i]-8*phi1[Nx*(j-1)+i]+phi1[Nx*(j-2)+i])/(12*hy);
    }
  }

//calcoliamo il flusso del campo elettrico uscente da un rettangolo chiuso che circonda la prima piastra
//mi basta calcolare flusso attorno a prima piastra perchè quello della seconda sarà in modulo uguale per simmetria
//integrale metodo Simpson
int j1, i1, j2, i2; //indici di coordinate vertici rettangolo su cui calcolare flusso
j1 = floor((Py/2.)*(Ny-1.)/Ly);
j2 = floor((Py+Hc+0.5*(Ly-Py-Hc))*(Ny-1.)/Ly);
i1 = floor((Px/2.)*(Nx-1.)/Lx);
i2 = floor((Px+Lc+0.5*(Zx-Px-Lc))*(Nx-1.)/Lx);

double Esbx = -E[1][Nx*j1+i1-1]; //prodotto scalare E*u con E uno step più a sx di i1 e u uscente da base sotto rettangolo
double Edbx = -E[1][Nx*j1+i2+1]; //u uscente base sotto rettangolo

double Esby = -E[0][Nx*(j1-1)+i1]; //prodotto scalare E*u con E uno step più in basso di j1 e u uscente da lato sinistro rettangolo
double Esay = -E[0][Nx*(j2+1)+i1];

double Edby = E[0][Nx*(j1-1)+i2]; //u uscente da lato destro
double Eday = E[0][Nx*(j2+1)+i2];

double Esax = E[1][Nx*j2+i1-1]; //u uscente da base sopra rettangolo
double Edax = E[1][Nx*j2+i2+1];

//flusso attraverso base bassa
double Intg = 0;
 Intg = (1./3.)*hx*(Esbx+Edbx);
for(int i=i1; i<i2+1; i++){
          Intg+=((i%2+1)*2)/3.*(-E[1][Nx*j1+i])*hx;
}
double Flusso1 = Intg;
cout << "Flusso1 " << Flusso1 << " Vm" << endl;
Flusso1 = abs(Flusso1);

Intg =0;
//flusso attraverso lato sinistro
 Intg = (1./3.)*hy*(Esby+Esay);
for(int j=j1; j<j2+1; j++){
          Intg+=((j%2+1)*2)/3.*(-E[0][Nx*j+i1])*hy;
}
double Flusso2 = Intg;
cout << "Flusso2 " << Flusso2 << " Vm" << endl;
Flusso2 = abs(Flusso2);

Intg =0;
//flusso attraverso base alta
 Intg= (1./3.)*hx*(Esax+Edax);
for(int i=i1; i<i2+1; i++){
          Intg+=((i%2+1)*2)/3.*(E[1][Nx*j2+i])*hx;
}
double Flusso3 = Intg;
cout << "Flusso3 " << Flusso3 << " Vm" << endl;
Flusso3 = abs(Flusso3);

Intg =0;
//flusso attraverso lato destro
 Intg = (1./3.)*hy*(Edby+Eday);
for(int j=j1; j<j2+1; j++){
          Intg+=((j%2+1)*2)/3.*(E[0][Nx*j+i2])*hy;
}
double Flusso4 = Intg;
cout << "Flusso4 " << Flusso4 << " Vm" << endl;
Flusso4 = abs(Flusso4);

//flusso corrisponde a carica interna/cost dielettrica nel vuoto

double Q = (Flusso1+Flusso2+Flusso3+Flusso4)*(8.8541878e-12);
cout << "Carica contenuta in una piastra: " << Q << " Coulomb" << endl;
double Cap = 0;
Cap = (double) Q/(2.*V);
cout << "Capacità reale " <<  Cap << " Farad" << endl;
double Capid = 0;
Capid = 8.8541878e-12*Hc/(Zx-Px-Lc);
cout << "Capacità ideale " <<  Capid << " Farad" << endl;
double Rap = 0;
Rap = Cap/Capid;
cout << "Rapporto capacità reale e ideale R " <<  Rap << endl;



    //scrive grafici
    fcharge.open("carica.dat",ios::out);
    fphi.open("phi.dat",ios::out);
    ffield.open("campo.dat",ios::out);
    ffield.precision(5);
    fcharge.precision(5);
    fphi.precision(5);
    iCount=0;
    for(int j=0;j<Ny;j++){
        double y;
        y=Ly*j/(Ny-1.);
        for(int i=0;i<Nx;i++){
            double x;
            x=Lx*i/(Nx-1.);
            fcharge << x << "  " << y << "  " << ftilde[iCount]*(-2.*(1./pow(hx,2)+1./pow(hy,2))) << '\n';
            fphi << x << "  " << y << "  " << -phi1[iCount] << '\n';
            iCount++;
        }
        fphi << '\n';
        fcharge << '\n';
    }

    for(int j=0;j<Ny;j++){
        double y;
        y=Ly*j/(Ny-1.);
        for(int i=floor((Px)*(Nx-1.)/Lx)-12;i<floor((Zx+Lc)*(Nx-1.)/Lx)+12;i++){
            double x;
            x=Lx*i/(Nx-1.);
            ffield << x << " " << y << " " << E[0][j*Nx+i] << " " << E[1][j*Nx+i] << '\n';
        }
    }

    fcharge.close();
    fphi.close();
    ffield.close();

    delete [] phi0;
    delete [] phi1;
    delete [] phi2;
    delete [] ftilde;
    delete [] mphi;
    //delete [] Emat;
    //delete [] E;


}

Answer 1

About mixing languages

I see your program is intended for Italian users. But what language are the developers speaking? The comments are also all in Italian, but the names of the variables are a mix between Italian (for example, Nstampa, Flusso1) and English (dist_type, iCount).

Be consistent in which language you are using. Of course, it doesn't help that the keywords of the C++ language, as well as the names of functions and types in the standard library, are all in English. So I would recommend that you write all your code in English. If you expect only Italian developers to ever work on this project, then using Italian in the comments is fine. And the text written to std::cout should of course be in the language your users expect, so if that is Italian, that is fine too.

Avoid short variable names

Short variable names save some typing, but it's hard to tell what they mean exactly. If I see Lx used in the code, I'm not sure the L stands for length, but even if I know that, what length? The length of the capacitor itself? Or one of the cells? A better variable name would be cell_size_x.

Why are there three arrays phi0, phi1 and phi2? What does phi mean anyway? Flux? Phase? Angle? The comments where they are declared don't mention that at all.

If you are using some formulas from a book or a paper, and want to use the same symbols in your code as they used in those formulas, that's a valid reason to keep the variable names short, but in that case it is essential that you add a link to the book or paper in your code.

Group related information together in structs

I see you create 5 arrays holding information about each cell in the simulation. Consider creating a struct instead that holds all the information about a single cell, then you only have to create a single array:

struct Cell {
    double phi0;
    double phi1;
    double phi2;
    double ftilde;
    double mphi;
};

Cell* cells = new Cell[Nx * Ny];
…
for (int i = 0; i < Nx * Ny; i++) {
    cells[i].phi1 = cells[i].ftilde;
    cells[i].phi0 = cells[i].ftilde;
    cells[i].phi2 = 0;
}

You also have several 2D vectors. It is helpful to create a struct for that as well. For example:

struct Vec2 {
    double x;
    double y;
};

Vec2 cell_size; // replaces Lx and Ly

Those become even more useful if you add operator overloads, so you can do math on them like you do with regular ints and doubles. Creating this yourself can be a lot of work though, so I recommend you choose a library that implements vector types for you, for example GLM, but there are many others.

Avoid manual memory allocation

Instead of allocating an array with new, and having to remember to delete it, use std::vector instead. You have a memory leak, because you commented out the delete statements for Emat and E.

There are other benefits to using std::vector: you can use a range-for statement in some places. Consider:

struct Cell {…};

std::vector<Cell> cells(Nx * Ny);
…
for (auto& cell: cells) {
    cell.phi1 = cell.ftilde;
    cell.phi0 = cell.ftilde;
    cell.phi2 = 0;
}

Give constants names

I see the number 8.8541878e-12 used several times in your code. Instead of writing that out every time, define a constant for it, for example:

constexpr double e0 = 8.854187812813e-12; // Vacuum permitivity in Farads/meter

Note that e0 resembles the mathematical notation \$\epsilon_0\$ that is commonly used for this constant, and the comment explains it in more detail in case someone not familiar with this is reading the code.

Split the code into multiple functions

You already noticed yourself that the code doesn't look great, and that it could benefit from functions. A good way to start is to rewrite main() into something that describes the things you need to do at a very high level. For example, you have the following steps in your program:

• Read parameters
• Initialize the data
• Run the simulation
• Process the results
• Write the results to disk

You could write your main() like so:

struct State {
    // All the variables needed to run the simulation
    …
}

static void read_parameters(State& state) {
    cout << "Lato cella lungo x\n";
    cin >> state.Lx;
    …
}

static void initialize_state(State& state) {
    …
    state.cells.resize(state.Nx * state.Ny);
    …
    for (auto& cell: state.cells) {
        cell.phi1 = cell.ftilde;
        cell.phi0 = cell.ftilde;
        cell.phi2 = 0;
    }
    …
}

static void run_simulation(State& state) {
    …
}

static void process_results(State& state) {
    …
}

static void write_results(State &state) {
    std::ofstream fcharge("carica.dat");
    std::ofstream fphi("phi.dat");
    std::ofstream ffield("campo.dat");
    …
}

int main() {
    State state;
    read_parameters(state);
    initialize_state(state);
    run_simulation(state);
    process_results(state);
    write_results(state);
}

Some of the new functios will also be very large, and you could split them up further. A good rule of thumb is that functions should not be larger than one screen. Also, in the example, all the state used in the program is in a struct State and passed to every function. However, not every function needs access to everything, so struct State can probably be split up into multiple structs, so you limit what information is available to each function. You will probably have to iterate this process a few times to get to a situation where you have many small, concise functions and types that together form a cohesive program.

Another benefit of creating functions and structs is that these will have names that describe what they do. This helps document the code, and reduces the need to add comments.