Improved silver-muller conditions

[VadimRomansky]

The problem and the context

The standard SM conditions supposed that electric field should be zero, otherwise they create some unphysical oscillation even in simple case of constant field.

Current Sm conditions are:

( *B[1] )( iB_[1], j )
= Alpha_ * ( *E[2] )( iB_[0] , j )
+ Beta_ * ( ( *B[1] )( iB_[1]-sign_, j )-B_val[1][j] )
+ Gamma_ * b1[j]
+ Delta_ * ( ( B[0] )( iB_[0], j+1 )-B_val[0][j+1] )
+ Epsilon_ ( ( *B[0] )( iB_[0], j )-B_val[0][j ] )
+ B_val[1][j];

The first term in this equation creates problems if E = const != 0. So I suggest to use E_val similary to E_val and rewrite conditions as
( *B[1] )( iB_[1], j )
= Alpha_ * (( *E[2] )( iB_[0] , j ) - E_val[2])
+ Beta_ * ( ( *B[1] )( iB_[1]-sign_, j )-B_val[1][j] )
+ Gamma_ * b1[j]
+ Delta_ * ( ( B[0] )( iB_[0], j+1 )-B_val[0][j+1] )
+ Epsilon_ ( ( *B[0] )( iB_[0], j )-B_val[0][j ] )
+ B_val[1][j];
where only difference from constant level will matter

Failed attempts

I tried to introduced it myself and found some strange things like
void ElectroMagnBC2D_SM::save_fields( Field *my_field, Patch *patch )
{
Field2D *field2D=static_cast<Field2D *>( my_field );

if( patch->isBoundary( i_boundary_ ) ) {
    
    unsigned int xyz = 0;
    if( field2D->name=="Bx" ) {
        xyz = 0;
    } else if( field2D->name=="By" ) {
        xyz = 1;
    } else if( field2D->name=="Bz" ) {
        xyz = 2;
    }
    
    if( axis0_ == 0 ) {
        for( unsigned int j=0; j<B_val[xyz].size(); j++ ) {
            B_val[xyz][j]=( *field2D )( iB_[xyz], j );
        }
    } else {
        for( unsigned int i=0; i<B_val[xyz].size(); i++ ) {
            B_val[xyz][i]=( *field2D )( i, iB_[xyz] );
        }
    }.....

What happens when we get here with E field? xyz = 0 and B_val is incorrect?

[mccoys]

Note: surround code parts with ```

[mccoys]

You could make a new variable E_val and another variable similar to xyz

[VadimRomansky]

Hi all! It seemes that i managed to do it - now Ey is constant, but after finishing work i have exception in destructor.
Can anyone check - have i done it right
ElectroMagnBC2D_SM.cpp.txt
ElectroMagnBC2D_SM.h.txt
?

[mccoys]

I suggest you fork the Smilei repository in your own account, then modify the files there, so that it is easier to follow your modifications.

[mccoys]

Please also show the error you get.

Note: you probably have something wrong in the definition of your variable iE_ as the E field arrays do not have the same shape as B fields. Anyways you should not need a new variable, as iB_ also serves for the E field in the apply functions. For instance Alpha_ * ( *E[2] )( iB_[0] , j )

[VadimRomansky]

I suggest you fork the Smilei repository in your own account, then modify the files there, so that it is easier to follow your modifications.

already done

exception in destructor magically disappeared
If i use iB instead of iE i have problems in line

"
if (axis0_ == 0) {
for (unsigned int j = 0; j < E_val[xyz].size(); j++) {
E_val[xyz][j] = (*field2D)(iB_[xyz], j);
//E_val[xyz][j] = 0;
}
} else {
for (unsigned int i = 0; i < E_val[xyz].size(); i++) {
E_val[xyz][i] = (*field2D)(i, iB_[xyz]);
//E_val[xyz][i] = 0;
}
}
"

 Applying external fields at time t = 0

Stack trace (most recent call last):
#6 Object "[0xffffffffffffffff]", at 0xffffffffffffffff, in
#5 Object "./build/smilei", at 0x55f08a669dbd, in _start
#4 Object "/lib/x86_64-linux-gnu/libc.so.6", at 0x7f284ef9f0b2, in __libc_start_main
#3 Object "./build/smilei", at 0x55f08a922c64, in main
#2 Object "./build/smilei", at 0x55f08a8528eb, in VectorPatch::saveExternalFields(Params&)
#1 Object "./build/smilei", at 0x55f08a708978, in ElectroMagn::saveExternalFields(Patch*)
#0 Object "./build/smilei", at 0x55f08a77ac05, in ElectroMagnBC2D_SM::save_fields(Field*, Patch*)
Segmentation fault (Address not mapped to object [0x31])
Segmentation fault (core dumped)

i don't fully understand how grid is arranged and which indexes i should use

[mccoys]

The Yee grid is staggered according to this: https://smileipic.github.io/Smilei/algorithms.html#time-and-space-discretization

Basically, we speak about primal and dual directions. Dual directions have 1 more point than the primal directions because the first point of the dual grid is at -1/2 compared to the first point of the primal grid at 0.

Each field can be primal in some directions, and dual in other directions.
Ex is dual, primal, primal
Ey is primal, dual, primal
Ez is primal, primal, dual
Bx is primal, dual, dual
By is dual, primal, dual
Bz is dual, dual, primal

But for something simpler, look at those lines:

if( axis0_ == 0 ) {
    Alpha_  *   ( *E[2] )( iB_[0]      , j )
} else {
    Alpha_ *   ( *E[2] )( j, iB_[1]        )
}
if( axis0_ == 0 ) {
    Alpha_ *   ( *E[1] )( iB_[0]      , j )
} else {
    Alpha_ *   ( *E[0] )( j, iB_[1]       )
}

You can see that what matters is the value of axis_. If ==0, then use iB_[0], otherwise use iB_[1].

[VadimRomansky]

Thanks a lot! will fix it now

[VadimRomansky]

did as you said, looks like works