open TELEMAC-MASCARET The mathematically superior suite of solvers

Dear Jean-Michel,

Thank you very much for helping me with fortran coding for tracers value at closed wall. I have run it in Serial and it run. When I looked at results at those points, then I am not sure if ATABOL/BTABOL values were applied. I think I need to reduce them to 0.1.
Just to check if I understood you correctly: for ADR(1...)=2D boundary point number. So for ADR() I used the 2D boundary points, but not the global numbers, is it right?

If it is correct, then looks like it is working with hardcoded boundary point numbers in Serial. Could you please give me some guides on how to make it working in Parallel, when you get a minute. Thank you very much. I really appreciate your help.

Kind Regards!

Violeta

Hello Violeta,

Yes, your use of boundary points numbers is correct. Here are some explanations and bits of code for parallelism.

Basically you now need to find what is the array ADR on every processor, keeping in mind that not all the 6 points will be on a given processor.
If a point is not in a processor, its ADR is 0.

You will have to rely on array KNOLG, because there is no equivalent for boundary nodes, so you need to know the global number of your 6 points (taken in the boundary conditions ASCII file):

274:399
275:391
276:392
277:396
278:398
279:402

Now a tentative algorithm :

Initialising ADR to 0.

Finding ADR:

DO IPTFR=1,NPTFR2
IF(MESH2D%KNOLG%I(NBOR2%I(IPTFR)).EQ.399) ADR(1)=IPTFR
IF(MESH2D%KNOLG%I(NBOR2%I(IPTFR)).EQ.391) ADR(2)=IPTFR
IF(MESH2D%KNOLG%I(NBOR2%I(IPTFR)).EQ.392) ADR(3)=IPTFR
IF(MESH2D%KNOLG%I(NBOR2%I(IPTFR)).EQ.396) ADR(4)=IPTFR
IF(MESH2D%KNOLG%I(NBOR2%I(IPTFR)).EQ.398) ADR(5)=IPTFR
IF(MESH2D%KNOLG%I(NBOR2%I(IPTFR)).EQ.402) ADR(6)=IPTFR
ENDDO

then your previous loop on II, but with a test IF(ADR(II).GT.0) THEN...

DO II=1,6
IF(ADR(II).GT.0) THEN
IPTFR=ADR(II)
....
....
ENDIF
ENDDO


All this is to be put in the following structure:

IF(NCSIZE.GT.1) THEN
the new algorithm
ELSE
your previous scalar implementation
ENDIF

I hope this is not too wrong, not tested!

With best regards,

Jean-Michel Hervouet

Thank you so much Jean-Michel. I could not do it myself. I will use the code and see how it goes. Thank you very much. Sure after many tries it will work, as always:).
I will keep you posted.Thanks.

Kind Regards!

Violeta

Dear Jean-Michel,

Thank you for helping me with this case. After many trials, the code has started to work in parallel. But I am not sure if i get the right results. I used values of ATABOL and BTABOL as 0.1 on all layers, just to see if it is working. Eventually I wish to achieve 0 degree Temperature, and 0 psu for Salinity on those nodes on all layers, constant in time. So basically to represent it as a glacier wall but without any discharge. Here is the latest code from bord3d. I have attached the fortran file and pic of vertical section on those nodes. The tracer values at those nodes are same as everywhere in the domain. But good news is that it is not crashing in Parallel. Is there anything I can do to make it working? Thanks a mil for your help and time.

Kind Regards!

Violeta


IF(NCSIZE.GT.1)THEN
DO IPTFR=1,NPTFR2
IF(MESH2D%KNOLG%I(NBOR2%I(IPTFR)).EQ.399)ADR(1)=IPTFR
IF(MESH2D%KNOLG%I(NBOR2%I(IPTFR)).EQ.391)ADR(1)=IPTFR
IF(MESH2D%KNOLG%I(NBOR2%I(IPTFR)).EQ.392)ADR(1)=IPTFR
IF(MESH2D%KNOLG%I(NBOR2%I(IPTFR)).EQ.396)ADR(1)=IPTFR
IF(MESH2D%KNOLG%I(NBOR2%I(IPTFR)).EQ.398)ADR(1)=IPTFR
IF(MESH2D%KNOLG%I(NBOR2%I(IPTFR)).EQ.402)ADR(1)=IPTFR
ENDDO
IF(ADR(II).GT.0)THEN
IPTFR=ADR(II)
DO II=1,6
IPTFR = ADR(II)
IPOIN2 = NBOR2%I(IPTFR)
IF(NTRAC.NE.0)THEN
DO ITRAC=1,NTRAC
DO IPLAN = 1,NPLAN
IBORD = (IPLAN-1)*NPTFR2 + IPTFR
ATABOL%ADR(ITRAC)%P%R((IPLAN-1)*NPTFR2+IPTFR) =0.1D0
BTABOL%ADR(ITRAC)%P%R((IPLAN-1)*NPTFR2+IPTFR) =0.1D0
ENDDO
ENDDO
ENDIF
ENDDO
ENDIF
ELSE
ADR(1)=274
ADR(2)=275
ADR(3)=276
ADR(4)=277
ADR(5)=278
ADR(6)=279
DO II=1,6
IPTFR = ADR(II)
IPOIN2 = NBOR2%I(IPTFR)
IF(NTRAC.NE.0)THEN
DO ITRAC=1,NTRAC
DO IPLAN = 1,NPLAN
IBORD = (IPLAN-1)*NPTFR2 + IPTFR
ATABOL%ADR(ITRAC)%P%R((IPLAN-1)*NPTFR2+IPTFR) =0.1D0
BTABOL%ADR(ITRAC)%P%R((IPLAN-1)*NPTFR2+IPTFR) =0.1D0
ENDDO
ENDDO
ENDIF
ENDDO
ENDIF
ENDIF

Hello,

I see several problems:

In the parallel part, there should be before all:
DO II=1,6
ADR(II)=0
ENDDO

and then the sequence:

IF(MESH2D%KNOLG%I(NBOR2%I(IPTFR)).EQ.399)ADR(1)=IPTFR
IF(MESH2D%KNOLG%I(NBOR2%I(IPTFR)).EQ.391)ADR(1)=IPTFR
IF(MESH2D%KNOLG%I(NBOR2%I(IPTFR)).EQ.392)ADR(1)=IPTFR
IF(MESH2D%KNOLG%I(NBOR2%I(IPTFR)).EQ.396)ADR(1)=IPTFR
IF(MESH2D%KNOLG%I(NBOR2%I(IPTFR)).EQ.398)ADR(1)=IPTFR
IF(MESH2D%KNOLG%I(NBOR2%I(IPTFR)).EQ.402)ADR(1)=IPTFR

should be:

IF(MESH2D%KNOLG%I(NBOR2%I(IPTFR)).EQ.399)ADR(1)=IPTFR
IF(MESH2D%KNOLG%I(NBOR2%I(IPTFR)).EQ.391)ADR(2)=IPTFR
IF(MESH2D%KNOLG%I(NBOR2%I(IPTFR)).EQ.392)ADR(3)=IPTFR
IF(MESH2D%KNOLG%I(NBOR2%I(IPTFR)).EQ.396)ADR(4)=IPTFR
IF(MESH2D%KNOLG%I(NBOR2%I(IPTFR)).EQ.398)ADR(5)=IPTFR
IF(MESH2D%KNOLG%I(NBOR2%I(IPTFR)).EQ.402)ADR(6)=IPTFR

then:

IF(ADR(II).GT.0)THEN
IPTFR=ADR(II)
DO II=1,6
IPTFR = ADR(II)
IPOIN2 = NBOR2%I(IPTFR)
IF(NTRAC.NE.0)THEN
DO ITRAC=1,NTRAC
DO IPLAN = 1,NPLAN
IBORD = (IPLAN-1)*NPTFR2 + IPTFR
ATABOL%ADR(ITRAC)%P%R((IPLAN-1)*NPTFR2+IPTFR) =0.1D0
BTABOL%ADR(ITRAC)%P%R((IPLAN-1)*NPTFR2+IPTFR) =0.1D0
ENDDO
ENDDO
ENDIF
ENDDO
ENDIF

should be:

DO II=1,6
IF(ADR(II).GT.0)THEN
IPTFR = ADR(II)
IF(NTRAC.NE.0)THEN
DO ITRAC=1,NTRAC
DO IPLAN = 1,NPLAN
IBORD = (IPLAN-1)*NPTFR2 + IPTFR
ATABOL%ADR(ITRAC)%P%R(IBORD) =0.1D0
BTABOL%ADR(ITRAC)%P%R(IBORD) =0.1D0
ENDDO
ENDDO
ENDIF
ENDIF
ENDDO

Then I do not understand why ATABOL = BTABOL

If we want a flux law at boundary in the form:

NU * dT/dn=a*(T-T0) it gives :

ATABOL = a
BTABOL = -a*T0

and if I am not wrong with n pointing to exterior of domain the value a must be negative.

I hope this will improve the results,

Regards,

JMH

Dear Jean-Michel,

Thank you so much for your help. It is running now. I will see results tomorrow. Fingers crossed. Thanks again.

Kind Regards!

Violeta

Dear Jean-Michel,

Thank you for your patience. The new code works fine, except I do not see the desired Temperature at the points I picked. I have tried different values of ATABOL and BTABOL, but without success. Maybe it is that I do not understand how these coefficients work on the vertical planes. At those points I have different values of Sal/Temp on every horizontal level. So it means that each plane must have different values for ATABOL and BTABOL, isn't it? And then it has to be different for Salinity and Temperature as well? Sorry to be a pain. But I think we are close now to get the right results:) Thank you.

Kind Regards!

Violeta

Hello Violeta,

I think I see the point that I overlooked so far, in subroutine LIMI3D, ATABOL and BTABOL are coded as 0 for optimisation, so in bord3d.f when you fill them you have to say that they are no longer 0, by adding, for example just before the RETURN command:

IF(NTRAC.NE.0)THEN
DO ITRAC =1,NTRAC
ATABOL%ADR(ITRAC)%P%TYPR='Q'
ATABOL%ADR(ITRAC)%P%TYPR='Q'
ENDDO
ENDIF

Now you should see an effect. The arrays act on the lateral boundaries. As I explained in my previous post they allow to give a flux across the lateral boundaries, expressing the derivative of the tracer along a direction normal to the boundary as a function of an exterior tracer that I call T0, this flux depends on the turbulent diffusion NU. The variational formulation of the diffusion terms in the tracer equation, after an integration by parts, give an integral along the boundary of NU*dT/dn, and you have to give the value of NU*DT/dn in the form ATABOL*T + BTABOL, this is done so just to split into implicit and explicit terms, which are treated differently. The physics of fluxes through a solid boundary tells you that this flux is proportional to the difference between the value of T on either sides of the boundary, hence something like a*(T-T0), a being negative if the vector n points outward.

Fingers crossed again,

JMH

Dear Jean-Michel,

Thank you so much for code and for detailed explanations. It is much appreciated. It is running now, so hopefully to produce desired results. Thanks again. I will keep you posted. Fingers crossed:)

Kind Regards!

Violeta

Dear Jean-Michel,

Thanks again for your explanations in the previous post. The addition of below

IF(NTRAC.NE.0)THEN
DO ITRAC =1,NTRAC
ATABOL%ADR(ITRAC)%P%TYPR='Q'
ATABOL%ADR(ITRAC)%P%TYPR='Q'
ENDDO
ENDIF

did not stop the simulation but it produced the

GRACJG (BIEF) : EXCEEDING MAXIMUM ITERATIONS: 500 RELATIVE PRECISION: NaN

error right at the beginning. I came across this error before few times and there always were different fixes for that. Last time I had to use the modified by yourself the streamline.f subroutine and it worked fine since for all my cases. I do need to use methods of characteristics as my stratification is not stable with other. Maybe there is something else needs to be done in streamline.f so model can read these boundary nodes in parallel. Here is the post I am referring to. I have attached the error file and the latest fortran file which contain the modified streamline subroutine. Thank you for your help, as always.

www.openmascaret.org/index.php/assistanc...500-relative-precisi

Kind Regards!

Violeta